Fork me on GitHub

source: git/modules/JetFakeParticle.cc@ c749957

ImprovedOutputFile Timing llp
Last change on this file since c749957 was 410f3a2, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

fixed charge variable in JetFakeParticle module

  • Property mode set to 100644
File size: 5.6 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class JetFakeParticle
21 *
22 * Converts jet into particle with some PID,
23 * according to parametrized probability.
24 *
25 * \author M. Selvaggi - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/JetFakeParticle.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootClassifier.h"
38
39#include "TMath.h"
40#include "TString.h"
41#include "TFormula.h"
42#include "TRandom3.h"
43#include "TObjArray.h"
44#include "TDatabasePDG.h"
45#include "TLorentzVector.h"
46
47#include <algorithm>
48#include <stdexcept>
49#include <iostream>
50#include <sstream>
51
52
53using namespace std;
54
55//------------------------------------------------------------------------------
56
57JetFakeParticle::JetFakeParticle() :
58 fItInputArray(0)
59{
60}
61
62//------------------------------------------------------------------------------
63
64JetFakeParticle::~JetFakeParticle()
65{
66}
67
68//------------------------------------------------------------------------------
69
70void JetFakeParticle::Init()
71{
72 TFakeMap::iterator itEfficiencyMap;
73 ExRootConfParam param;
74 DelphesFormula *formula;
75 Int_t i, size, pdgCode;
76
77 // read efficiency formulas
78 param = GetParam("EfficiencyFormula");
79 size = param.GetSize();
80
81 fEfficiencyMap.clear();
82
83 for(i = 0; i < size/2; ++i)
84 {
85 formula = new DelphesFormula;
86 formula->Compile(param[i*2 + 1].GetString());
87 pdgCode = param[i*2].GetInt();
88
89 if(TMath::Abs(pdgCode) != 11 && TMath::Abs(pdgCode) != 13 && TMath::Abs(pdgCode) != 22)
90 {
91 throw runtime_error("Jets can only fake into electrons, muons or photons. Other particles are not authorized.");
92 }
93
94 fEfficiencyMap[param[i*2].GetInt()] = formula;
95 }
96
97 // set default efficiency formula
98 itEfficiencyMap = fEfficiencyMap.find(0);
99 if(itEfficiencyMap == fEfficiencyMap.end())
100 {
101 formula = new DelphesFormula;
102 formula->Compile("0.0");
103
104 fEfficiencyMap[0] = formula;
105 }
106
107 // import input array
108
109 fInputArray = ImportArray(GetString("InputArray", "FastJetFinder/jets"));
110 fItInputArray = fInputArray->MakeIterator();
111
112 // create output array
113
114 fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "fakeElectrons"));
115 fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "fakeMuons"));
116 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "fakePhotons"));
117 fJetOutputArray = ExportArray(GetString("JetOutputArray", "jets"));
118
119}
120
121//------------------------------------------------------------------------------
122
123void JetFakeParticle::Finish()
124{
125 if(fItInputArray) delete fItInputArray;
126
127 TFakeMap::iterator itEfficiencyMap;
128 DelphesFormula *formula;
129 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
130 {
131 formula = itEfficiencyMap->second;
132 if(formula) delete formula;
133 }
134}
135
136//------------------------------------------------------------------------------
137
138void JetFakeParticle::Process()
139{
140 Candidate *candidate, *fake = 0;
141 Double_t pt, eta, phi, e;
142 TFakeMap::iterator itEfficiencyMap;
143 DelphesFormula *formula;
144 Int_t pdgCodeOut;
145
146 Double_t p, r, rs, total;
147
148 fItInputArray->Reset();
149 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
150 {
151 const TLorentzVector &candidatePosition = candidate->Position;
152 const TLorentzVector &candidateMomentum = candidate->Momentum;
153 eta = candidatePosition.Eta();
154 phi = candidatePosition.Phi();
155 pt = candidateMomentum.Pt();
156 e = candidateMomentum.E();
157
158 r = gRandom->Uniform();
159 total = 0.0;
160 fake = 0;
161
162 // loop over map for this jet
163 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
164 {
165 formula = itEfficiencyMap->second;
166 pdgCodeOut = itEfficiencyMap->first;
167
168 p = formula->Eval(pt, eta, phi, e);
169
170 if(total <= r && r < total + p)
171 {
172 fake = static_cast<Candidate*>(candidate->Clone());
173
174 // convert jet
175
176 if(TMath::Abs(pdgCodeOut) == 11 || TMath::Abs(pdgCodeOut) == 13)
177 {
178 if(candidate->Charge != 0)
179 {
180 fake->Charge = candidate->Charge/TMath::Abs(candidate->Charge);
181 }
182 else
183 {
184 rs = gRandom->Uniform();
185 fake->Charge = (rs < 0.5) ? -1 : 1;
186
187 }
188 }
189
190 if(TMath::Abs(pdgCodeOut) == 22) fake->PID = 22;
191
192 if(TMath::Abs(pdgCodeOut) == 11) fElectronOutputArray->Add(fake);
193 if(TMath::Abs(pdgCodeOut) == 13) fMuonOutputArray->Add(fake);
194 if(TMath::Abs(pdgCodeOut) == 22) fPhotonOutputArray->Add(fake);
195
196 break;
197 }
198
199 total += p;
200 }
201
202 if(!fake) fJetOutputArray->Add(candidate);
203 }
204}
205
206//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.