Fork me on GitHub

source: git/modules/JetFlavourAssociation.cc@ edf10ba

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

add JetFlavourAssociation
  • Property mode set to 100644
File size: 12.7 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 JetFlavourAssociation
21 *
22 * Find origin of jet && evaluate jet flavor
23 *
24 * \author P. Demin - UCL, Louvain-la-Neuve
25 *
26 */
27
28#include "modules/JetFlavourAssociation.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 PartonClassifierLHEF: public ExRootClassifier
90{
91public:
92
93 PartonClassifierLHEF() {}
94 Int_t GetCategory(TObject *object);
95 Double_t fEtaMax, fPTMin;
96};
97
98Int_t PartonClassifierLHEF::GetCategory(TObject *object)
99{
100 // select parton in the parton list
101
102 Candidate *parton = static_cast<Candidate *>(object);
103 const TLorentzVector &momentum = parton->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(parton->PID);
110 if(parton->Status == -1) return -1;
111 if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
112 if(parton->Status != 1) return -1; // if status 3 return
113
114 return 0;
115}
116
117//------------------------------------------------------------------------------
118
119JetFlavourAssociation::JetFlavourAssociation() :
120 fClassifier(0), fFilter(0),
121 fItPartonInputArray(0), fItPartonInputArrayLHEF(0),
122 fItJetInputArray(0), fItParticleInputArray(0)
123{
124 fClassifier = new PartonClassifier;
125 fClassifierLHEF = new PartonClassifierLHEF;
126}
127
128//------------------------------------------------------------------------------
129
130JetFlavourAssociation::~JetFlavourAssociation()
131{
132 if(fClassifier) delete fClassifier;
133 if(fClassifierLHEF) delete fClassifierLHEF;
134}
135
136//------------------------------------------------------------------------------
137
138void JetFlavourAssociation::Init()
139{
140 ExRootConfParam param;
141
142 fDeltaR = GetDouble("DeltaR", 0.5);
143
144 fClassifier->fPTMin = GetDouble("PartonPTMin", 0.);
145 fClassifier->fEtaMax = GetDouble("PartonEtaMax",2.5);
146
147 fClassifierLHEF->fPTMin = GetDouble("PartonPTMin", 0.);
148 fClassifierLHEF->fEtaMax = GetDouble("PartonEtaMax",2.5);
149
150 // import input array(s)
151 fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
152 fItPartonInputArray = fPartonInputArray->MakeIterator();
153
154 fPartonInputArrayLHEF = ImportArray(GetString("PartonInputArrayLHEF", "Delphes/partonsLHEF"));
155 fItPartonInputArrayLHEF = fPartonInputArrayLHEF->MakeIterator();
156
157 fFilter = new ExRootFilter(fPartonInputArray);
158 fFilterLHEF = new ExRootFilter(fPartonInputArrayLHEF);
159
160 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
161 fItJetInputArray = fJetInputArray->MakeIterator();
162
163 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
164 fItParticleInputArray = fParticleInputArray->MakeIterator();
165
166}
167
168//------------------------------------------------------------------------------
169
170void JetFlavourAssociation::Finish()
171{
172 if(fFilter) delete fFilter;
173 if(fFilterLHEF) delete fFilterLHEF;
174
175 if(fItJetInputArray) delete fItJetInputArray;
176 if(fItParticleInputArray) delete fItParticleInputArray;
177 if(fItPartonInputArray) delete fItPartonInputArray;
178 if(fItPartonInputArrayLHEF) delete fItPartonInputArrayLHEF;
179}
180
181//------------------------------------------------------------------------------
182
183void JetFlavourAssociation::Process(){
184
185 Candidate *jet;
186 TObjArray *partonArray;
187 TObjArray *partonArrayLHEF;
188
189 // select quark && gluons
190 fFilter->Reset();
191 partonArray = fFilter->GetSubArray(fClassifier, 0); // get the filtered parton array
192
193 if(partonArray == 0) return;
194 TIter itPartonArray(partonArray);
195
196 fFilterLHEF->Reset();
197 partonArrayLHEF = fFilterLHEF->GetSubArray(fClassifierLHEF, 0); // get the filtered parton array
198
199 if(partonArrayLHEF == 0) return;
200 TIter itPartonLHEFArray(partonArrayLHEF);
201
202 // loop over all input jets
203 fItJetInputArray->Reset();
204 while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
205 {
206 // get standard flavor
207 GetAlgoFlavour(jet, itPartonArray, itPartonLHEFArray);
208 GetPhysicsFlavour(jet, itPartonArray, itPartonLHEFArray);
209 }
210}
211
212//------------------------------------------------------------------------------
213// Standard definition of jet flavor in
214// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/JetPartonMatcher.cc?v=CMSSW_7_3_0_pre1
215
216void JetFlavourAssociation::GetAlgoFlavour(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
217{
218 float maxPt = 0;
219 float minDr = 1000;
220 bool isGoodParton = true;
221 int daughterCounter = 0;
222 Candidate *parton, *partonLHEF;
223 Candidate *tempParton = 0, *tempPartonNearest = 0, *tempPartonHighestPt = 0;
224 int pdgCode, pdgCodeMax = -1;
225
226 itPartonArray.Reset();
227 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
228 {
229 // default delphes method
230 pdgCode = TMath::Abs(parton->PID);
231 if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
232 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
233 {
234 if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
235 }
236
237 isGoodParton = true;
238 itPartonLHEFArray.Reset();
239 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
240 {
241 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 &&
242 parton->PID == partonLHEF->PID &&
243 partonLHEF->Charge == parton->Charge)
244 {
245 isGoodParton = false;
246 break;
247 }
248
249 if(!isGoodParton) continue;
250
251 // check the daugheter
252 daughterCounter = 0;
253 if(parton->D1 != -1 || parton->D2 != -1)
254 {
255 // partons are only quarks || gluons
256 int daughterFlavour1 = -1;
257 int daughterFlavour2 = -1;
258 if(parton->D1 != -1) daughterFlavour1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
259 if(parton->D2 != -1) daughterFlavour2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
260 if((daughterFlavour1 == 1 || daughterFlavour1 == 2 || daughterFlavour1 == 3 || daughterFlavour1 == 4 || daughterFlavour1 == 5 || daughterFlavour1 == 21)) daughterCounter++;
261 if((daughterFlavour2 == 1 || daughterFlavour2 == 2 || daughterFlavour2 == 3 || daughterFlavour2 == 4 || daughterFlavour1 == 5 || daughterFlavour2 == 21)) daughterCounter++;
262 }
263 if(daughterCounter > 0) continue;
264 if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
265 {
266 if(jet->Momentum.DeltaR(parton->Momentum) < minDr)
267 {
268 minDr = jet->Momentum.DeltaR(parton->Momentum);
269 tempPartonNearest = parton;
270 }
271
272 // if not yet found && pdgId is a c, take as c
273 if(TMath::Abs(parton->PID) == 4) tempParton = parton;
274 if(TMath::Abs(parton->PID) == 5) tempParton = parton;
275 if(parton->Momentum.Pt() > maxPt)
276 {
277 maxPt = parton->Momentum.Pt();
278 tempPartonHighestPt = parton;
279 }
280 }
281 }
282 }
283
284 jet->FlavorHeaviest = tempParton ? TMath::Abs(tempParton->PID) : 0;
285 jet->FlavorHighestPt = tempPartonHighestPt ? TMath::Abs(tempPartonHighestPt->PID) : 0;
286 jet->FlavorNearest2 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
287 if(!tempParton) tempParton = tempPartonHighestPt;
288 jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
289
290 if(pdgCodeMax == 0) pdgCodeMax = 21;
291 if(pdgCodeMax == -1) pdgCodeMax = 0;
292
293 jet->FlavorDefault = pdgCodeMax;
294}
295
296//------------------------------------------------------------------------------
297
298void JetFlavourAssociation::GetPhysicsFlavour(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
299{
300 float minDr = 1000;
301 int partonCounter = 0;
302 float biggerConeSize = 0.7;
303 float dist;
304 bool isGoodCandidate;
305 int contaminatingFlavour = 0;
306 int motherCounter = 0;
307 Candidate *parton, *partonLHEF, *mother1, *mother2;
308 Candidate *tempParton = 0, *tempPartonNearest = 0;
309 vector<Candidate *> contaminations;
310 vector<Candidate *>::iterator itContaminations;
311
312 contaminations.clear();
313
314 itPartonLHEFArray.Reset();
315 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
316 {
317 dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
318 if(partonLHEF->Status == 1 && dist < minDr)
319 {
320 tempPartonNearest = partonLHEF;
321 minDr = dist;
322 }
323
324 if(partonLHEF->Status == 1 && dist <= fDeltaR)
325 {
326 tempParton = partonLHEF;
327 partonCounter++;
328 }
329 }
330
331 itPartonArray.Reset();
332 itPartonLHEFArray.Reset();
333 while((parton = static_cast<Candidate *>(itPartonArray.Next())))
334 {
335 dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
336 isGoodCandidate = true;
337 while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
338 {
339 if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 &&
340 parton->PID == partonLHEF->PID &&
341 partonLHEF->Charge == parton->Charge)
342 {
343 isGoodCandidate = false;
344 break;
345 }
346 }
347
348 if(!isGoodCandidate) continue;
349
350 if(parton->D1 != -1 || parton->D2 != -1)
351 {
352 if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
353 if(dist < biggerConeSize) contaminations.push_back(parton);
354 }
355 }
356
357 jet->FlavorNearest3 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
358
359 if(partonCounter != 1)
360 {
361 jet->FlavorPhysics = 0;
362 }
363 else if(contaminations.size() == 0)
364 {
365 jet->FlavorPhysics = TMath::Abs(tempParton->PID);
366 }
367 else if(contaminations.size() > 0)
368 {
369 jet->FlavorPhysics = TMath::Abs(tempParton->PID);
370
371 for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
372 {
373 parton = *itContaminations;
374 contaminatingFlavour = TMath::Abs(parton->PID);
375 motherCounter = 0;
376 if(parton->M1 != -1) motherCounter++;
377 if(parton->M2 != -1) motherCounter++;
378
379 if(parton->M1 != -1)
380 {
381 mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
382 if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
383 }
384 if(parton->M2 != -1)
385 {
386 mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
387 if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
388 }
389 // mother is the initialParton --> OK
390 if(TMath::Abs(tempParton->PID) == 4)
391 {
392 // keep association --> the initialParton is a c --> the contaminated parton is a c
393 if(contaminatingFlavour == 4) continue;
394 jet->FlavorPhysics = 0; // all the other cases reject!
395 break;
396 }
397 }
398 }
399}
Note: See TracBrowser for help on using the repository browser.