ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/GeneratorTools.cc
Revision: 1.1
Committed: Sat Mar 27 13:38:44 2010 UTC (15 years, 1 month ago) by sixie
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a
Log Message: 
Add Generator Tools

File Contents

# User Rev Content
1 sixie 1.1 // $Id $
2    
3     #include "MitPhysics/Utils/interface/GeneratorTools.h"
4     #include "MitCommon/MathTools/interface/MathUtils.h"
5     #include "MitAna/DataCont/interface/ObjArray.h"
6     #include "MitAna/DataUtil/interface/Debug.h"
7    
8     using namespace mithep;
9     using namespace std;
10    
11     //--------------------------------------------------------------------------------------------------
12     void GeneratorTools::PrintHepMCTable(const mithep::Collection<mithep::MCParticle> *particles,
13     Bool_t showDaughters = true, int suppressEntriesAfterThisIndex = -1) {
14    
15     //Print entire HepMC Decay Table
16    
17     for (UInt_t d=0; d<particles->GetEntries(); d++) {
18     if (suppressEntriesAfterThisIndex == -1 || int(d) <= suppressEntriesAfterThisIndex) {
19     cout << "Particle " << d << " : " << particles->At(d)->PdgId() << " "
20     << particles->At(d)->Status() << " "
21     << particles->At(d)->Pt() << " "
22     << particles->At(d)->Eta() << " " << particles->At(d)->Phi()
23     << " IsParton:" << particles->At(d)->IsParton()
24     << " Decays at " << particles->At(d)->DecayVertex().Rho() << " ";
25     if (particles->At(d)->HasMother()) {
26     cout << " Mother: " << particles->At(d)->Mother()->PdgId() << " "
27     << particles->At(d)->Mother()->Status() << " "
28     << particles->At(d)->Mother()->Pt() << " "
29     << particles->At(d)->Mother()->Eta() << " "
30     << particles->At(d)->Mother()->Phi() << " ";
31     }
32     cout << " " << particles->At(d)->NDaughters() << " Daughters : ";
33     if (showDaughters) {
34     for (UInt_t e=0;e<particles->At(d)->NDaughters();e++) {
35     if (e <= 4 || particles->At(d)->Daughter(e)->Pt() > 5.0 ) {
36     cout << " | "
37     << particles->At(d)->Daughter(e)->PdgId() << " "
38     << particles->At(d)->Daughter(e)->Status() << " "
39     << particles->At(d)->Daughter(e)->Pt() << " "
40     << particles->At(d)->Daughter(e)->Eta() << " "
41     << particles->At(d)->Daughter(e)->Phi() << " "
42     << " ";
43     } else {
44     cout << " ... ";
45     }
46     }
47     }
48     cout << endl;
49     }
50     }
51     return;
52     }
53    
54     //--------------------------------------------------------------------------------------------------
55     void GeneratorTools::PrintNearbyParticles(const mithep::Collection<mithep::MCParticle> *particles,
56     Double_t eta, Double_t phi, Double_t deltaR = 0.3) {
57    
58     //Print Particles in the HepMC table which are within the given DeltaR distance
59    
60     mithep::ObjArray<mithep::MCParticle> *descendants = new mithep::ObjArray<mithep::MCParticle>;
61    
62     //Print all particles and their descendants that are within the dR Cone.
63     for (UInt_t d=0; d<particles->GetEntries(); d++) {
64     Double_t dR = MathUtils::DeltaR(phi,eta,particles->At(d)->Phi(),particles->At(d)->Eta());
65    
66     //print out the particle if it is inside the dR cone, or if it is a descendant
67     if (dR < deltaR || descendants->HasObject(particles->At(d))) {
68     cout << "Particle " << d << " : " << particles->At(d)->PdgId() << " "
69     << particles->At(d)->Status() << " "
70     << particles->At(d)->Pt() << " "
71     << particles->At(d)->Eta() << " " << particles->At(d)->Phi()
72     << " IsParton:" << particles->At(d)->IsParton()
73     << " Decays at " << particles->At(d)->DecayVertex().Rho()
74     << " " << particles->At(d)->NDaughters() << " Daughters : ";
75    
76     for (UInt_t e=0;e<particles->At(d)->NDaughters();e++) {
77     descendants->Add(particles->At(d)->Daughter(e));
78     if (e <= 4 || particles->At(d)->Daughter(e)->Pt() > 5.0) {
79     cout << " | "
80     << particles->At(d)->Daughter(e)->PdgId() << " "
81     << particles->At(d)->Daughter(e)->Status() << " "
82     << particles->At(d)->Daughter(e)->Pt() << " "
83     << particles->At(d)->Daughter(e)->Eta() << " "
84     << particles->At(d)->Daughter(e)->Phi() << " "
85     << " ";
86     }
87     }
88     cout << " ... ";
89     cout << endl;
90     }
91     }
92     return;
93     }
94    
95     //***********************************************************************************************
96     //Match Electron to a SimTrack
97     //***********************************************************************************************
98     const MCParticle* GeneratorTools::MatchElectronToSimParticle(const mithep::Collection<mithep::MCParticle> *particles,
99     const mithep::Track *eleTrack, Bool_t isGsfTrack, Int_t printDebugLevel,
100     Bool_t printHepMCTable) {
101    
102     //Match Clean Electron to a SimTrack and then trace back up the mother chain.
103     const MCParticle *match = NULL;
104    
105     Double_t matchDeltaPtOverPt = 0.0;
106     Double_t secondaryMatchDeltaPtOverPt = 0.0;
107     Double_t minProductionVertexRho = 10000.0;
108     Double_t secondaryMatchMinDR = 5000.0;
109     Double_t minDR = 5000.0;
110     const MCParticle *matchedSimChargedParticle = NULL;
111     const MCParticle *secondaryMatchedSimChargedParticle = NULL;
112    
113     if (printDebugLevel >= 5) cout << "\nSimulated particles near Electron\n";
114     for(UInt_t j=0;j<particles->GetEntries();j++) {
115     Bool_t isStable = false;
116     Bool_t isTrackable = false;
117     if(particles->At(j)->NDaughters() == 0 || particles->At(j)->DecayVertex().Rho() > 10.0)
118     isStable = true;
119    
120     if (particles->At(j)->Mother() && particles->At(j)->Mother()->DecayVertex().Rho() <= 40)
121     isTrackable = true;
122    
123     if (particles->At(j)->IsSimulated()) {
124     Double_t DR = MathUtils::DeltaR(eleTrack->Phi(), eleTrack->Eta(), particles->At(j)->Phi(),
125     particles->At(j)->Eta());
126    
127     Double_t dPtOverPt = (eleTrack->Pt() - particles->At(j)->Pt())/eleTrack->Pt();
128    
129     //remove Neutrals
130     if (particles->At(j)->Charge() != 0) {
131     if( DR < 0.3) {
132     if (printDebugLevel >= 5) {
133     cout << "ChargedSimulatedParticle: " << j << " : " << fabs(dPtOverPt) << " : " << DR
134     << " " << particles->At(j)->PdgId()
135     << " " << particles->At(j)->Status() << " " << particles->At(j)->Pt()
136     << " " << particles->At(j)->Eta() << " " << particles->At(j)->Phi()
137     << " produced at ";
138    
139     if (particles->At(j)->Mother())
140     cout << particles->At(j)->Mother()->DecayVertex().Rho();
141     else
142     cout << "0.0";
143     cout << " decays at ";
144     if (particles->At(j)->NDaughters() > 0)
145     cout << particles->At(j)->DecayVertex().Rho();
146     else
147     cout << " Stable ";
148     cout << endl;
149     }
150     }
151    
152     //regular match
153     if (DR < minDR && fabs(dPtOverPt) < 0.5 && DR < 0.3 && isStable && isTrackable ) {
154     minDR = DR;
155     matchedSimChargedParticle = particles->At(j);
156     }
157    
158     Double_t productionVertexRho = 20000.0;
159     if (particles->At(j)->Mother())
160     productionVertexRho = particles->At(j)->Mother()->DecayVertex().Rho();
161    
162     //perform secondary match for GSF track only electrons
163     if ( DR < 0.3 && particles->At(j)->AbsPdgId() == 11 && particles->At(j)->Pt() > 3.0
164     && productionVertexRho < minProductionVertexRho ) {
165     minProductionVertexRho = productionVertexRho;
166     secondaryMatchedSimChargedParticle = particles->At(j);
167     }
168     } //is charged particle
169     } //is simulated
170     }// for all particles
171    
172     if (printDebugLevel >= 5) {
173     if (matchedSimChargedParticle && minDR < 0.3) {
174     cout << "Match : " << matchDeltaPtOverPt << " : " << minDR
175     << " " << matchedSimChargedParticle->PdgId()
176     << " " << matchedSimChargedParticle->Status() << " " << matchedSimChargedParticle->Pt()
177     << " " << matchedSimChargedParticle->Eta() << " " << matchedSimChargedParticle->Phi()
178     << " produced at ";
179     if (matchedSimChargedParticle->Mother())
180     cout << matchedSimChargedParticle->Mother()->DecayVertex().Rho();
181     else
182     cout << "0.0";
183     cout << " decays at ";
184     if (matchedSimChargedParticle->NDaughters() > 0)
185     cout << matchedSimChargedParticle->DecayVertex().Rho();
186     else
187     cout << " Stable ";
188     cout << endl;
189     }
190     }
191    
192     if (!(matchedSimChargedParticle && minDR < 0.3)) {
193     if (secondaryMatchedSimChargedParticle && isGsfTrack) {
194     matchedSimChargedParticle = secondaryMatchedSimChargedParticle;
195     matchDeltaPtOverPt = secondaryMatchDeltaPtOverPt;
196     minDR = secondaryMatchMinDR;
197     }
198    
199     if (printDebugLevel >= 5) {
200     if (matchedSimChargedParticle && minDR < 0.3) {
201     cout << "secondaryMatch : " << matchDeltaPtOverPt << " : " << minDR
202     << " " << matchedSimChargedParticle->PdgId()
203     << " " << matchedSimChargedParticle->Status() << " " << matchedSimChargedParticle->Pt()
204     << " " << matchedSimChargedParticle->Eta() << " " << matchedSimChargedParticle->Phi()
205     << " produced at ";
206     if (matchedSimChargedParticle->Mother())
207     cout << matchedSimChargedParticle->Mother()->DecayVertex().Rho();
208     else
209     cout << "0.0";
210     cout << " decays at ";
211     if (matchedSimChargedParticle->NDaughters() > 0)
212     cout << matchedSimChargedParticle->DecayVertex().Rho();
213     else
214     cout << " Stable ";
215     cout << endl;
216     }
217     }
218     }
219    
220     if (matchedSimChargedParticle && minDR < 0.3)
221     match = matchedSimChargedParticle;
222    
223     //MDB(kAnalysis, 6) PrintNearbyParticles(particles, ele->Eta(), ele->Phi(), 0.3);
224     MDB(kAnalysis, 7) if (printHepMCTable) PrintHepMCTable(particles);
225    
226     return match;
227     }
228    
229     //***********************************************************************************************
230     //Return the meson ancestor given the input matched simulated particle. Intended to be the
231     //particle corresponding to the large x fragmentation resulting in an electron fake.
232     //***********************************************************************************************
233     const MCParticle* GeneratorTools::FindElectronFakeAncestor(const mithep::MCParticle *matchedSimChargedParticle) {
234    
235     if (!matchedSimChargedParticle)
236     return 0;
237    
238     //If it's an electron, try
239     //1) to match to a conversion, then find ancestor pi0
240     //2) find B / D meson if it's semi-leptonic decay
241     //3) if it's a prompt photon, return null.
242     if (matchedSimChargedParticle->AbsPdgId() == 11) {
243    
244     //if no mother then return null
245     if (!matchedSimChargedParticle->Mother()) {
246     MDB(kAnalysis, 5) cout << "no mother\n";
247     return 0;
248     }
249    
250    
251     //while look to look for photon mother
252     Bool_t FoundPhotonMother = false;
253     const MCParticle *tmpParticle = matchedSimChargedParticle;
254     while (tmpParticle->Mother() && tmpParticle->Mother()->PdgId() != 22 ) {
255     tmpParticle = tmpParticle->Mother();
256     }
257     if (tmpParticle->Mother() && tmpParticle->Mother()->PdgId() == 22) {
258     FoundPhotonMother = true;
259     }
260    
261     //conversions
262     if (FoundPhotonMother) {
263     //find the meson ancestor.
264     const MCParticle *tmpParticle = matchedSimChargedParticle;
265     while (tmpParticle->Mother()
266     &&
267     tmpParticle->Mother()->PdgId() != 111
268     &&
269     tmpParticle->Mother()->PdgId() != 221
270     &&
271     tmpParticle->Mother()->PdgId() != 223
272     &&
273     !(tmpParticle->Mother()->PdgId() == 22 && tmpParticle->Mother()->Mother() &&
274     tmpParticle->Mother()->Mother()->IsParton())
275     &&
276     tmpParticle->Mother()->AbsPdgId() != 24
277     &&
278     tmpParticle->Mother()->AbsPdgId() != 23
279     &&
280     tmpParticle->Mother()->AbsPdgId() != 92
281     &&
282     !tmpParticle->Mother()->IsParton()
283     ) {
284     MDB(kAnalysis, 5) {
285     cout << "find : " << tmpParticle->Mother()->PdgId() << " " << tmpParticle->Mother()->Status()
286     << " " << tmpParticle->Mother()->Pt() << " " << tmpParticle->Mother()->Eta() << " "
287     << tmpParticle->Mother()->Phi() << endl;
288     }
289     tmpParticle = tmpParticle->Mother();
290     }
291     if (tmpParticle->Mother()->PdgId() == 111
292     ||tmpParticle->Mother()->PdgId() == 221
293     ||tmpParticle->Mother()->PdgId() == 223
294     ) {
295     return tmpParticle->Mother();
296     } else if (tmpParticle->Mother()->PdgId() == 22
297     || tmpParticle->Mother()->PdgId() == 23
298     || tmpParticle->Mother()->PdgId() == 24
299     ) {
300     //prompt photon or FSR photon or Brem photon
301     return 0;
302     } else {
303     return 0;
304     }
305     //end if sim electron mother was a photon
306     } else {
307     //now look for semi-leptonic decays
308     const MCParticle *tmpParticle = matchedSimChargedParticle;
309     while (tmpParticle->Mother()
310     &&
311     !((tmpParticle->Mother()->AbsPdgId() >= 400 &&
312     tmpParticle->Mother()->AbsPdgId() < 600) ||
313     (tmpParticle->Mother()->AbsPdgId() >= 4000 &&
314     tmpParticle->Mother()->AbsPdgId() < 6000)
315     )
316     ) {
317     MDB(kAnalysis, 5) {
318     cout << "find : " << tmpParticle->Mother()->PdgId() << " " << tmpParticle->Mother()->Status()
319     << " " << tmpParticle->Mother()->Pt() << " " << tmpParticle->Mother()->Eta() << " "
320     << tmpParticle->Mother()->Phi() << endl;
321     }
322     tmpParticle = tmpParticle->Mother();
323     }
324     if (tmpParticle->Mother() && tmpParticle->Mother()->HasDaughter(11)) {
325     return tmpParticle->Mother();
326     }
327     }
328     //end if sim particle was an electron
329     } else if (matchedSimChargedParticle->AbsPdgId() == 211 ) {
330     return matchedSimChargedParticle;
331     } else if (matchedSimChargedParticle->AbsPdgId() == 321) {
332     return matchedSimChargedParticle;
333     } else {
334     return 0;
335     }
336     return 0;
337     }
338    
339    
340     //***********************************************************************************************
341     //Categorizes Electron Fakes based on the input matched simulated particle
342     //***********************************************************************************************
343     Int_t GeneratorTools::CategorizeFakeElectron(const mithep::MCParticle *matchedSimChargedParticle) {
344    
345     //Categorizes Electron Fakes based on the input matched simulated particle
346     Int_t FakeCategory = 0;
347    
348     if (!matchedSimChargedParticle)
349     return FakeCategory;
350    
351     //Try to match to conversion
352     if (matchedSimChargedParticle->AbsPdgId() == 11) {
353     const MCParticle *mother = matchedSimChargedParticle->Mother();
354    
355     if (!mother) {
356     MDB(kAnalysis, 5) cout << "no mother\n";
357     FakeCategory = 7;
358     } else {
359     //try to find an ancestor that was a pi0 or a heavy meson
360     MDB(kAnalysis, 5) cout << "find pi0 or heavy meson ancestor\n";
361     const MCParticle *tmpParticle = matchedSimChargedParticle;
362     while (tmpParticle->Mother()
363     &&
364     tmpParticle->Mother()->PdgId() != 111
365     &&
366     tmpParticle->Mother()->PdgId() != 221
367     &&
368     tmpParticle->Mother()->PdgId() != 223
369     &&
370     !(tmpParticle->Mother()->PdgId() == 22 && tmpParticle->Mother()->Mother() &&
371     tmpParticle->Mother()->Mother()->IsParton())
372     &&
373     !((tmpParticle->Mother()->AbsPdgId() >= 400 &&
374     tmpParticle->Mother()->AbsPdgId() < 600) ||
375     (tmpParticle->Mother()->AbsPdgId() >= 4000 &&
376     tmpParticle->Mother()->AbsPdgId() < 6000)
377     )
378     &&
379     tmpParticle->Mother()->AbsPdgId() != 24
380     &&
381     tmpParticle->Mother()->AbsPdgId() != 23
382     &&
383     tmpParticle->Mother()->AbsPdgId() != 92
384     &&
385     !tmpParticle->Mother()->IsParton()
386     ) {
387     MDB(kAnalysis, 5) {
388     cout << "find : " << tmpParticle->Mother()->PdgId() << " " << tmpParticle->Mother()->Status()
389     << " " << tmpParticle->Mother()->Pt() << " " << tmpParticle->Mother()->Eta() << " "
390     << tmpParticle->Mother()->Phi() << endl;
391     }
392     tmpParticle = tmpParticle->Mother();
393     }
394     const MCParticle *Ancestor = tmpParticle->Mother();
395     if (Ancestor) {
396    
397     MDB(kAnalysis, 5) {
398     cout << "Found ancestor: " << Ancestor->PdgId() << " "
399     << Ancestor->Status() << " " << Ancestor->Pt()
400     << " " << Ancestor->Eta() << " "
401     << Ancestor->Phi() << " HasElectronDaughter: " << Ancestor->HasDaughter(11)
402     << endl;
403     }
404    
405     if (Ancestor->PdgId() == 111 || Ancestor->PdgId() == 221 || Ancestor->PdgId() == 223) {
406     FakeCategory = 1; //pi0->photon->conversion
407     } else if (Ancestor->PdgId() == 22) {
408     FakeCategory = 6; //prompt photon -> conversion
409     } else if ((Ancestor->AbsPdgId() >= 400 &&
410     Ancestor->AbsPdgId() < 600) ||
411     (Ancestor->AbsPdgId() >= 4000 &&
412     Ancestor->AbsPdgId() < 6000) &&
413     Ancestor->HasDaughter(11)) {
414     FakeCategory = 4; //heavy flavor
415     } else if (Ancestor->AbsPdgId() == 24 || Ancestor->AbsPdgId() == 23) {
416     FakeCategory = 5;
417     } else {
418     FakeCategory = 7;
419     }
420     } else {
421     FakeCategory = 7;
422     }
423     }
424     } else if (matchedSimChargedParticle->AbsPdgId() == 211 ) {
425     FakeCategory = 2;
426     } else if (matchedSimChargedParticle->AbsPdgId() == 321) {
427     FakeCategory = 3;
428     } else {
429     FakeCategory = 7;
430     }
431     return FakeCategory;
432     }
433    
434    
435     //***********************************************************************************************
436     //Match Muon to a SimTrack
437     //***********************************************************************************************
438     const MCParticle* GeneratorTools::MatchMuonToSimParticle(
439     const mithep::Collection<mithep::MCParticle> *particles,const mithep::Track *muonTrack,
440     Bool_t isTrackerTrack, Int_t printDebugLevel, Bool_t printHepMCTable) {
441    
442     //Match Muon to a SimTrack
443     const MCParticle *match = NULL;
444    
445     Double_t minDR = 5000.0;
446     Double_t matchDeltaPtOverPt = 0.0;
447     Double_t secondaryMatchMinDR = 5000.0;
448     Double_t secondaryMatchDeltaPtOverPt = 0.0;
449    
450     const MCParticle *matchedSimChargedParticle = NULL;
451     const MCParticle *secondaryMatchedSimChargedParticle = NULL;
452    
453     if (printDebugLevel >= 5) cout << "\nSimulated particles near Electron\n";
454     for(UInt_t j=0;j<particles->GetEntries();j++) {
455     Bool_t isStable = false;
456     Bool_t isTrackable = false;
457    
458     if (particles->At(j)->IsSimulated()) {
459     Double_t DR = MathUtils::DeltaR(muonTrack->Phi(), muonTrack->Eta(), particles->At(j)->Phi(),
460     particles->At(j)->Eta());
461    
462     Double_t dPtOverPt = (muonTrack->Pt() - particles->At(j)->Pt())/muonTrack->Pt();
463    
464     //remove Neutrals
465     if (particles->At(j)->Charge() != 0) {
466     if( DR < 0.3) {
467     if (printDebugLevel >= 5) {
468     cout << "ChargedSimulatedParticle: " << j << " : " << fabs(dPtOverPt) << " : " << DR
469     << " " << particles->At(j)->PdgId()
470     << " " << particles->At(j)->Status() << " " << particles->At(j)->Pt()
471     << " " << particles->At(j)->Eta() << " " << particles->At(j)->Phi()
472     << " produced at ";
473    
474     if (particles->At(j)->Mother())
475     cout << particles->At(j)->Mother()->DecayVertex().Rho();
476     else
477     cout << "0.0";
478     cout << " decays at ";
479     if (particles->At(j)->NDaughters() > 0)
480     cout << particles->At(j)->DecayVertex().Rho();
481     else
482     cout << " Stable ";
483     cout << endl;
484     }
485     }
486    
487     //regular match. match to a muon
488     if (DR < minDR && DR < 0.3 && particles->At(j)->AbsPdgId() == 13
489     && particles->At(j)->Pt() > 5.0
490     && (particles->At(j)->Mother() &&
491     particles->At(j)->Mother()->DecayVertex().Rho() <= 150
492     )
493     ) {
494     minDR = DR;
495     matchedSimChargedParticle = particles->At(j);
496     matchDeltaPtOverPt = dPtOverPt;
497     }
498    
499     //perform secondary match for non muons
500     if (DR < secondaryMatchMinDR && DR < 0.3 && fabs(dPtOverPt) < 0.5 ) {
501     secondaryMatchMinDR = DR;
502     secondaryMatchedSimChargedParticle = particles->At(j);
503     secondaryMatchDeltaPtOverPt = dPtOverPt;
504     }
505    
506     } //is charged particle
507     } //is simulated
508     }// for all particles
509    
510     if (printDebugLevel >= 5) {
511     if (matchedSimChargedParticle && minDR < 0.3) {
512     cout << "Match : " << matchDeltaPtOverPt << " : " << minDR
513     << " " << matchedSimChargedParticle->PdgId()
514     << " " << matchedSimChargedParticle->Status() << " " << matchedSimChargedParticle->Pt()
515     << " " << matchedSimChargedParticle->Eta() << " " << matchedSimChargedParticle->Phi()
516     << " produced at ";
517     if (matchedSimChargedParticle->Mother())
518     cout << matchedSimChargedParticle->Mother()->DecayVertex().Rho();
519     else
520     cout << "0.0";
521     cout << " decays at ";
522     if (matchedSimChargedParticle->NDaughters() > 0)
523     cout << matchedSimChargedParticle->DecayVertex().Rho();
524     else
525     cout << " Stable ";
526     cout << endl;
527     }
528     }
529    
530     if (!matchedSimChargedParticle) {
531     if (secondaryMatchedSimChargedParticle) {
532     matchedSimChargedParticle = secondaryMatchedSimChargedParticle;
533     minDR = secondaryMatchMinDR;
534     matchDeltaPtOverPt = secondaryMatchDeltaPtOverPt;
535     }
536    
537     if (printDebugLevel >= 5) {
538     if (matchedSimChargedParticle && minDR < 0.3) {
539     cout << "secondaryMatch : " << matchDeltaPtOverPt << " : " << minDR
540     << " " << matchedSimChargedParticle->PdgId()
541     << " " << matchedSimChargedParticle->Status() << " " << matchedSimChargedParticle->Pt()
542     << " " << matchedSimChargedParticle->Eta() << " " << matchedSimChargedParticle->Phi()
543     << " produced at ";
544     if (matchedSimChargedParticle->Mother())
545     cout << matchedSimChargedParticle->Mother()->DecayVertex().Rho();
546     else
547     cout << "0.0";
548     cout << " decays at ";
549     if (matchedSimChargedParticle->NDaughters() > 0)
550     cout << matchedSimChargedParticle->DecayVertex().Rho();
551     else
552     cout << " Stable ";
553     cout << endl;
554     }
555     }
556     }
557    
558     if (matchedSimChargedParticle && minDR < 0.3)
559     match = matchedSimChargedParticle;
560    
561     return match;
562     }
563    
564    
565     //***********************************************************************************************
566     //Categorizes Electron Fakes based on the input matched simulated particle
567     //***********************************************************************************************
568     Int_t GeneratorTools::CategorizeFakeMuon(const mithep::MCParticle *matchedSimChargedParticle) {
569    
570     //Categorizes Electron Fakes based on the input matched simulated particle
571     Int_t FakeCategory = 0;
572    
573     if (!matchedSimChargedParticle)
574     return FakeCategory;
575    
576     //look for a muon first
577     if (matchedSimChargedParticle->AbsPdgId() == 13) {
578     const MCParticle *mother = matchedSimChargedParticle->Mother();
579    
580     if (!mother) {
581     MDB(kAnalysis, 5) cout << "no mother\n";
582     FakeCategory = 7;
583     } else {
584     //try to find an ancestor that was a pi+,K+ or a heavy meson
585     const MCParticle *tmpParticle = matchedSimChargedParticle;
586     while (tmpParticle->Mother()
587     &&
588     tmpParticle->Mother()->AbsPdgId() != 24
589     &&
590     tmpParticle->Mother()->AbsPdgId() != 23
591     &&
592     tmpParticle->Mother()->AbsPdgId() != 211
593     &&
594     tmpParticle->Mother()->AbsPdgId() != 321
595     &&
596     tmpParticle->Mother()->AbsPdgId() != 113
597     &&
598     !((tmpParticle->Mother()->AbsPdgId() >= 400 &&
599     tmpParticle->Mother()->AbsPdgId() < 600) ||
600     (tmpParticle->Mother()->AbsPdgId() >= 4000 &&
601     tmpParticle->Mother()->AbsPdgId() < 6000)
602     )
603     ) {
604     MDB(kAnalysis, 5) {
605     cout << "find : " << tmpParticle->Mother()->PdgId() << " " << tmpParticle->Mother()->Status()
606     << " " << tmpParticle->Mother()->Pt() << " " << tmpParticle->Mother()->Eta() << " "
607     << tmpParticle->Mother()->Phi() << endl;
608     }
609     tmpParticle = tmpParticle->Mother();
610     }
611     const MCParticle *Ancestor = tmpParticle->Mother();
612    
613     if (Ancestor) {
614     MDB(kAnalysis, 5) {
615     cout << "Found ancestor: " << Ancestor->PdgId() << " "
616     << Ancestor->Status() << " " << Ancestor->Pt()
617     << " " << Ancestor->Eta() << " "
618     << Ancestor->Phi() << " HasElectronDaughter: " << Ancestor->HasDaughter(13)
619     << endl;
620     }
621    
622     if (Ancestor->AbsPdgId() == 211) {
623     FakeCategory = 1; //pi+ decay in flight
624     } else if (Ancestor->AbsPdgId() == 321) {
625     FakeCategory = 2; //K+ decay in flight
626     } else if ((Ancestor->AbsPdgId() >= 400 &&
627     Ancestor->AbsPdgId() < 600) ||
628     (Ancestor->AbsPdgId() >= 4000 &&
629     Ancestor->AbsPdgId() < 6000) &&
630     Ancestor->HasDaughter(13)) {
631     FakeCategory = 4; //heavy flavor semileptonic decay
632     } else if ( Ancestor->AbsPdgId() == 24 || Ancestor->AbsPdgId() == 23 ) {
633     FakeCategory = 7; //comes from real muon. bad global fit direction
634     } else {
635     FakeCategory = 3; //muon from something else
636     }
637     } else {
638     FakeCategory = -1; //no ancestor found
639     }
640     }
641     } else if (matchedSimChargedParticle->AbsPdgId() == 211
642     || matchedSimChargedParticle->AbsPdgId() == 321
643     || matchedSimChargedParticle->AbsPdgId() == 2212 ) {
644     FakeCategory = 5; //pion punchthrough
645     } else {
646     FakeCategory = 6; //other punchthrough?
647     }
648     return FakeCategory;
649     }
650    
651    
652    
653     //--------------------------------------------------------------------------------------------------
654     TString GeneratorTools::ConvertPdgIdToName(Int_t pdgId) {
655    
656     //Convert the pdgID code into a particle name
657    
658     TString answer = "";
659    
660     switch ( pdgId )
661     {
662     case 11:
663     answer = "e-";
664     break;
665     case -11:
666     answer = "e+";
667     break;
668     case 12:
669     answer = "nu_e";
670     break;
671     case -12:
672     answer = "nu_e bar";
673     break;
674     case 13:
675     answer = "mu-";
676     break;
677     case -13:
678     answer = "mu+";
679     break;
680     case 15:
681     answer = "tau-";
682     break;
683     case -15:
684     answer = "tau+";
685     break;
686     case 21:
687     answer = "gluon";
688     break;
689     case 22:
690     answer = "photon";
691     break;
692     case 23:
693     answer = "Z";
694     break;
695     case 24:
696     answer = "W+";
697     break;
698     case -24:
699     answer = "W-";
700     break;
701     case 25:
702     answer = "Higgs";
703     break;
704     case 111:
705     answer = "Pi0";
706     break;
707     case 113:
708     answer = "Rho0(770)";
709     break;
710     case 115:
711     answer = "a_2^0(1320)";
712     break;
713     case 117:
714     answer = "rho_3^0(1690)";
715     break;
716     case 119:
717     answer = "a_4^0(2040)";
718     break;
719     case 130:
720     answer = "K_L^0";
721     break;
722    
723    
724    
725     case 211:
726     answer = "Pi+";
727     break;
728     case -211:
729     answer = "Pi-";
730     break;
731     case 213:
732     answer = "Rho+";
733     break;
734     case -213:
735     answer = "Rho-";
736     break;
737     case 221:
738     answer = "Eta";
739     break;
740     case 223:
741     answer = "omega(782)";
742     break;
743     case 310:
744     answer = "K_S^0";
745     break;
746     case 311:
747     answer = "K0";
748     break;
749     case 321:
750     answer = "K+";
751     break;
752     case -321:
753     answer = "K-";
754     break;
755     case 313:
756     answer = "K0*(892)";
757     break;
758     case 323:
759     answer = "K+*(892)";
760     break;
761     case -323:
762     answer = "K-*(892)";
763     break;
764     case 315:
765     answer = "K_2^0*(1430)";
766     break;
767     case 325:
768     answer = "K_2^+*(1430)";
769     break;
770     case -325:
771     answer = "K_2^-*(1430)";
772     break;
773     case 317:
774     answer = "K_3^0*(1780)";
775     break;
776     case 327:
777     answer = "K_3^+*(1780)";
778     break;
779     case -327:
780     answer = "K_3^-*(1780)";
781     break;
782     case 319:
783     answer = "K_4^0*(2045)";
784     break;
785     case 329:
786     answer = "K_4^+*(2045)";
787     break;
788     case -329:
789     answer = "K_4^-*(2045)";
790     break;
791     case 331:
792     answer = "Eta'(958)";
793     break;
794     case 333:
795     answer = "phi'(1020)";
796     break;
797    
798    
799     case 411:
800     answer = "D+";
801     break;
802     case -411:
803     answer = "D-";
804     break;
805     case 421:
806     answer = "D0";
807     break;
808     case -421:
809     answer = "D0 bar";
810     break;
811     case 413:
812     answer = "D^+*(2010)";
813     break;
814     case -413:
815     answer = "D^-*(2010)";
816     break;
817     case 423:
818     answer = "D^0*(2010)";
819     break;
820     case -423:
821     answer = "D^0*(2010) bar";
822     break;
823     case 415:
824     answer = "D_2^+*(2460)";
825     break;
826     case -415:
827     answer = "D_2^-*(2460)";
828     break;
829     case 425:
830     answer = "D_2^0*(2460)";
831     break;
832     case -425:
833     answer = "D_2^0*(2460) bar";
834     break;
835     case 431:
836     answer = "D_s^+";
837     break;
838     case -431:
839     answer = "D_s^-";
840     break;
841     case 433:
842     answer = "D_s^+*";
843     break;
844     case -433:
845     answer = "D_s^-*";
846     break;
847     case 435:
848     answer = "D_s2^+*";
849     break;
850     case -435:
851     answer = "D_s2^-*";
852     break;
853    
854     case 511:
855     answer = "B0";
856     break;
857     case -511:
858     answer = "B0 bar";
859     break;
860     case -521:
861     answer = "B-";
862     break;
863     case 521:
864     answer = "B+";
865     break;
866     case 513:
867     answer = "B0*";
868     break;
869     case -513:
870     answer = "B0* bar";
871     break;
872     case -523:
873     answer = "B-*";
874     break;
875     case 523:
876     answer = "B+*";
877     break;
878     case 515:
879     answer = "B_2^0*";
880     break;
881     case -515:
882     answer = "B_2^0* bar";
883     break;
884     case -525:
885     answer = "B_2^-*";
886     break;
887     case 525:
888     answer = "B_2^+*";
889     break;
890     case 531:
891     answer = "B_s^0";
892     break;
893     case -531:
894     answer = "B_s^0 bar";
895     break;
896     case 533:
897     answer = "B_s^0*";
898     break;
899     case -533:
900     answer = "B_s^0* bar";
901     break;
902     case 535:
903     answer = "B_s2^0*";
904     break;
905     case -535:
906     answer = "B_s2^0* bar";
907     break;
908     case 541:
909     answer = "B_c^+";
910     break;
911     case -541:
912     answer = "B_c^-";
913     break;
914     case 543:
915     answer = "B_c^+*";
916     break;
917     case -543:
918     answer = "B_c^-*";
919     break;
920     case 545:
921     answer = "B_c2^+*";
922     break;
923     case -545:
924     answer = "B_c2^-*";
925     break;
926    
927     case 441:
928     answer = "eta_c(1S)";
929     break;
930     case 443:
931     answer = "JPsi(1S)";
932     break;
933     case 445:
934     answer = "chi_c2(1P)";
935     break;
936     case 551:
937     answer = "eta_b(1S)";
938     break;
939     case 553:
940     answer = "Upsilon(1S)";
941     break;
942     case 555:
943     answer = "chi_b2(1P)";
944     break;
945     case 557:
946     answer = "Upsilon(1D)";
947     break;
948    
949     case 2212:
950     answer = "proton";
951     break;
952     case 2112:
953     answer = "neutron";
954     break;
955     case 2224:
956     answer = "Delta++";
957     break;
958     case 2214:
959     answer = "Delta+";
960     break;
961     case 2114:
962     answer = "Delta0";
963     break;
964     case 1114:
965     answer = "Delta-";
966     break;
967     case 3122:
968     answer = "Lambda";
969     break;
970     case 3222:
971     answer = "Sigma+";
972     break;
973     case 3212:
974     answer = "Sigma0";
975     break;
976     case 3112:
977     answer = "Sigma-";
978     break;
979     case 3224:
980     answer = "Sigma+*";
981     break;
982     case 3214:
983     answer = "Sigma0*";
984     break;
985     case 3114:
986     answer = "Sigma-*";
987     break;
988     case 3322:
989     answer = "Xi0";
990     break;
991     case 3312:
992     answer = "Xi-";
993     break;
994     case 3324:
995     answer = "Xi0*";
996     break;
997     case 3314:
998     answer = "Xi-*";
999     break;
1000     case 3334:
1001     answer = "Omega-";
1002     break;
1003    
1004     default:
1005     answer = "";
1006     }
1007     return answer;
1008     }