Fork me on GitHub

source: git/modules/JetFlavorAssociation.cc@ 5d2481f

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

merge BTagging and BTaggingCMS
  • Property mode set to 100644
File size: 12.2 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 bool isGoodParton = true;
219 int daughterCounter = 0;
220 Candidate *parton, *partonLHEF;
221 Candidate *tempParton = 0, *tempPartonHighestPt = 0;
222 int pdgCode, pdgCodeMax = -1;
223
224 itPartonArray.Reset();
225 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
226 {
227 // default delphes method
228 pdgCode = TMath::Abs(parton->PID);
229 if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
230 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
231 {
232 if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
233 }
234
235 isGoodParton = true;
236 itPartonLHEFArray.Reset();
237 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
238 {
239 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 &&
240 parton->PID == partonLHEF->PID &&
241 partonLHEF->Charge == parton->Charge)
242 {
243 isGoodParton = false;
244 break;
245 }
246
247 if(!isGoodParton) continue;
248
249 // check the daugheter
250 daughterCounter = 0;
251 if(parton->D1 != -1 || parton->D2 != -1)
252 {
253 // partons are only quarks || gluons
254 int daughterFlavor1 = -1;
255 int daughterFlavor2 = -1;
256 if(parton->D1 != -1) daughterFlavor1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
257 if(parton->D2 != -1) daughterFlavor2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
258 if((daughterFlavor1 == 1 || daughterFlavor1 == 2 || daughterFlavor1 == 3 || daughterFlavor1 == 4 || daughterFlavor1 == 5 || daughterFlavor1 == 21)) daughterCounter++;
259 if((daughterFlavor2 == 1 || daughterFlavor2 == 2 || daughterFlavor2 == 3 || daughterFlavor2 == 4 || daughterFlavor1 == 5 || daughterFlavor2 == 21)) daughterCounter++;
260 }
261 if(daughterCounter > 0) continue;
262 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
263 {
264 // if not yet found && pdgId is a c, take as c
265 if(TMath::Abs(parton->PID) == 4) tempParton = parton;
266 if(TMath::Abs(parton->PID) == 5) tempParton = parton;
267 if(parton->Momentum.Pt() > maxPt)
268 {
269 maxPt = parton->Momentum.Pt();
270 tempPartonHighestPt = parton;
271 }
272 }
273 }
274 }
275
276 if(!tempParton) tempParton = tempPartonHighestPt;
277 jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
278
279 if(pdgCodeMax == 0) pdgCodeMax = 21;
280 if(pdgCodeMax == -1) pdgCodeMax = 0;
281
282 jet->Flavor = pdgCodeMax;
283}
284
285//------------------------------------------------------------------------------
286
287void JetFlavorAssociation::GetPhysicsFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
288{
289 int partonCounter = 0;
290 float biggerConeSize = 0.7;
291 float dist;
292 bool isGoodCandidate;
293 int contaminatingFlavor = 0;
294 int motherCounter = 0;
295 Candidate *parton, *partonLHEF, *mother1, *mother2;
296 Candidate *tempParton = 0;
297 vector<Candidate *> contaminations;
298 vector<Candidate *>::iterator itContaminations;
299
300 contaminations.clear();
301
302 itPartonLHEFArray.Reset();
303 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
304 {
305 dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
306
307 if(partonLHEF->Status == 1 && dist <= fDeltaR)
308 {
309 tempParton = partonLHEF;
310 partonCounter++;
311 }
312 }
313
314 itPartonArray.Reset();
315 itPartonLHEFArray.Reset();
316 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
317 {
318 dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
319 isGoodCandidate = true;
320 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
321 {
322 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 &&
323 parton->PID == partonLHEF->PID &&
324 partonLHEF->Charge == parton->Charge)
325 {
326 isGoodCandidate = false;
327 break;
328 }
329 }
330
331 if(!isGoodCandidate) continue;
332
333 if(parton->D1 != -1 || parton->D2 != -1)
334 {
335 if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
336 if(dist < biggerConeSize) contaminations.push_back(parton);
337 }
338 }
339
340 if(partonCounter != 1)
341 {
342 jet->FlavorPhys = 0;
343 }
344 else if(contaminations.size() == 0)
345 {
346 jet->FlavorPhys = TMath::Abs(tempParton->PID);
347 }
348 else if(contaminations.size() > 0)
349 {
350 jet->FlavorPhys = TMath::Abs(tempParton->PID);
351
352 for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
353 {
354 parton = *itContaminations;
355 contaminatingFlavor = TMath::Abs(parton->PID);
356 motherCounter = 0;
357 if(parton->M1 != -1) motherCounter++;
358 if(parton->M2 != -1) motherCounter++;
359
360 if(parton->M1 != -1)
361 {
362 mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
363 if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
364 }
365 if(parton->M2 != -1)
366 {
367 mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
368 if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
369 }
370 // mother is the initialParton --> OK
371 if(TMath::Abs(tempParton->PID) == 4)
372 {
373 // keep association --> the initialParton is a c --> the contaminated parton is a c
374 if(contaminatingFlavor == 4) continue;
375 jet->FlavorPhys = 0; // all the other cases reject!
376 break;
377 }
378 }
379 }
380}
Note: See TracBrowser for help on using the repository browser.