Fork me on GitHub

source: git/modules/JetFlavorAssociation.cc@ 1d9974e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 1d9974e was 1d9974e, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

fix default input arrays
  • Property mode set to 100644
File size: 12.9 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
154 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
155 fItParticleInputArray = fParticleInputArray->MakeIterator();
156
157 fParticleLHEFInputArray = ImportArray(GetString("ParticleLHEFInputArray", "Delphes/allParticlesLHEF"));
158 fItParticleLHEFInputArray = fParticleLHEFInputArray->MakeIterator();
159
160 fPartonFilter = new ExRootFilter(fPartonInputArray);
161 fParticleLHEFFilter = new ExRootFilter(fParticleLHEFInputArray);
162
163 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
164 fItJetInputArray = fJetInputArray->MakeIterator();
165}
166
167//------------------------------------------------------------------------------
168
169void JetFlavorAssociation::Finish()
170{
171 if(fPartonFilter) delete fPartonFilter;
172 if(fParticleLHEFFilter) delete fParticleLHEFFilter;
173
174 if(fItJetInputArray) delete fItJetInputArray;
175 if(fItParticleLHEFInputArray) delete fItParticleLHEFInputArray;
176 if(fItParticleInputArray) delete fItParticleInputArray;
177 if(fItPartonInputArray) delete fItPartonInputArray;
178}
179
180//------------------------------------------------------------------------------
181
182void JetFlavorAssociation::Process(){
183
184 Candidate *jet;
185 TObjArray *partonArray;
186 TObjArray *partonLHEFArray;
187
188 // select quark && gluons
189 fPartonFilter->Reset();
190 partonArray = fPartonFilter->GetSubArray(fPartonClassifier, 0); // get the filtered parton array
191
192 if(partonArray == 0) return;
193 TIter itPartonArray(partonArray);
194
195 fParticleLHEFFilter->Reset();
196 partonLHEFArray = fParticleLHEFFilter->GetSubArray(fParticleLHEFClassifier, 0); // get the filtered parton array
197
198 if(partonLHEFArray == 0) return;
199 TIter itPartonLHEFArray(partonLHEFArray);
200
201 // loop over all input jets
202 fItJetInputArray->Reset();
203 while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
204 {
205 // get standard flavor
206 GetAlgoFlavor(jet, itPartonArray, itPartonLHEFArray);
207 GetPhysicsFlavor(jet, itPartonArray, itPartonLHEFArray);
208 }
209}
210
211//------------------------------------------------------------------------------
212// Standard definition of jet flavor in
213// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/JetPartonMatcher.cc?v=CMSSW_7_3_0_pre1
214
215void JetFlavorAssociation::GetAlgoFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
216{
217 float maxPt = 0;
218 float minDr = 1000;
219 bool isGoodParton = true;
220 int daughterCounter = 0;
221 Candidate *parton, *partonLHEF;
222 Candidate *tempParton = 0, *tempPartonNearest = 0, *tempPartonHighestPt = 0;
223 int pdgCode, pdgCodeMax = -1;
224
225 itPartonArray.Reset();
226 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
227 {
228 // default delphes method
229 pdgCode = TMath::Abs(parton->PID);
230 if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
231 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
232 {
233 if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
234 }
235
236 isGoodParton = true;
237 itPartonLHEFArray.Reset();
238 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
239 {
240 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 &&
241 parton->PID == partonLHEF->PID &&
242 partonLHEF->Charge == parton->Charge)
243 {
244 isGoodParton = false;
245 break;
246 }
247
248 if(!isGoodParton) continue;
249
250 // check the daugheter
251 daughterCounter = 0;
252 if(parton->D1 != -1 || parton->D2 != -1)
253 {
254 // partons are only quarks || gluons
255 int daughterFlavor1 = -1;
256 int daughterFlavor2 = -1;
257 if(parton->D1 != -1) daughterFlavor1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
258 if(parton->D2 != -1) daughterFlavor2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
259 if((daughterFlavor1 == 1 || daughterFlavor1 == 2 || daughterFlavor1 == 3 || daughterFlavor1 == 4 || daughterFlavor1 == 5 || daughterFlavor1 == 21)) daughterCounter++;
260 if((daughterFlavor2 == 1 || daughterFlavor2 == 2 || daughterFlavor2 == 3 || daughterFlavor2 == 4 || daughterFlavor1 == 5 || daughterFlavor2 == 21)) daughterCounter++;
261 }
262 if(daughterCounter > 0) continue;
263 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
264 {
265 if(jet->Momentum.DeltaR(parton->Momentum) < minDr)
266 {
267 minDr = jet->Momentum.DeltaR(parton->Momentum);
268 tempPartonNearest = parton;
269 }
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 jet->FlavorHeaviest = tempParton ? TMath::Abs(tempParton->PID) : 0;
284 jet->FlavorHighestPt = tempPartonHighestPt ? TMath::Abs(tempPartonHighestPt->PID) : 0;
285 jet->FlavorNearest2 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
286 if(!tempParton) tempParton = tempPartonHighestPt;
287 jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
288
289 if(pdgCodeMax == 0) pdgCodeMax = 21;
290 if(pdgCodeMax == -1) pdgCodeMax = 0;
291
292 jet->FlavorDefault = pdgCodeMax;
293}
294
295//------------------------------------------------------------------------------
296
297void JetFlavorAssociation::GetPhysicsFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
298{
299 float minDr = 1000;
300 int partonCounter = 0;
301 float biggerConeSize = 0.7;
302 float dist;
303 bool isGoodCandidate;
304 int contaminatingFlavor = 0;
305 int motherCounter = 0;
306 Candidate *parton, *partonLHEF, *mother1, *mother2;
307 Candidate *tempParton = 0, *tempPartonNearest = 0;
308 vector<Candidate *> contaminations;
309 vector<Candidate *>::iterator itContaminations;
310
311 contaminations.clear();
312
313 itPartonLHEFArray.Reset();
314 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
315 {
316 dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
317 if(partonLHEF->Status == 1 && dist < minDr)
318 {
319 tempPartonNearest = partonLHEF;
320 minDr = dist;
321 }
322
323 if(partonLHEF->Status == 1 && dist <= fDeltaR)
324 {
325 tempParton = partonLHEF;
326 partonCounter++;
327 }
328 }
329
330 itPartonArray.Reset();
331 itPartonLHEFArray.Reset();
332 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
333 {
334 dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
335 isGoodCandidate = true;
336 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
337 {
338 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 &&
339 parton->PID == partonLHEF->PID &&
340 partonLHEF->Charge == parton->Charge)
341 {
342 isGoodCandidate = false;
343 break;
344 }
345 }
346
347 if(!isGoodCandidate) continue;
348
349 if(parton->D1 != -1 || parton->D2 != -1)
350 {
351 if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
352 if(dist < biggerConeSize) contaminations.push_back(parton);
353 }
354 }
355
356 jet->FlavorNearest3 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
357
358 if(partonCounter != 1)
359 {
360 jet->FlavorPhysics = 0;
361 }
362 else if(contaminations.size() == 0)
363 {
364 jet->FlavorPhysics = TMath::Abs(tempParton->PID);
365 }
366 else if(contaminations.size() > 0)
367 {
368 jet->FlavorPhysics = TMath::Abs(tempParton->PID);
369
370 for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
371 {
372 parton = *itContaminations;
373 contaminatingFlavor = TMath::Abs(parton->PID);
374 motherCounter = 0;
375 if(parton->M1 != -1) motherCounter++;
376 if(parton->M2 != -1) motherCounter++;
377
378 if(parton->M1 != -1)
379 {
380 mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
381 if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
382 }
383 if(parton->M2 != -1)
384 {
385 mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
386 if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
387 }
388 // mother is the initialParton --> OK
389 if(TMath::Abs(tempParton->PID) == 4)
390 {
391 // keep association --> the initialParton is a c --> the contaminated parton is a c
392 if(contaminatingFlavor == 4) continue;
393 jet->FlavorPhysics = 0; // all the other cases reject!
394 break;
395 }
396 }
397 }
398}
Note: See TracBrowser for help on using the repository browser.