Fork me on GitHub

source: git/modules/JetFlavorAssociation.cc@ b30245b

Last change on this file since b30245b was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files
  • Property mode set to 100644
File size: 12.3 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 || 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/** \class JetFlavorAssociation
20 *
21 * Find origin of jet && evaluate jet flavor
22 *
23 * \author P. Demin - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "modules/JetFlavorAssociation.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesFormula.h"
32
33#include "ExRootAnalysis/ExRootClassifier.h"
34#include "ExRootAnalysis/ExRootFilter.h"
35#include "ExRootAnalysis/ExRootResult.h"
36
37#include "TDatabasePDG.h"
38#include "TFormula.h"
39#include "TLorentzVector.h"
40#include "TMath.h"
41#include "TObjArray.h"
42#include "TRandom3.h"
43#include "TString.h"
44
45#include <algorithm>
46#include <iostream>
47#include <sstream>
48#include <stdexcept>
49
50using namespace std;
51
52//------------------------------------------------------------------------------
53
54class PartonClassifier : public ExRootClassifier
55{
56public:
57 PartonClassifier() {}
58 Int_t GetCategory(TObject *object);
59 Double_t fEtaMax, fPTMin;
60};
61
62//------------------------------------------------------------------------------
63// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/PartonSelector.cc
64
65Int_t PartonClassifier::GetCategory(TObject *object)
66{
67 // select parton in the parton list
68
69 Candidate *parton = static_cast<Candidate *>(object);
70 const TLorentzVector &momentum = parton->Momentum;
71 Int_t pdgCode;
72
73 // inside the eta && momentum range (be a little bit larger that the tracking coverage
74 if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
75
76 pdgCode = TMath::Abs(parton->PID);
77
78 if(parton->Status == -1) return -1;
79 if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
80 if(parton->Status == 3 || parton->Status == 2) return 0; // if status 3 return
81
82 return 0;
83}
84
85//------------------------------------------------------------------------------
86
87class ParticleLHEFClassifier : public ExRootClassifier
88{
89public:
90 ParticleLHEFClassifier() {}
91 Int_t GetCategory(TObject *object);
92 Double_t fEtaMax, fPTMin;
93};
94
95Int_t ParticleLHEFClassifier::GetCategory(TObject *object)
96{
97 // select parton in the parton list
98
99 Candidate *particleLHEF = static_cast<Candidate *>(object);
100 const TLorentzVector &momentum = particleLHEF->Momentum;
101 Int_t pdgCode;
102
103 // inside the eta && momentum range (be a little bit larger that the tracking coverage
104 if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
105
106 pdgCode = TMath::Abs(particleLHEF->PID);
107 if(particleLHEF->Status == -1) return -1;
108 if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
109 if(particleLHEF->Status != 1) return -1; // if status 3 return
110
111 return 0;
112}
113
114//------------------------------------------------------------------------------
115
116JetFlavorAssociation::JetFlavorAssociation() :
117 fPartonClassifier(0), fPartonFilter(0), fParticleLHEFFilter(0),
118 fItPartonInputArray(0), fItParticleInputArray(0),
119 fItParticleLHEFInputArray(0), fItJetInputArray(0)
120{
121 fPartonClassifier = new PartonClassifier;
122 fParticleLHEFClassifier = new ParticleLHEFClassifier;
123}
124
125//------------------------------------------------------------------------------
126
127JetFlavorAssociation::~JetFlavorAssociation()
128{
129 if(fPartonClassifier) delete fPartonClassifier;
130 if(fParticleLHEFClassifier) delete fParticleLHEFClassifier;
131}
132
133//------------------------------------------------------------------------------
134
135void JetFlavorAssociation::Init()
136{
137 ExRootConfParam param;
138
139 fDeltaR = GetDouble("DeltaR", 0.5);
140
141 fPartonClassifier->fPTMin = GetDouble("PartonPTMin", 0.0);
142 fPartonClassifier->fEtaMax = GetDouble("PartonEtaMax", 2.5);
143
144 fParticleLHEFClassifier->fPTMin = GetDouble("PartonPTMin", 0.0);
145 fParticleLHEFClassifier->fEtaMax = GetDouble("PartonEtaMax", 2.5);
146
147 // import input array(s)
148 fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
149 fItPartonInputArray = fPartonInputArray->MakeIterator();
150 fPartonFilter = new ExRootFilter(fPartonInputArray);
151
152 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
153 fItParticleInputArray = fParticleInputArray->MakeIterator();
154
155 try
156 {
157 fParticleLHEFInputArray = ImportArray(GetString("ParticleLHEFInputArray", "Delphes/allParticlesLHEF"));
158 }
159 catch(runtime_error &e)
160 {
161 fParticleLHEFInputArray = 0;
162 }
163
164 if(fParticleLHEFInputArray)
165 {
166 fItParticleLHEFInputArray = fParticleLHEFInputArray->MakeIterator();
167 fParticleLHEFFilter = new ExRootFilter(fParticleLHEFInputArray);
168 }
169
170 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
171 fItJetInputArray = fJetInputArray->MakeIterator();
172}
173
174//------------------------------------------------------------------------------
175
176void JetFlavorAssociation::Finish()
177{
178 if(fPartonFilter) delete fPartonFilter;
179 if(fParticleLHEFFilter) delete fParticleLHEFFilter;
180
181 if(fItJetInputArray) delete fItJetInputArray;
182 if(fItParticleLHEFInputArray) delete fItParticleLHEFInputArray;
183 if(fItParticleInputArray) delete fItParticleInputArray;
184 if(fItPartonInputArray) delete fItPartonInputArray;
185}
186
187//------------------------------------------------------------------------------
188
189void JetFlavorAssociation::Process()
190{
191
192 Candidate *jet;
193 TObjArray *partonArray = 0;
194 TObjArray *partonLHEFArray = 0;
195
196 // select quark and gluons
197 fPartonFilter->Reset();
198 partonArray = fPartonFilter->GetSubArray(fPartonClassifier, 0); // get the filtered parton array
199 if(partonArray == 0) return;
200
201 if(fParticleLHEFInputArray)
202 {
203 fParticleLHEFFilter->Reset();
204 partonLHEFArray = fParticleLHEFFilter->GetSubArray(fParticleLHEFClassifier, 0); // get the filtered parton array
205 }
206 // loop over all input jets
207 fItJetInputArray->Reset();
208 while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
209 {
210 // get standard flavor
211 GetAlgoFlavor(jet, partonArray, partonLHEFArray);
212 if(fParticleLHEFInputArray) GetPhysicsFlavor(jet, partonArray, partonLHEFArray);
213 }
214}
215
216//------------------------------------------------------------------------------
217// Standard definition of jet flavor in
218// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/JetPartonMatcher.cc?v=CMSSW_7_3_0_pre1
219
220void JetFlavorAssociation::GetAlgoFlavor(Candidate *jet, TObjArray *partonArray, TObjArray *partonLHEFArray)
221{
222 float maxPt = 0;
223 int daughterCounter = 0;
224 Candidate *parton, *partonLHEF;
225 Candidate *tempParton = 0, *tempPartonHighestPt = 0;
226 int pdgCode, pdgCodeMax = -1;
227
228 TIter itPartonArray(partonArray);
229 TIter itPartonLHEFArray(partonLHEFArray);
230
231 itPartonArray.Reset();
232 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
233 {
234 // default delphes method
235 pdgCode = TMath::Abs(parton->PID);
236 if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
237 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
238 {
239 if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
240 }
241
242 if(!fParticleLHEFInputArray) continue;
243
244 itPartonLHEFArray.Reset();
245 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
246 {
247 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 && parton->PID == partonLHEF->PID && partonLHEF->Charge == parton->Charge)
248 {
249 break;
250 }
251
252 // check the daughter
253 daughterCounter = 0;
254 if(parton->D1 != -1 || parton->D2 != -1)
255 {
256 // partons are only quarks || gluons
257 int daughterFlavor1 = -1;
258 int daughterFlavor2 = -1;
259 if(parton->D1 != -1) daughterFlavor1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
260 if(parton->D2 != -1) daughterFlavor2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
261 if((daughterFlavor1 == 1 || daughterFlavor1 == 2 || daughterFlavor1 == 3 || daughterFlavor1 == 4 || daughterFlavor1 == 5 || daughterFlavor1 == 21)) daughterCounter++;
262 if((daughterFlavor2 == 1 || daughterFlavor2 == 2 || daughterFlavor2 == 3 || daughterFlavor2 == 4 || daughterFlavor2 == 5 || daughterFlavor2 == 21)) daughterCounter++;
263 }
264 if(daughterCounter > 0) continue;
265 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
266 {
267 // if not yet found && pdgId is a c, take as c
268 if(TMath::Abs(parton->PID) == 4) tempParton = parton;
269 if(TMath::Abs(parton->PID) == 5) tempParton = parton;
270 if(parton->Momentum.Pt() > maxPt)
271 {
272 maxPt = parton->Momentum.Pt();
273 tempPartonHighestPt = parton;
274 }
275 }
276 }
277 }
278
279 if(!tempParton) tempParton = tempPartonHighestPt;
280 jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
281
282 if(pdgCodeMax == 0) pdgCodeMax = 21;
283 if(pdgCodeMax == -1) pdgCodeMax = 0;
284
285 jet->Flavor = pdgCodeMax;
286}
287
288//------------------------------------------------------------------------------
289
290void JetFlavorAssociation::GetPhysicsFlavor(Candidate *jet, TObjArray *partonArray, TObjArray *partonLHEFArray)
291{
292 int partonCounter = 0;
293 float biggerConeSize = 0.7;
294 float dist;
295 bool isGoodCandidate;
296 int contaminatingFlavor = 0;
297 int motherCounter = 0;
298 Candidate *parton, *partonLHEF, *mother1, *mother2;
299 Candidate *tempParton = 0;
300 vector<Candidate *> contaminations;
301 vector<Candidate *>::iterator itContaminations;
302
303 TIter itPartonArray(partonArray);
304 TIter itPartonLHEFArray(partonLHEFArray);
305
306 contaminations.clear();
307
308 itPartonLHEFArray.Reset();
309 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
310 {
311 dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
312
313 if(partonLHEF->Status == 1 && dist <= fDeltaR)
314 {
315 tempParton = partonLHEF;
316 partonCounter++;
317 }
318 }
319
320 itPartonArray.Reset();
321 itPartonLHEFArray.Reset();
322 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
323 {
324 dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
325 isGoodCandidate = true;
326 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
327 {
328 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 && parton->PID == partonLHEF->PID && partonLHEF->Charge == parton->Charge)
329 {
330 isGoodCandidate = false;
331 break;
332 }
333 }
334
335 if(!isGoodCandidate) continue;
336
337 if(parton->D1 != -1 || parton->D2 != -1)
338 {
339 if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
340 if(dist < biggerConeSize) contaminations.push_back(parton);
341 }
342 }
343
344 if(partonCounter != 1)
345 {
346 jet->FlavorPhys = 0;
347 }
348 else if(contaminations.size() == 0)
349 {
350 jet->FlavorPhys = TMath::Abs(tempParton->PID);
351 }
352 else if(contaminations.size() > 0)
353 {
354 jet->FlavorPhys = TMath::Abs(tempParton->PID);
355
356 for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
357 {
358 parton = *itContaminations;
359 contaminatingFlavor = TMath::Abs(parton->PID);
360 motherCounter = 0;
361 if(parton->M1 != -1) motherCounter++;
362 if(parton->M2 != -1) motherCounter++;
363
364 if(parton->M1 != -1)
365 {
366 mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
367 if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
368 }
369 if(parton->M2 != -1)
370 {
371 mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
372 if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
373 }
374 // mother is the initialParton --> OK
375 if(TMath::Abs(tempParton->PID) == 4)
376 {
377 // keep association --> the initialParton is a c --> the contaminated parton is a c
378 if(contaminatingFlavor == 4) continue;
379 jet->FlavorPhys = 0; // all the other cases reject!
380 break;
381 }
382 }
383 }
384}
Note: See TracBrowser for help on using the repository browser.