ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/GeneratorTools.cc
Revision: 1.3
Committed: Fri Jul 8 17:54:27 2011 UTC (13 years, 9 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.2: +3 -2 lines
Log Message: 
*** empty log message ***

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 ceballos 1.2 ((Ancestor->AbsPdgId() >= 4000 &&
412 sixie 1.1 Ancestor->AbsPdgId() < 6000) &&
413 ceballos 1.2 Ancestor->HasDaughter(11))) {
414 sixie 1.1 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 fabstoec 1.3 // removed byb Fabian... they are NOT used
456     // Bool_t isStable = false;
457     // Bool_t isTrackable = false;
458 sixie 1.1
459     if (particles->At(j)->IsSimulated()) {
460     Double_t DR = MathUtils::DeltaR(muonTrack->Phi(), muonTrack->Eta(), particles->At(j)->Phi(),
461     particles->At(j)->Eta());
462    
463     Double_t dPtOverPt = (muonTrack->Pt() - particles->At(j)->Pt())/muonTrack->Pt();
464    
465     //remove Neutrals
466     if (particles->At(j)->Charge() != 0) {
467     if( DR < 0.3) {
468     if (printDebugLevel >= 5) {
469     cout << "ChargedSimulatedParticle: " << j << " : " << fabs(dPtOverPt) << " : " << DR
470     << " " << particles->At(j)->PdgId()
471     << " " << particles->At(j)->Status() << " " << particles->At(j)->Pt()
472     << " " << particles->At(j)->Eta() << " " << particles->At(j)->Phi()
473     << " produced at ";
474    
475     if (particles->At(j)->Mother())
476     cout << particles->At(j)->Mother()->DecayVertex().Rho();
477     else
478     cout << "0.0";
479     cout << " decays at ";
480     if (particles->At(j)->NDaughters() > 0)
481     cout << particles->At(j)->DecayVertex().Rho();
482     else
483     cout << " Stable ";
484     cout << endl;
485     }
486     }
487    
488     //regular match. match to a muon
489     if (DR < minDR && DR < 0.3 && particles->At(j)->AbsPdgId() == 13
490     && particles->At(j)->Pt() > 5.0
491     && (particles->At(j)->Mother() &&
492     particles->At(j)->Mother()->DecayVertex().Rho() <= 150
493     )
494     ) {
495     minDR = DR;
496     matchedSimChargedParticle = particles->At(j);
497     matchDeltaPtOverPt = dPtOverPt;
498     }
499    
500     //perform secondary match for non muons
501     if (DR < secondaryMatchMinDR && DR < 0.3 && fabs(dPtOverPt) < 0.5 ) {
502     secondaryMatchMinDR = DR;
503     secondaryMatchedSimChargedParticle = particles->At(j);
504     secondaryMatchDeltaPtOverPt = dPtOverPt;
505     }
506    
507     } //is charged particle
508     } //is simulated
509     }// for all particles
510    
511     if (printDebugLevel >= 5) {
512     if (matchedSimChargedParticle && minDR < 0.3) {
513     cout << "Match : " << matchDeltaPtOverPt << " : " << minDR
514     << " " << matchedSimChargedParticle->PdgId()
515     << " " << matchedSimChargedParticle->Status() << " " << matchedSimChargedParticle->Pt()
516     << " " << matchedSimChargedParticle->Eta() << " " << matchedSimChargedParticle->Phi()
517     << " produced at ";
518     if (matchedSimChargedParticle->Mother())
519     cout << matchedSimChargedParticle->Mother()->DecayVertex().Rho();
520     else
521     cout << "0.0";
522     cout << " decays at ";
523     if (matchedSimChargedParticle->NDaughters() > 0)
524     cout << matchedSimChargedParticle->DecayVertex().Rho();
525     else
526     cout << " Stable ";
527     cout << endl;
528     }
529     }
530    
531     if (!matchedSimChargedParticle) {
532     if (secondaryMatchedSimChargedParticle) {
533     matchedSimChargedParticle = secondaryMatchedSimChargedParticle;
534     minDR = secondaryMatchMinDR;
535     matchDeltaPtOverPt = secondaryMatchDeltaPtOverPt;
536     }
537    
538     if (printDebugLevel >= 5) {
539     if (matchedSimChargedParticle && minDR < 0.3) {
540     cout << "secondaryMatch : " << matchDeltaPtOverPt << " : " << minDR
541     << " " << matchedSimChargedParticle->PdgId()
542     << " " << matchedSimChargedParticle->Status() << " " << matchedSimChargedParticle->Pt()
543     << " " << matchedSimChargedParticle->Eta() << " " << matchedSimChargedParticle->Phi()
544     << " produced at ";
545     if (matchedSimChargedParticle->Mother())
546     cout << matchedSimChargedParticle->Mother()->DecayVertex().Rho();
547     else
548     cout << "0.0";
549     cout << " decays at ";
550     if (matchedSimChargedParticle->NDaughters() > 0)
551     cout << matchedSimChargedParticle->DecayVertex().Rho();
552     else
553     cout << " Stable ";
554     cout << endl;
555     }
556     }
557     }
558    
559     if (matchedSimChargedParticle && minDR < 0.3)
560     match = matchedSimChargedParticle;
561    
562     return match;
563     }
564    
565    
566     //***********************************************************************************************
567     //Categorizes Electron Fakes based on the input matched simulated particle
568     //***********************************************************************************************
569     Int_t GeneratorTools::CategorizeFakeMuon(const mithep::MCParticle *matchedSimChargedParticle) {
570    
571     //Categorizes Electron Fakes based on the input matched simulated particle
572     Int_t FakeCategory = 0;
573    
574     if (!matchedSimChargedParticle)
575     return FakeCategory;
576    
577     //look for a muon first
578     if (matchedSimChargedParticle->AbsPdgId() == 13) {
579     const MCParticle *mother = matchedSimChargedParticle->Mother();
580    
581     if (!mother) {
582     MDB(kAnalysis, 5) cout << "no mother\n";
583     FakeCategory = 7;
584     } else {
585     //try to find an ancestor that was a pi+,K+ or a heavy meson
586     const MCParticle *tmpParticle = matchedSimChargedParticle;
587     while (tmpParticle->Mother()
588     &&
589     tmpParticle->Mother()->AbsPdgId() != 24
590     &&
591     tmpParticle->Mother()->AbsPdgId() != 23
592     &&
593     tmpParticle->Mother()->AbsPdgId() != 211
594     &&
595     tmpParticle->Mother()->AbsPdgId() != 321
596     &&
597     tmpParticle->Mother()->AbsPdgId() != 113
598     &&
599     !((tmpParticle->Mother()->AbsPdgId() >= 400 &&
600     tmpParticle->Mother()->AbsPdgId() < 600) ||
601     (tmpParticle->Mother()->AbsPdgId() >= 4000 &&
602     tmpParticle->Mother()->AbsPdgId() < 6000)
603     )
604     ) {
605     MDB(kAnalysis, 5) {
606     cout << "find : " << tmpParticle->Mother()->PdgId() << " " << tmpParticle->Mother()->Status()
607     << " " << tmpParticle->Mother()->Pt() << " " << tmpParticle->Mother()->Eta() << " "
608     << tmpParticle->Mother()->Phi() << endl;
609     }
610     tmpParticle = tmpParticle->Mother();
611     }
612     const MCParticle *Ancestor = tmpParticle->Mother();
613    
614     if (Ancestor) {
615     MDB(kAnalysis, 5) {
616     cout << "Found ancestor: " << Ancestor->PdgId() << " "
617     << Ancestor->Status() << " " << Ancestor->Pt()
618     << " " << Ancestor->Eta() << " "
619     << Ancestor->Phi() << " HasElectronDaughter: " << Ancestor->HasDaughter(13)
620     << endl;
621     }
622    
623     if (Ancestor->AbsPdgId() == 211) {
624     FakeCategory = 1; //pi+ decay in flight
625     } else if (Ancestor->AbsPdgId() == 321) {
626     FakeCategory = 2; //K+ decay in flight
627     } else if ((Ancestor->AbsPdgId() >= 400 &&
628     Ancestor->AbsPdgId() < 600) ||
629 ceballos 1.2 ((Ancestor->AbsPdgId() >= 4000 &&
630 sixie 1.1 Ancestor->AbsPdgId() < 6000) &&
631 ceballos 1.2 Ancestor->HasDaughter(13))) {
632 sixie 1.1 FakeCategory = 4; //heavy flavor semileptonic decay
633     } else if ( Ancestor->AbsPdgId() == 24 || Ancestor->AbsPdgId() == 23 ) {
634     FakeCategory = 7; //comes from real muon. bad global fit direction
635     } else {
636     FakeCategory = 3; //muon from something else
637     }
638     } else {
639     FakeCategory = -1; //no ancestor found
640     }
641     }
642     } else if (matchedSimChargedParticle->AbsPdgId() == 211
643     || matchedSimChargedParticle->AbsPdgId() == 321
644     || matchedSimChargedParticle->AbsPdgId() == 2212 ) {
645     FakeCategory = 5; //pion punchthrough
646     } else {
647     FakeCategory = 6; //other punchthrough?
648     }
649     return FakeCategory;
650     }
651    
652    
653    
654     //--------------------------------------------------------------------------------------------------
655     TString GeneratorTools::ConvertPdgIdToName(Int_t pdgId) {
656    
657     //Convert the pdgID code into a particle name
658    
659     TString answer = "";
660    
661     switch ( pdgId )
662     {
663     case 11:
664     answer = "e-";
665     break;
666     case -11:
667     answer = "e+";
668     break;
669     case 12:
670     answer = "nu_e";
671     break;
672     case -12:
673     answer = "nu_e bar";
674     break;
675     case 13:
676     answer = "mu-";
677     break;
678     case -13:
679     answer = "mu+";
680     break;
681     case 15:
682     answer = "tau-";
683     break;
684     case -15:
685     answer = "tau+";
686     break;
687     case 21:
688     answer = "gluon";
689     break;
690     case 22:
691     answer = "photon";
692     break;
693     case 23:
694     answer = "Z";
695     break;
696     case 24:
697     answer = "W+";
698     break;
699     case -24:
700     answer = "W-";
701     break;
702     case 25:
703     answer = "Higgs";
704     break;
705     case 111:
706     answer = "Pi0";
707     break;
708     case 113:
709     answer = "Rho0(770)";
710     break;
711     case 115:
712     answer = "a_2^0(1320)";
713     break;
714     case 117:
715     answer = "rho_3^0(1690)";
716     break;
717     case 119:
718     answer = "a_4^0(2040)";
719     break;
720     case 130:
721     answer = "K_L^0";
722     break;
723    
724    
725    
726     case 211:
727     answer = "Pi+";
728     break;
729     case -211:
730     answer = "Pi-";
731     break;
732     case 213:
733     answer = "Rho+";
734     break;
735     case -213:
736     answer = "Rho-";
737     break;
738     case 221:
739     answer = "Eta";
740     break;
741     case 223:
742     answer = "omega(782)";
743     break;
744     case 310:
745     answer = "K_S^0";
746     break;
747     case 311:
748     answer = "K0";
749     break;
750     case 321:
751     answer = "K+";
752     break;
753     case -321:
754     answer = "K-";
755     break;
756     case 313:
757     answer = "K0*(892)";
758     break;
759     case 323:
760     answer = "K+*(892)";
761     break;
762     case -323:
763     answer = "K-*(892)";
764     break;
765     case 315:
766     answer = "K_2^0*(1430)";
767     break;
768     case 325:
769     answer = "K_2^+*(1430)";
770     break;
771     case -325:
772     answer = "K_2^-*(1430)";
773     break;
774     case 317:
775     answer = "K_3^0*(1780)";
776     break;
777     case 327:
778     answer = "K_3^+*(1780)";
779     break;
780     case -327:
781     answer = "K_3^-*(1780)";
782     break;
783     case 319:
784     answer = "K_4^0*(2045)";
785     break;
786     case 329:
787     answer = "K_4^+*(2045)";
788     break;
789     case -329:
790     answer = "K_4^-*(2045)";
791     break;
792     case 331:
793     answer = "Eta'(958)";
794     break;
795     case 333:
796     answer = "phi'(1020)";
797     break;
798    
799    
800     case 411:
801     answer = "D+";
802     break;
803     case -411:
804     answer = "D-";
805     break;
806     case 421:
807     answer = "D0";
808     break;
809     case -421:
810     answer = "D0 bar";
811     break;
812     case 413:
813     answer = "D^+*(2010)";
814     break;
815     case -413:
816     answer = "D^-*(2010)";
817     break;
818     case 423:
819     answer = "D^0*(2010)";
820     break;
821     case -423:
822     answer = "D^0*(2010) bar";
823     break;
824     case 415:
825     answer = "D_2^+*(2460)";
826     break;
827     case -415:
828     answer = "D_2^-*(2460)";
829     break;
830     case 425:
831     answer = "D_2^0*(2460)";
832     break;
833     case -425:
834     answer = "D_2^0*(2460) bar";
835     break;
836     case 431:
837     answer = "D_s^+";
838     break;
839     case -431:
840     answer = "D_s^-";
841     break;
842     case 433:
843     answer = "D_s^+*";
844     break;
845     case -433:
846     answer = "D_s^-*";
847     break;
848     case 435:
849     answer = "D_s2^+*";
850     break;
851     case -435:
852     answer = "D_s2^-*";
853     break;
854    
855     case 511:
856     answer = "B0";
857     break;
858     case -511:
859     answer = "B0 bar";
860     break;
861     case -521:
862     answer = "B-";
863     break;
864     case 521:
865     answer = "B+";
866     break;
867     case 513:
868     answer = "B0*";
869     break;
870     case -513:
871     answer = "B0* bar";
872     break;
873     case -523:
874     answer = "B-*";
875     break;
876     case 523:
877     answer = "B+*";
878     break;
879     case 515:
880     answer = "B_2^0*";
881     break;
882     case -515:
883     answer = "B_2^0* bar";
884     break;
885     case -525:
886     answer = "B_2^-*";
887     break;
888     case 525:
889     answer = "B_2^+*";
890     break;
891     case 531:
892     answer = "B_s^0";
893     break;
894     case -531:
895     answer = "B_s^0 bar";
896     break;
897     case 533:
898     answer = "B_s^0*";
899     break;
900     case -533:
901     answer = "B_s^0* bar";
902     break;
903     case 535:
904     answer = "B_s2^0*";
905     break;
906     case -535:
907     answer = "B_s2^0* bar";
908     break;
909     case 541:
910     answer = "B_c^+";
911     break;
912     case -541:
913     answer = "B_c^-";
914     break;
915     case 543:
916     answer = "B_c^+*";
917     break;
918     case -543:
919     answer = "B_c^-*";
920     break;
921     case 545:
922     answer = "B_c2^+*";
923     break;
924     case -545:
925     answer = "B_c2^-*";
926     break;
927    
928     case 441:
929     answer = "eta_c(1S)";
930     break;
931     case 443:
932     answer = "JPsi(1S)";
933     break;
934     case 445:
935     answer = "chi_c2(1P)";
936     break;
937     case 551:
938     answer = "eta_b(1S)";
939     break;
940     case 553:
941     answer = "Upsilon(1S)";
942     break;
943     case 555:
944     answer = "chi_b2(1P)";
945     break;
946     case 557:
947     answer = "Upsilon(1D)";
948     break;
949    
950     case 2212:
951     answer = "proton";
952     break;
953     case 2112:
954     answer = "neutron";
955     break;
956     case 2224:
957     answer = "Delta++";
958     break;
959     case 2214:
960     answer = "Delta+";
961     break;
962     case 2114:
963     answer = "Delta0";
964     break;
965     case 1114:
966     answer = "Delta-";
967     break;
968     case 3122:
969     answer = "Lambda";
970     break;
971     case 3222:
972     answer = "Sigma+";
973     break;
974     case 3212:
975     answer = "Sigma0";
976     break;
977     case 3112:
978     answer = "Sigma-";
979     break;
980     case 3224:
981     answer = "Sigma+*";
982     break;
983     case 3214:
984     answer = "Sigma0*";
985     break;
986     case 3114:
987     answer = "Sigma-*";
988     break;
989     case 3322:
990     answer = "Xi0";
991     break;
992     case 3312:
993     answer = "Xi-";
994     break;
995     case 3324:
996     answer = "Xi0*";
997     break;
998     case 3314:
999     answer = "Xi-*";
1000     break;
1001     case 3334:
1002     answer = "Omega-";
1003     break;
1004    
1005     default:
1006     answer = "";
1007     }
1008     return answer;
1009     }