Fork me on GitHub

source: git/modules/JetFlavorAssociation.cc@ e5ea42e

Last change on this file since e5ea42e was d1942df, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

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