AliPhysics  608b256 (608b256)
AliMCAnalysisUtils.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TMath.h>
18 #include <TList.h>
19 #include "TParticle.h"
20 #include "TDatabasePDG.h"
21 #include "TVector3.h"
22 
23 //---- ANALYSIS system ----
24 #include "AliAnalysisManager.h"
25 #include "AliInputEventHandler.h"
26 #include "AliMCEvent.h"
27 #include "AliGenPythiaEventHeader.h"
28 #include "AliVParticle.h"
29 #include "AliLog.h"
30 
31 #include "AliMCAnalysisUtils.h"
32 
34 ClassImp(AliMCAnalysisUtils) ;
36 
37 //________________________________________
39 //________________________________________
41 TObject(),
42 fCurrentEvent(-1),
43 fDebug(0),
44 fJetsList(new TList),
45 fMCGenerator(kPythia),
46 fMCGeneratorString("PYTHIA"),
47 fDaughMom(), fDaughMom2(),
48 fMotherMom(), fGMotherMom(),
49 fPyGenHead(NULL) ,
50 fPyGenName(""),
51 fPyProcessName(""),
52 fPyProcess(-1),
53 fPyFirstParticle(0),
54 fPyVersion(0),
55 fMinPartonicParent(5),
56 fMaxPartonicParent(8)
57 {}
58 
59 //_______________________________________
61 //_______________________________________
63 {
64  if (fJetsList)
65  {
66  fJetsList->Clear();
67  delete fJetsList ;
68  }
69 }
70 
71 //_____________________________________________________________________________________________
74 //_____________________________________________________________________________________________
76  const AliMCEvent* mcevent,
77  Int_t & ancPDG, Int_t & ancStatus,
78  TLorentzVector & momentum, TVector3 & prodVertex)
79 {
80  if ( index1 < 0 || index2 < 0 )
81  {
82  ancPDG = -10000;
83  ancStatus = -10000;
84  momentum.SetXYZT(0,0,0,0);
85  prodVertex.SetXYZ(-10,-10,-10);
86  //printf("\t Negative index (%d, %d)\n",index1,index2);
87  }
88 
89  Int_t label1[1000];
90  Int_t label2[1000];
91  label1[0]= index1;
92  label2[0]= index2;
93  Int_t counter1 = 0;
94  Int_t counter2 = 0;
95 
96  if(label1[0]==label2[0])
97  {
98  //printf("\t Already the same label: %d\n",label1[0]);
99  counter1=1;
100  counter2=1;
101  }
102  else
103  {
104  Int_t label=label1[0];
105 
106  while(label > -1 && counter1 < 999)
107  {
108  counter1++;
109  AliVParticle * mom = mcevent->GetTrack(label);
110  if(mom)
111  {
112  label = mom->GetMother() ;
113  label1[counter1]=label;
114  }
115  //printf("\t 1) counter %d, mom label %d, pdg %d\n", counter1,label, mom->PdgCode());
116  }
117 
118  //printf("Org label2=%d,\n",label2[0]);
119  label=label2[0];
120 
121  while(label > -1 && counter2 < 999)
122  {
123  counter2++;
124  AliVParticle * mom = mcevent->GetTrack(label);
125  if(mom)
126  {
127  label = mom->GetMother() ;
128  label2[counter2]=label;
129  }
130  //printf("\t 2) counter %d, mom label %d, pdg %d \n", counter2,label, mom->PdgCode());
131  }
132  }//First labels not the same
133 
134  if((counter1==999 || counter2==999))
135  AliWarning(Form("Genealogy too large, generations: cluster1: %d, cluster2= %d", counter1, counter2));
136 
137  //printf("CheckAncestor:\n");
138 
139  //Int_t commonparents = 0;
140  Int_t ancLabel = -1;
141  //printf("counters %d %d \n",counter1, counter2);
142  for (Int_t c1 = 0; c1 < counter1; c1++)
143  {
144  for (Int_t c2 = 0; c2 < counter2; c2++)
145  {
146  if ( label1[c1]==label2[c2] && label1[c1]>-1 &&
147  ancLabel < label1[c1] ) // make sure to take the first common parent, not needed since array is ordered, just in case
148  {
149  ancLabel = label1[c1];
150  //commonparents++;
151 
152  AliVParticle * mom = mcevent->GetTrack(label1[c1]);
153 
154  if (mom)
155  {
156  ancPDG = mom->PdgCode();
157  ancStatus = mom->MCStatusCode();
158  momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
159  prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
160  //printf("Ancestor label %d PDG %d, status %d\n",ancLabel,ancPDG,ancStatus);
161  }
162 
163  //First ancestor found, end the loops
164  //counter1=0;
165  //counter2=0;
166  }//Ancestor found
167  }//second cluster loop
168  }//first cluster loop
169 
170  //printf("common ancestors %d\n",commonparents);
171 
172  if(ancLabel < 0)
173  {
174  //printf("No ancestor found!\n");
175  ancPDG = -10000;
176  ancStatus = -10000;
177  momentum.SetXYZT(0,0,0,0);
178  prodVertex.SetXYZ(-10,-10,-10);
179  }
180 
181  return ancLabel;
182 }
183 
184 //____________________________________________________________________________________________________
192 //_____________________________________________________________________________________________________
193 Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, AliMCEvent* mcevent,
194  TString selectHeaderName, Float_t clusE)
195 {
196  Int_t labels [] = { label };
197  UShort_t edep[] = { 0 };
198 
199  return CheckOrigin(labels, edep, 1, mcevent, selectHeaderName, clusE);
200 }
201 
202 //__________________________________________________________________________________________
220 //__________________________________________________________________________________________
221 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t *labels, const UShort_t * edepFrac, Int_t nlabels,
222  AliMCEvent* mcevent, TString selectHeaderName,
223  Float_t clusE, const TObjArray* arrayCluster)
224 {
225  if( nlabels <= 0 )
226  {
227  AliWarning("No MC labels available, please check!!!");
228  return kMCBadLabel;
229  }
230 
231  if ( !mcevent )
232  {
233  AliDebug(1,"MCEvent is not available, check analysis settings in configuration file, do nothing with MC!!");
234  return -1;
235  }
236 
237  Int_t tag = 0;
238  Int_t nprimaries = mcevent->GetNumberOfTracks();
239 
240  // Bad label
241  if ( labels[0] < 0 || labels[0] >= nprimaries )
242  {
243  if(labels[0] < 0)
244  AliWarning(Form("*** bad label ***: label %d", labels[0]));
245 
246  if(labels[0] >= nprimaries)
247  AliWarning(Form("*** large label ***: label %d, n tracks %d", labels[0], nprimaries));
248 
249  SetTagBit(tag,kMCBadLabel);
250 
251  return tag;
252  } // Bad label
253 
255  //
256  if(fMCGenerator == kPythia)
257  {
258  CheckAndGetPythiaEventHeader(mcevent,selectHeaderName);
259 
260 // printf("CheckOrigin() - Name <%s>, include name? <%s>, Pythia version %d,"
261 // " process %d, process name <%s>, first pythia particle %d\n",
262 // fPyGenName.Data(), selectHeaderName.Data(), fPyVersion,
263 // fPyProcess, fPyProcessName.Data(),fPyFirstParticle);
264 
265  // Select the index range where the prompt photon partonic parent can be
266  // it depends on the pythia version
267  }
268  //
270 
271  // Most significant particle contributing to the cluster
272  Int_t label=labels[0];
273 
274  // Mother
275  AliVParticle * mom = mcevent->GetTrack(label);
276  Int_t iMom = label;
277  Int_t mPdgSign = mom->PdgCode();
278  Int_t mPdg = TMath::Abs(mPdgSign);
279  Int_t mStatus = mom->MCStatusCode() ;
280  Int_t iParent = mom->GetMother() ;
281 
282  //if(label < 8 && fMCGenerator != kBoxLike) AliDebug(1,Form("Mother is parton %d\n",iParent));
283 
284  //GrandParent
285  AliVParticle * firstGparent = NULL ;
286  AliVParticle * parent = NULL ;
287  Int_t pPdg =-1;
288  Int_t pStatus =-1;
289  if(iParent >= 0)
290  {
291  parent = mcevent->GetTrack(iParent);
292  pPdg = TMath::Abs(parent->PdgCode());
293  pStatus = parent->MCStatusCode();
294  }
295  else AliDebug(1,Form("Parent with label %d",iParent));
296 
297  AliDebug(2,"Cluster most contributing mother and its parent:");
298  AliDebug(2,Form("\t Mother label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d",
299  iMom , mPdg, mStatus, mom->IsPrimary() , mom->IsPhysicalPrimary()));
300  AliDebug(2,Form("\t Parent label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d",
301  iParent, pPdg, pStatus, parent?parent->IsPrimary():-1, parent?parent->IsPhysicalPrimary():-1));
302 
303  //Check if mother is converted, if not, get the first non converted mother
304  if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus==0)
305  {
307 
308  // Check if the mother is photon or electron with status not stable
309  while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
310  {
311  // Mother
312  iMom = mom->GetMother();
313 
314  if(iMom < 0)
315  {
316  AliInfo(Form("pdg = %d, mother = %d, skip",pPdg,iMom));
317  break;
318  }
319 
320  mom = mcevent->GetTrack(iMom);
321  mPdgSign = mom->PdgCode();
322  mPdg = TMath::Abs(mPdgSign);
323  mStatus = mom->MCStatusCode() ;
324  iParent = mom->GetMother() ;
325  //if(label < 8 ) AliDebug(1, Form("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent));
326 
327  // GrandParent
328  if(iParent >= 0 && parent)
329  {
330  parent = mcevent->GetTrack(iParent);
331  pPdg = TMath::Abs(parent->PdgCode());
332  pStatus = parent->MCStatusCode();
333  }
334  // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
335  // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
336 
337  }//while
338 
339  AliDebug(2,"Converted photon/electron:");
340  AliDebug(2,Form("\t Mother label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d"
341  ,iMom , mPdg, mStatus, mom->IsPrimary() , mom->IsPhysicalPrimary()));
342  AliDebug(2,Form("\t Parent label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d"
343  ,iParent, pPdg, pStatus, parent?parent->IsPrimary():-1, parent?parent->IsPhysicalPrimary():-1));
344 
345  } // mother and parent are electron or photon and have status 0 and parent is photon or electron
346  else if((mPdg == 22 || mPdg == 11) && mStatus==0)
347  {
348  // Still a conversion but only one electron/photon generated. Just from hadrons
349  if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
350  pPdg == 2212 || pPdg == 130 || pPdg == 13 )
351  {
353  iMom = mom->GetMother();
354 
355  if(iMom < 0)
356  {
357  AliInfo(Form("pdg = %d, mother = %d, skip",pPdg,iMom));
358  }
359  else
360  {
361  mom = mcevent->GetTrack(iMom);
362  mPdgSign = mom->PdgCode();
363  mPdg = TMath::Abs(mPdgSign);
364  mStatus = mom->MCStatusCode() ;
365 
366  AliDebug(2,"Converted hadron:");
367  AliDebug(2,Form("\t Mother label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d",
368  iMom, mPdg, mStatus, mom->IsPrimary(), mom->IsPhysicalPrimary()));
369  }
370  } // hadron converted
371 
372  //Comment for next lines, we do not check the parent of the hadron for the moment.
373  //iParent = mom->GetMother() ;
374  //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
375 
376  //GrandParent
377  //if(iParent >= 0){
378  // parent = mcevent->GetTrack(iParent);
379  // pPdg = TMath::Abs(parent->PdgCode());
380  //}
381  }
382 
383  //printf("Final mother mPDG %d\n",mPdg);
384 
385  // conversion into electrons/photons checked
386 
387  //first check for typical charged particles
388  if (mPdg == 13) SetTagBit(tag,kMCMuon);
389  else if(mPdg == 211) SetTagBit(tag,kMCPion);
390  else if(mPdg == 321) SetTagBit(tag,kMCKaon);
391  else if(mPdgSign == 2212) SetTagBit(tag,kMCProton);
392  else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron);
393  else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
394  else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
395 
396  //check for pi0 and eta (shouldn't happen unless their decays were turned off)
397  else if(mPdg == 111)
398  {
399  SetTagBit(tag,kMCPi0Decay);
400 
401  AliDebug(2,"First mother is directly pi0, not decayed by generator");
402 
403  CheckOverlapped2GammaDecay(labels, edepFrac, nlabels, iMom, clusE, mcevent, tag); //set to kMCPi0 if 2 gammas in same cluster
404  }
405  else if(mPdg == 221)
406  {
407  SetTagBit(tag,kMCEtaDecay);
408 
409  AliDebug(2,"First mother is directly eta, not decayed by generator");
410 
411  CheckOverlapped2GammaDecay(labels, edepFrac, nlabels, iMom, clusE, mcevent, tag); //set to kMCEta if 2 gammas in same cluster
412  }
413  //Photons
414  else if(mPdg == 22)
415  {
416  SetTagBit(tag,kMCPhoton);
417 
418  if ( pPdg == 111 )
419  {
420  SetTagBit(tag,kMCPi0Decay);
421 
422  AliDebug(2,"Generator pi0 decay photon");
423 
424  // Set to kMCPi0 if 2 gammas in same cluster
425  CheckOverlapped2GammaDecay(labels, edepFrac, nlabels, iParent, clusE, mcevent, tag);
426 
427  // In case it did not merge, check if the decay companion is lost
429  {
430  CheckLostDecayPair(arrayCluster,iMom, iParent, mcevent, tag);
431  }
432 
433  //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
434  }
435  else if ( pPdg == 221 )
436  {
437  SetTagBit(tag, kMCEtaDecay);
438 
439  AliDebug(2,"Generator eta decay photon");
440 
441  // Set to kMCEta if 2 gammas in same cluster
442  CheckOverlapped2GammaDecay(labels, edepFrac, nlabels, iParent, clusE, mcevent, tag);
443 
444  // In case it did not merge, check if the decay companion is lost
446  CheckLostDecayPair(arrayCluster,iMom, iParent, mcevent, tag);
447  }
448  else if ( pPdg > 100 )
449  {
451 
452  AliDebug(2,Form("Generator decay photon from parent pdg %d",pPdg));
453  }
454  else if( mom->IsPhysicalPrimary() && mStatus == 1 &&
455  ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle
456  {
457  Int_t firstpPdg =-1 ;
458  Int_t firstPhotonLabel =-1 ;
459  Int_t gparentlabel =-1 ;
460  Bool_t ok = kFALSE;
461 
462 // printf("=== Photon Label %d, Mom label %d, iParent %d, pPdg %d, pT %2.3f; "
463 // "process %d (%s), min parent %d, max parent %d, pythia v%d\n",
464 // label,iMom, iParent, pPdg, mom->Pt(), fPyProcess,
465 // fPyProcessName.Data(), fMinPartonicParent,fMaxPartonicParent, fPyVersion);
466 
467  //PrintAncestry(mcevent,iMom,10);
468 
469  if ( fPyVersion == 6 || fMCGenerator == kHerwig ) // to be checked for herwig
470  {
471  if( iParent-fPyFirstParticle < fMaxPartonicParent &&
473  {
474  //outgoing partons
475  if ( pPdg == 22 ) { SetTagBit(tag,kMCPrompt) ; /*printf("\t \t Prompt6\n")*/ ; }
476  else { SetTagBit(tag,kMCFragmentation); /*printf("\t \t Fragment66\n");*/ }
477  }//Outgoing partons
478  else if( iParent <= fMinPartonicParent ) // Pythia6
479  {
480  SetTagBit(tag, kMCISR); //Initial state radiation
481  //printf("\t \t ISR\n");
482  }
483  else if ( pPdg < 23 )
484  {
486 
487  AliDebug(2,Form("Generator fragmentation photon from parent pdg %d",firstpPdg));
488  //printf("\t \t Fragment6\n");
489  }
490  else
491  {
492  AliDebug(2,Form("Generator physical primary (pythia/herwig) unknown photon from parent pdg %d",firstpPdg));
493  //printf("\t \t Unknown\n");
494  SetTagBit(tag,kMCUnknown);
495  }
496  }
497  else if ( fPyVersion == 8 )
498  {
499  GetFirstMotherWithPDG(iMom,22,mcevent, ok, firstPhotonLabel,gparentlabel);
500 // printf("\t First Gen %d, First photon: %d, parent %d; first photon - first gen %d, parent-firstGen %d\n",
501 // fPyFirstParticle, firstPhotonLabel,gparentlabel,
502 // firstPhotonLabel-fPyFirstParticle, gparentlabel-fPyFirstParticle);
503  firstGparent = mcevent->GetTrack(firstPhotonLabel);
504  firstpPdg = TMath::Abs(firstGparent->PdgCode());
505 
506  if( firstPhotonLabel-fPyFirstParticle < fMaxPartonicParent &&
507  firstPhotonLabel-fPyFirstParticle > fMinPartonicParent )
508  {
509  SetTagBit(tag,kMCPrompt);
510  //printf("\t \t Prompt8\n");
511  }
512  else
513  {
515 
516  AliDebug(2,Form("Generator fragmentation photon from parent pdg %d",firstpPdg));
517  //printf("\t \t Fragment8\n");
518  }
519  }
520  } // Physical primary
521  else
522  {
523  AliDebug(2,Form("Generator unknown unphysical photon from parent pdg %d",pPdg));
524 
525  SetTagBit(tag,kMCUnknown);
526  }
527 // //Old Herwig selection for ESDs, maybe should be considered.
528 // if(fMCGenerator == kHerwig)
529 // {
530 // if(pStatus < 197)
531 // {//Not decay
532 // while(1)
533 // {
534 // if(parent)
535 // {
536 // if(parent->GetMother()<=5) break;
537 // iParent = parent->GetMother();
538 // parent = mcevent->GetTrack(iParent);
539 // pStatus = parent->MCStatusCode();
540 // pPdg = TMath::Abs(parent->PdgCode());
541 // } else break;
542 // }//Look for the parton
543 //
544 // if(iParent < 8 && iParent > 5)
545 // {
546 // if(pPdg == 22) SetTagBit(tag,kMCPrompt);
547 // else SetTagBit(tag,kMCFragmentation);
548 // }
549 // else SetTagBit(tag,kMCISR);//Initial state radiation
550 // }//Not decay
551 // else SetTagBit(tag,kMCUnknown);
552 // }//HERWIG
553 
554  } // Mother Photon
555 
556  // Electron check. Where did that electron come from?
557  else if(mPdg == 11)
558  {
559  //electron
560  if(pPdg == 11 && parent)
561  {
562  Int_t iGrandma = parent->GetMother();
563  if(iGrandma >= 0)
564  {
565  AliVParticle * gma = mcevent->GetTrack(iGrandma);
566  Int_t gPdg = TMath::Abs(gma->PdgCode());
567 
568  if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
569  else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
570  }
571  }
572 
573  SetTagBit(tag,kMCElectron);
574 
575  AliDebug(1,"Checking ancestors of electrons");
576 
577  if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay
578  else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay
579  else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
580  else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
581  {
582  //c-hadron decay check
583  if(parent)
584  {
585  Int_t iGrandma = parent->GetMother();
586  if(iGrandma >= 0)
587  {
588  AliVParticle * gma = mcevent->GetTrack(iGrandma); //charm's mother
589  Int_t gPdg = TMath::Abs(gma->PdgCode());
590  if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
591  else SetTagBit(tag,kMCEFromC); //c-hadron decay
592  }
593  else SetTagBit(tag,kMCEFromC); //c-hadron decay
594  }//parent
595  } else
596  { //prompt or other decay
597 
598  TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
599  TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
600 
601  AliDebug(1,Form("Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)",foo->GetName(), pPdg,foo1->GetName(),mPdg));
602 
603  if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
604  else SetTagBit(tag,kMCOtherDecay);
605  }
606  }// electron check
607  // cluster was made by something else
608  else
609  {
610  AliDebug(1,Form("\t Setting kMCUnknown for cluster with pdg = %d, Parent pdg = %d",mPdg,pPdg));
611  SetTagBit(tag,kMCUnknown);
612  }
613 
614  return tag;
615 }
616 
617 //_________________________________________________________________________________________
628 //_________________________________________________________________________________________
629 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const UShort_t *edepFrac,
630  Int_t nlabels, Int_t mesonIndex, Float_t clusE,
631  const AliMCEvent* mcevent, Int_t & tag)
632 {
633  if(labels[0] < 0 || labels[0] > mcevent->GetNumberOfTracks() || nlabels <= 1)
634  {
635  AliDebug(2,Form("Exit : label[0] %d, n primaries %d, nlabels %d",labels[0],mcevent->GetNumberOfTracks(), nlabels));
636  return;
637  }
638 
639 
640  AliVParticle * meson = mcevent->GetTrack(mesonIndex);
641  Int_t mesonPdg = meson->PdgCode();
642  if ( mesonPdg != 111 && mesonPdg != 221 )
643  {
644  AliWarning(Form("Wrong pi0/eta PDG : %d",mesonPdg));
645  return;
646  }
647 
648  AliDebug(2,Form("pdg %d, label %d, ndaughters %d", mesonPdg, mesonIndex, meson->GetNDaughters()));
649 
650 // Bool_t convOrg = CheckTagBit(tag,kMCConversion);
651 // printf("pdg %d, label %d, ndaughters %d, E prim %2.2f, E clus %2.2f, conv? %d\n",
652 // mesonPdg, mesonIndex, meson->GetNDaughters(),meson->E(),clusE, convOrg);
653 
654  // Reject very low energy mesons unlikely to merge
655  // Pi0s merge from ~5 GeV and Eta's from ~20 GeV
656  if ( ( meson->E() < 5 && mesonPdg == 111 ) ||
657  ( meson->E() < 20 && mesonPdg == 221 ))
658  {
659  AliDebug(3,Form("Too low generated meson Energy to decay overlap: pdg %d, E prim %2.2f\n",
660  mesonPdg,meson->E()));
661 
662  return;
663  }
664 
665  // Check energy of input cluster and meson, if too low cluster energy,
666  // discard overlap, it will be likely some weird large conversion
667  if ( clusE/meson->E() < 0.75 )
668  {
669  AliDebug(3,Form("Too low generated vs reconstructed meson Energy : pdg %d, E prim %2.2f, E clus %2.2f\n",
670  mesonPdg,meson->E(),clusE));
671 
672  return;
673  }
674 
675  // Get the daughters
676  if ( meson->GetNDaughters() != 2 )
677  {
678 // if(meson->GetNDaughters()>2)
679 // {
680 // printf("xx More than 2 daughters <%d> for pdg %d with E %2.2f, pT %2.2f:\n",
681 // meson->GetNDaughters(),meson->PdgCode(),meson->E(),meson->Pt());
682 // for(Int_t idaugh = 0; idaugh < meson->GetNDaughters(); idaugh++)
683 // {
684 // if(meson->GetDaughterLabel(idaugh) < 0)
685 // {
686 // printf("\t no daughther with index %d:\n",idaugh);
687 // continue;
688 // }
689 // AliVParticle * daugh = mcevent->GetTrack(meson->GetDaughterLabel(idaugh));
690 // printf("\t daughther pdg %d, E %2.2f, pT %2.2f:\n",daugh->PdgCode(),daugh->E(),daugh->Pt());
691 // }
692 // }
693 
694  AliDebug(2,Form("Not overlapped. Number of daughters is %d, not 2",meson->GetNDaughters()));
695 
696  return;
697  }
698 
699  Int_t iPhoton0 = meson->GetDaughterLabel(0);
700  Int_t iPhoton1 = meson->GetDaughterLabel(1);
701 
702  //if((iPhoton0 == -1) || (iPhoton1 == -1))
703  //{
704  // if(fDebug > 2)
705  // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
706  // return;
707  //}
708 
709  AliVParticle *photon0 = mcevent->GetTrack(iPhoton0);
710  AliVParticle *photon1 = mcevent->GetTrack(iPhoton1);
711 
712  // Check if both daughters are photons
713  if ( photon0->PdgCode() != 22 && photon1->PdgCode() != 22 )
714  {
715  AliWarning(Form("Not overlapped. PDG: daughter 1 = %d, of daughter 2 = %d",photon0->PdgCode(),photon1->PdgCode()));
716 
717  return;
718  }
719 
720  AliDebug(2,Form("Daughter labels : photon0 = %d, photon1 = %d",iPhoton0,iPhoton1));
721 
722 
723  // Avoid cases with large opening angle, usually weird conversions, with one electron/positron
724  // falling into the other photon cluster
725  fDaughMom .SetPxPyPzE(photon0->Px(),photon0->Py(),photon0->Pz(),photon0->E());
726  fDaughMom2.SetPxPyPzE(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());
727  Double_t angle = fDaughMom.Angle(fDaughMom2.Vect());
728 
729  // Reject pairs with opening angle smaller than 4 EMCal cells size
730  // careful in case of PHOS, but it should be safe enough
731  if ( angle/0.0143 > 4 )
732  {
733  AliDebug(2,Form("Not an overlapped pair, opening angle %f rad, %2.2f EMCal cells",angle,angle/0.0143));
734 
735  return;
736  }
737 
738  // Check if both photons contribute to the cluster
739  Bool_t okPhoton0 = kFALSE;
740  Bool_t okPhoton1 = kFALSE;
741 
742  AliDebug(3,"Labels loop:");
743 
744  Bool_t conversion = kFALSE;
745 
746 // TLorentzVector d0, d1;
747 // photon0->Momentum(d0);
748 // photon1->Momentum(d1);
749 // Float_t openAngle = d0.Angle(d1.Vect());
750 
751  for(Int_t i = 0; i < nlabels; i++)
752  {
753  AliDebug(3, Form("\t label %d/%d: %d, ok? %d, %d", i, nlabels, labels[i], okPhoton0, okPhoton1));
754 
755  if ( labels[i] < 0 ) continue;
756 
757  // If we already found both, break the loop
758  if ( okPhoton0 && okPhoton1 ) break;
759 
760  Int_t index = labels[i];
761  if (iPhoton0 == index)
762  {
763  okPhoton0 = kTRUE;
764  continue;
765  }
766  else if (iPhoton1 == index)
767  {
768  okPhoton1 = kTRUE;
769  continue;
770  }
771 
772  // Trace back the mother in case it was a conversion
773  if ( index >= mcevent->GetNumberOfTracks() )
774  {
775  AliWarning(Form("Particle index %d larger than size of list %d!!",index,mcevent->GetNumberOfTracks()));
776  continue;
777  }
778 
779  AliVParticle * daught = mcevent->GetTrack(index);
780 
781 // printf("\t parent index %d, pdg %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n",
782 // index,daught->PdgCode(),daught->E(),daught->Pt(),daught->Eta()/0.0143,daught->Phi()/0.0143);
783 
784  Int_t tmpindex = daught->GetMother();
785  AliDebug(3,Form("Conversion? : mother %d",tmpindex));
786 
787  while ( tmpindex >= 0 )
788  {
789  // MC particle of interest is the mother
790  AliDebug(3,Form("\t parent index %d",tmpindex));
791  daught = mcevent->GetTrack(tmpindex);
792  //printf("tmpindex %d\n",tmpindex);
793 
794 // if ( daught->Pt()>0.1 )
795 // printf("\t new parent index %d, pdg %d, E %2.2f, pT %2.2f, eta %2.2f, phi %2.2f\n",
796 // tmpindex,daught->PdgCode(),daught->E(),daught->Pt(),daught->Eta()/0.0143,daught->Phi()/0.0143);
797 
798  if (iPhoton0 == tmpindex)
799  {
800  conversion = kTRUE;
801  okPhoton0 = kTRUE;
802  break;
803  }
804  else if (iPhoton1 == tmpindex)
805  {
806  conversion = kTRUE;
807  okPhoton1 = kTRUE;
808  break;
809  }
810 
811  tmpindex = daught->GetMother();
812 
813  }//While to check if pi0/eta daughter was one of these contributors to the cluster
814 
815  //if(i == 0 && (!okPhoton0 && !okPhoton1)) AliDebug(1,"Something happens, first label should be from a photon decay!");
816 
817  }//loop on list of labels
818 
819  //If both photons contribute tag as the corresponding meson.
820  if ( okPhoton0 && okPhoton1 )
821  {
822  AliDebug(2,Form("%s OVERLAPPED DECAY",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()));
823 
824 // if ( meson->E() < 7 )
825 // {
826 // printf("pdg %d, label %d, ndaughters %d, E %2.2f, pT %2.2f\n",
827 // mesonPdg, mesonIndex, meson->GetNDaughters(),meson->E(),meson->Pt());
828 // printf("\t %s OVERLAPPED DECAY, conversion %d?\n",
829 // (TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName(),conversion);
830 //
831 // printf("\t Daughter labels : photon0 = %d, E %2.2f, pT %2.2f, photon1 = %d, E %2.2f, pT %2.2f; angle/cell %2.2f\n",
832 // iPhoton0,photon0->E(),photon0->Pt(),
833 // iPhoton1,photon1->E(),photon1->Pt(),
834 // angle/0.0143);
835 // }
836 
837  if(!CheckTagBit(tag,kMCConversion) && conversion)
838  {
839  AliDebug(2,"Second decay photon produced a conversion");
840  SetTagBit(tag,kMCConversion) ;
841  }
842 
843  if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
844  else SetTagBit(tag,kMCEta);
845 
846 // printf("\t Meson E %2.2f, pT %2.2f, Overlapped open angle %2.2f deg, cells %2.1f, conversion %d\n",
847 // meson->E(), meson->Pt(),openAngle*TMath::RadToDeg(),openAngle/0.015,conversion);
848 //
849 // Float_t phi = photon0->Phi()*TMath::RadToDeg();
850 // if(phi < 0) phi+=360;
851 // printf("\t \t daughther1: label %d, E %2.2f, pT %2.2f, eta%2.2f, phi %2.2f, status %d, phys prim %d:\n",
852 // iPhoton0, photon0->E(),photon0->Pt(),
853 // photon0->Eta(),phi,
854 // photon0->MCStatusCode(),photon0->IsPhysicalPrimary());
855 //
856 // phi = photon1->Phi()*TMath::RadToDeg();
857 // if(phi < 0) phi+=360;
858 // printf("\t \t daughther2: label %d, E %2.2f, pT %2.2f, eta%2.2f, phi %2.2f, status %d, phys prim %d:\n",
859 // iPhoton1, photon1->E(),photon1->Pt(),
860 // photon1->Eta(),phi,
861 // photon1->MCStatusCode(),photon1->IsPhysicalPrimary());
862  }
863 }
864 
865 //________________________________________________________________________________________________________
867 //________________________________________________________________________________________________________
868 void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster, Int_t iMom, Int_t iParent,
869  const AliMCEvent* mcevent, Int_t & tag)
870 {
871  if(!arrayCluster || iMom < 0 || iParent < 0|| !mcevent) return;
872 
873  AliVParticle * parent = mcevent->GetTrack(iParent);
874 
875  //printf("*** Check label %d with parent %d\n",iMom, iParent);
876 
877  if(parent->GetNDaughters()!=2)
878  {
880  //printf("\t ndaugh = %d\n",parent->GetNDaughters());
881  return ;
882  }
883 
884  Int_t pairLabel = -1;
885  if ( iMom != parent->GetDaughterLabel(0) ) pairLabel = parent->GetDaughterLabel(0);
886  else if( iMom != parent->GetDaughterLabel(1) ) pairLabel = parent->GetDaughterLabel(1);
887 
888  if(pairLabel<0)
889  {
890  //printf("\t pair Label not found = %d\n",pairLabel);
892  return ;
893  }
894 
895  //printf("\t *** find pair %d\n",pairLabel);
896 
897  for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
898  {
899  AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
900  //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
901 
902  for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
903  {
904  Int_t label = cluster->GetLabels()[ilab];
905 
906  //printf("\t \t label %d\n",label);
907 
908  if ( label==pairLabel )
909  {
910  //printf("\t \t Pair found\n");
912  return ;
913  }
914  else if ( label == iParent || label == iMom || label < 0 )
915  {
916  AliDebug(1,Form("Skip checking label %d, (iParent %d, iMom %d)",label,iParent,iMom));
917  continue;
918  }
919  else // check the ancestry
920  {
921  AliVParticle * mother = mcevent->GetTrack(label);
922 
923  if ( !mother )
924  {
925  AliInfo(Form("MC Mother not available for label %d",label));
926  continue;
927  }
928 
929  Int_t momPDG = TMath::Abs(mother->PdgCode());
930  if ( momPDG!=11 && momPDG!=22 ) continue;
931 
932  // Check if "mother" of entity is converted, if not, get the first non converted mother
933  Int_t iParentClus = mother->GetMother();
934  if(iParentClus < 0) continue;
935 
936  AliVParticle * parentClus = mcevent->GetTrack(iParentClus);
937  if(!parentClus) continue;
938 
939  Int_t parentClusPDG = TMath::Abs(parentClus->PdgCode());
940  Int_t parentClusStatus = parentClus->MCStatusCode();
941 
942  if ( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0 )
943  {
944  //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
945  continue;
946  }
947 
948  //printf("\t \t Conversion\n");
949 
950  // Check if the mother is photon or electron with status not stable
951  while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
952  {
953  //New Mother
954  label = iParentClus;
955  momPDG = parentClusPDG;
956 
957  iParentClus = parentClus->GetMother();
958  if ( iParentClus < 0 ) break;
959 
960  parentClus = mcevent->GetTrack(iParentClus);
961  if ( !parentClus ) break;
962 
963  parentClusPDG = TMath::Abs(parentClus->PdgCode());
964  parentClusStatus = parentClus->MCStatusCode() ;
965  }//while
966 
967  if ( (momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel) )
968  {
970  //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
971  return ;
972  }
973  else
974  {
975  //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
976  continue;
977  }
978 
979  } // last else, check the ancestry
980  } // label loop
981  } // cluster loop
982 
984 }
985 
986 //_____________________________________________________________________
990 //_____________________________________________________________________
991 TList * AliMCAnalysisUtils::GetJets(AliMCEvent* mcevent, Bool_t check)
992 {
993  if ( !fJetsList || !check ) return fJetsList ;
994 
995  if ( fJetsList ) fJetsList->Clear();
996  else fJetsList = new TList;
997 
998  Int_t nTriggerJets = 0;
999  Float_t tmpjet[]={0,0,0,0};
1000 
1001  // Get outgoing partons
1002  if(mcevent->GetNumberOfTracks() < 8) return fJetsList;
1003 
1004  AliVParticle * parton1 = mcevent->GetTrack(6);
1005  AliVParticle * parton2 = mcevent->GetTrack(7);
1006 
1007  AliDebug(2,Form("Parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1008  parton1->GetName(),parton1->Pt(),parton1->E(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()));
1009  AliDebug(2,Form("Parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1010  parton2->GetName(),parton2->Pt(),parton2->E(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()));
1011 
1012  // //Trace the jet from the mother parton
1013  // Float_t pt = 0;
1014  // Float_t pt1 = 0;
1015  // Float_t pt2 = 0;
1016  // Float_t e = 0;
1017  // Float_t e1 = 0;
1018  // Float_t e2 = 0;
1019  // TParticle * tmptmp = new TParticle;
1020  // for(Int_t i = 0; i< stack->GetNprimary(); i++){
1021  // tmptmp = stack->Particle(i);
1022 
1023  // if(tmptmp->GetStatusCode() == 1){
1024  // pt = tmptmp->Pt();
1025  // e = tmptmp->Energy();
1026  // Int_t imom = tmptmp->GetMother();
1027  // Int_t imom1 = 0;
1028  // //printf("1st imom %d\n",imom);
1029  // while(imom > 5){
1030  // imom1=imom;
1031  // tmptmp = stack->Particle(imom);
1032  // imom = tmptmp->GetMother();
1033  // //printf("imom %d \n",imom);
1034  // }
1035  // //printf("Last imom %d %d\n",imom1, imom);
1036  // if(imom1 == 6) {
1037  // pt1+=pt;
1038  // e1+=e;
1039  // }
1040  // else if (imom1 == 7){
1041  // pt2+=pt;
1042  // e2+=e; }
1043  // }// status 1
1044 
1045  // }// for
1046 
1047  // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
1048 
1049  //Get the jet, different way for different generator
1050  //PYTHIA
1051  if(fMCGenerator == kPythia)
1052  {
1053  TParticle * jet = 0x0;
1054  nTriggerJets = fPyGenHead->NTriggerJets();
1055  AliDebug(2,Form("PythiaEventHeader: Njets: %d",nTriggerJets));
1056 
1057  for(Int_t i = 0; i< nTriggerJets; i++)
1058  {
1059  fPyGenHead->TriggerJet(i, tmpjet);
1060 
1061  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1062 
1063  // Assign an outgoing parton as mother
1064  Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1065  Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1066 
1067  if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1068  else jet->SetFirstMother(6);
1069 
1070  //jet->Print();
1071  AliDebug(1,Form("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1072  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()));
1073  fJetsList->Add(jet);
1074  }
1075  }//Pythia triggered jets
1076  //HERWIG
1077  else if (fMCGenerator == kHerwig)
1078  {
1079  Int_t pdg = -1;
1080 
1081  // Check parton 1
1082  AliVParticle * tmp = parton1;
1083  if(parton1->PdgCode()!=22)
1084  {
1085  while(pdg != 94)
1086  {
1087  if(tmp->GetDaughterFirst()==-1) return fJetsList;
1088 
1089  tmp = mcevent->GetTrack(tmp->GetDaughterFirst());
1090  pdg = tmp->PdgCode();
1091  }//while
1092 
1093  // Add found jet to list
1094  TParticle *jet1 = new TParticle(94, 21, -1, -1, -1, -1, tmp->Px(),tmp->Py(),tmp->Pz(),tmp->E(), 0,0,0,0);//new TParticle(*tmp);
1095 
1096  jet1->SetFirstMother(6);
1097 
1098  fJetsList->Add(jet1);
1099 
1100  //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetDaughterFirst(), tmp->GetDaughterLast());
1101  //tmp = stack->Particle(tmp->GetDaughterFirst());
1102  //tmp->Print();
1103  //jet1->Print();
1104  AliDebug(1,Form("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1105  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()));
1106  }//not photon
1107 
1108  // Check parton 2
1109  pdg = -1;
1110  tmp = parton2;
1111  if(parton2->PdgCode()!=22)
1112  {
1113  while(pdg != 94)
1114  {
1115  if(tmp->GetDaughterFirst()==-1) return fJetsList;
1116 
1117  tmp = mcevent->GetTrack(tmp->GetDaughterFirst());
1118  pdg = tmp->PdgCode();
1119  }//while
1120 
1121  // Add found jet to list
1122  TParticle *jet2 = new TParticle(94, 21, -1, -1, -1, -1, tmp->Px(),tmp->Py(),tmp->Pz(),tmp->E(), 0,0,0,0);//new TParticle(*tmp);
1123 
1124  jet2->SetFirstMother(7);
1125 
1126  fJetsList->Add(jet2);
1127 
1128  //jet2->Print();
1129  AliDebug(2,Form("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1130  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()));
1131  //Int_t first = tmp->GetDaughterFirst();
1132  //Int_t last = tmp->GetDaughterLast();
1133  //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1134  // for(Int_t d = first ; d < last+1; d++){
1135  // tmp = stack->Particle(d);
1136  // if(i == tmp->GetMother())
1137  // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1138  // d,tmp->GetMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1139  // }
1140  //tmp->Print();
1141  }//not photon
1142  }//Herwig generated jets
1143 
1144  return fJetsList;
1145 }
1146 
1147 //______________________________________________________________________________
1149 //______________________________________________________________________________
1150 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
1151  const AliMCEvent* mcevent,
1152  Int_t & pdg, Int_t & status,
1153  Bool_t & ok, Int_t & daughlabel, TVector3 & prodVertex)
1154 {
1155  fDaughMom.SetPxPyPzE(0,0,0,0);
1156 
1157  if(!mcevent)
1158  {
1159  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
1160 
1161  ok = kFALSE;
1162  return fDaughMom;
1163  }
1164 
1165  Int_t nprimaries = mcevent->GetNumberOfTracks();
1166  if(label < 0 || label >= nprimaries)
1167  {
1168  ok = kFALSE;
1169  return fDaughMom;
1170  }
1171 
1172  AliVParticle * momP = mcevent->GetTrack(label);
1173  daughlabel = momP->GetDaughterLabel(idaugh);
1174 
1175  if(daughlabel < 0 || daughlabel >= nprimaries)
1176  {
1177  ok = kFALSE;
1178  return fDaughMom;
1179  }
1180 
1181  AliVParticle * daughP = mcevent->GetTrack(daughlabel);
1182  fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1183  pdg = daughP->PdgCode();
1184  status = daughP->MCStatusCode();
1185  prodVertex.SetXYZ(daughP->Xv(),daughP->Yv(),daughP->Zv());
1186 
1187  ok = kTRUE;
1188 
1189  return fDaughMom;
1190 }
1191 
1192 //______________________________________________________________________________________________________
1194 //______________________________________________________________________________________________________
1195 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent, Bool_t & ok)
1196 {
1197  Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1198 
1199  return GetMother(label,mcevent,pdg,status, ok,momlabel);
1200 }
1201 
1202 //_________________________________________________________________________________________
1204 //_________________________________________________________________________________________
1205 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent,
1206  Int_t & pdg, Int_t & status, Bool_t & ok)
1207 {
1208  Int_t momlabel = -1;
1209 
1210  return GetMother(label,mcevent,pdg,status, ok,momlabel);
1211 }
1212 
1213 //______________________________________________________________________________________________________
1215 //______________________________________________________________________________________________________
1216 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent,
1217  Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1218 {
1219  fMotherMom.SetPxPyPzE(0,0,0,0);
1220 
1221  if(!mcevent)
1222  {
1223  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1224 
1225  ok = kFALSE;
1226  return fMotherMom;
1227  }
1228 
1229  Int_t nprimaries = mcevent->GetNumberOfTracks();
1230  if(label < 0 || label >= nprimaries)
1231  {
1232  ok = kFALSE;
1233  return fMotherMom;
1234  }
1235 
1236  AliVParticle * momP = mcevent->GetTrack(label);
1237  fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1238  pdg = momP->PdgCode();
1239  status = momP->MCStatusCode();
1240  momlabel = momP->GetMother();
1241 
1242  ok = kTRUE;
1243 
1244  return fMotherMom;
1245 }
1246 
1247 //___________________________________________________________________________________
1249 //___________________________________________________________________________________
1251  const AliMCEvent* mcevent,
1252  Bool_t & ok, Int_t & momlabel)
1253 {
1254  fGMotherMom.SetPxPyPzE(0,0,0,0);
1255 
1256  if ( !mcevent )
1257  {
1258  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
1259 
1260  ok = kFALSE;
1261  return fGMotherMom;
1262  }
1263 
1264  Int_t nprimaries = mcevent->GetNumberOfTracks();
1265  if ( label < 0 || label >= nprimaries )
1266  {
1267  ok = kFALSE;
1268  return fGMotherMom;
1269  }
1270 
1271 
1272  AliVParticle * momP = mcevent->GetTrack(label);
1273 
1274  if(momP->PdgCode()==pdg)
1275  {
1276  AliDebug(2,"PDG of mother is already the one requested!");
1277  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1278 
1279  ok=kTRUE;
1280  return fGMotherMom;
1281  }
1282 
1283  Int_t grandmomLabel = momP->GetMother();
1284  Int_t grandmomPDG = -1;
1285  AliVParticle * grandmomP = 0x0;
1286 
1287  while (grandmomLabel >=0 )
1288  {
1289  grandmomP = mcevent->GetTrack(grandmomLabel);
1290  grandmomPDG = grandmomP->PdgCode();
1291  if(grandmomPDG==pdg)
1292  {
1293  //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1294  momlabel = grandmomLabel;
1295  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1296  break;
1297  }
1298 
1299  grandmomLabel = grandmomP->GetMother();
1300  }
1301 
1302  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found!",pdg));
1303 
1304  ok = kTRUE;
1305 
1306  return fGMotherMom;
1307 }
1308 
1309 //___________________________________________________________________________________
1312 //___________________________________________________________________________________
1314  const AliMCEvent* mcevent,
1315  Bool_t & ok, Int_t & momlabel,
1316  Int_t & gparentlabel)
1317 {
1318  fGMotherMom.SetPxPyPzE(0,0,0,0);
1319  momlabel = label;
1320 
1321  if ( !mcevent )
1322  {
1323  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
1324 
1325  ok = kFALSE;
1326  return fGMotherMom;
1327  }
1328 
1329  Int_t nprimaries = mcevent->GetNumberOfTracks();
1330  if ( label < 0 || label >= nprimaries )
1331  {
1332  ok = kFALSE;
1333  return fGMotherMom;
1334  }
1335 
1336  AliVParticle * momP = mcevent->GetTrack(label);
1337 
1338  Int_t grandmomLabel = momP->GetMother();
1339  gparentlabel = grandmomLabel;
1340  Int_t grandmomPDG = -1;
1341  AliVParticle * grandmomP = 0x0;
1342 
1343  while (grandmomLabel >=0 )
1344  {
1345  grandmomP = mcevent->GetTrack(grandmomLabel);
1346  grandmomPDG = grandmomP->PdgCode();
1347  if(grandmomPDG==pdg)
1348  {
1349  //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1350  momlabel = grandmomLabel;
1351  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1352  gparentlabel = grandmomP->GetMother();
1353  }
1354 
1355  grandmomLabel = grandmomP->GetMother();
1356  }
1357 
1358  ok = kTRUE;
1359 
1360  return fGMotherMom;
1361 }
1362 
1363 
1364 //______________________________________________________________________________________________
1366 //______________________________________________________________________________________________
1367 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliMCEvent* mcevent,
1368  Int_t & pdg, Int_t & status, Bool_t & ok,
1369  Int_t & grandMomLabel, Int_t & greatMomLabel)
1370 {
1371  fGMotherMom.SetPxPyPzE(0,0,0,0);
1372 
1373  if ( !mcevent )
1374  {
1375  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1376 
1377  ok = kFALSE;
1378  return fGMotherMom;
1379  }
1380 
1381  Int_t nprimaries = mcevent->GetNumberOfTracks();
1382  if ( label < 0 || label >= nprimaries )
1383  {
1384  ok = kFALSE;
1385  return fGMotherMom;
1386  }
1387 
1388  AliVParticle * momP = mcevent->GetTrack(label);
1389 
1390  grandMomLabel = momP->GetMother();
1391 
1392  AliVParticle * grandmomP = 0x0;
1393 
1394  if(grandMomLabel >=0 )
1395  {
1396  grandmomP = mcevent->GetTrack(grandMomLabel);
1397  pdg = grandmomP->PdgCode();
1398  status = grandmomP->MCStatusCode();
1399 
1400  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1401  greatMomLabel = grandmomP->GetMother();
1402 
1403  }
1404 
1405  ok = kTRUE;
1406 
1407  return fGMotherMom;
1408 }
1409 
1410 //_______________________________________________________________________________________________________________
1412 //_______________________________________________________________________________________________________________
1414  Float_t & asym, Float_t & angle, Bool_t & ok)
1415 {
1416  if(!mcevent)
1417  {
1418  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1419 
1420  ok = kFALSE;
1421  return;
1422  }
1423 
1424  Int_t nprimaries = mcevent->GetNumberOfTracks();
1425  if ( label < 0 || label >= nprimaries )
1426  {
1427  ok = kFALSE;
1428  return ;
1429  }
1430 
1431  AliVParticle * momP = mcevent->GetTrack(label);
1432 
1433  Int_t grandmomLabel = momP->GetMother();
1434  Int_t grandmomPDG = -1;
1435  AliVParticle * grandmomP = 0x0;
1436 
1437  while (grandmomLabel >=0 )
1438  {
1439  grandmomP = mcevent->GetTrack(grandmomLabel);
1440  grandmomPDG = grandmomP->PdgCode();
1441 
1442  if(grandmomPDG==pdg) break;
1443 
1444  grandmomLabel = grandmomP->GetMother();
1445  }
1446 
1447  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1448  {
1449  AliVParticle * d1 = mcevent->GetTrack(grandmomP->GetDaughterLabel(0));
1450  AliVParticle * d2 = mcevent->GetTrack(grandmomP->GetDaughterLabel(1));
1451 
1452  if(d1->PdgCode() == 22 && d1->PdgCode() == 22)
1453  {
1454  asym = (d1->E()-d2->E())/grandmomP->E();
1455  fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1456  fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1457  angle = fDaughMom.Angle(fDaughMom2.Vect());
1458  }
1459  }
1460  else
1461  {
1462  ok = kFALSE;
1463  AliInfo(Form("Mother with PDG %d, not found!",pdg));
1464  return;
1465  }
1466 
1467  ok = kTRUE;
1468 }
1469 
1470 //_________________________________________________________________________________________________
1472 //_________________________________________________________________________________________________
1473 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliMCEvent* mcevent, Bool_t & ok)
1474 {
1475  if ( !mcevent )
1476  {
1477  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1478 
1479  ok=kFALSE;
1480  return -1;
1481  }
1482 
1483  Int_t nprimaries = mcevent->GetNumberOfTracks();
1484  if ( label < 0 || label >= nprimaries )
1485  {
1486  ok = kFALSE;
1487  return -1;
1488  }
1489 
1490  AliVParticle * momP = mcevent->GetTrack(label);
1491 
1492  ok = kTRUE;
1493 
1494  return momP->GetNDaughters();
1495 }
1496 
1497 //_________________________________________________________________________________
1502 //_________________________________________________________________________________
1504  Int_t mctag, Int_t mesonLabel,
1505  AliMCEvent* mcevent,
1506  Int_t *overpdg, Int_t *overlabel)
1507 {
1508  Int_t ancPDG = 0, ancStatus = -1;
1509  TVector3 prodVertex;
1510  Int_t ancLabel = 0;
1511  Int_t noverlaps = 0;
1512  Bool_t ok = kFALSE;
1513 
1514  for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1515  {
1516  ancLabel = CheckCommonAncestor(label[0],label[ilab],mcevent,ancPDG,ancStatus,fMotherMom,prodVertex);
1517 
1518  //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1519  // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1520 
1521  Bool_t overlap = kFALSE;
1522 
1523  if ( ancLabel < 0 )
1524  {
1525  overlap = kTRUE;
1526  //printf("\t \t \t No Label = %d\n",ancLabel);
1527  }
1528  else if ( ( ancPDG==111 || ancPDG==221 ) &&
1529  ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta) ) &&
1530  ( (mesonLabel != ancLabel) && mesonLabel >=0 ) ) // in case the label is not provided check it is larger than 0
1531  {
1532  //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1533  overlap = kTRUE;
1534  }
1535  else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1536  {
1537  //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1538  overlap = kTRUE ;
1539  }
1540 
1541  if( !overlap ) continue ;
1542 
1543  // We have at least one overlap
1544 
1545  //printf("Overlap!!!!!!!!!!!!!!\n");
1546 
1547  noverlaps++;
1548 
1549  // What is the origin of the overlap?
1550  Bool_t mOK = 0, gOK = 0;
1551  Int_t mpdg = -999999, gpdg = -1;
1552  Int_t mstatus = -1, gstatus = -1;
1553  Int_t gLabel = -1, ggLabel = -1;
1554 
1555  GetMother (label[ilab],mcevent,mpdg,mstatus,mOK);
1556  fGMotherMom =
1557  GetGrandMother(label[ilab],mcevent,gpdg,gstatus,gOK, gLabel,ggLabel);
1558 
1559  //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
1560 
1561  if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
1562  ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
1563  gLabel >=0 )
1564  {
1565  Int_t labeltmp = gLabel;
1566  while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
1567  {
1568  mpdg=gpdg;
1569  fGMotherMom = GetGrandMother(labeltmp,mcevent,gpdg,gstatus,ok, gLabel,ggLabel);
1570  labeltmp=gLabel;
1571  }
1572  }
1573  overpdg [noverlaps-1] = mpdg;
1574  overlabel[noverlaps-1] = label[ilab];
1575  }
1576 
1577  return noverlaps ;
1578 }
1579 
1580 //________________________________________________________
1590 //________________________________________________________
1591 AliGenPythiaEventHeader * AliMCAnalysisUtils::CheckAndGetPythiaEventHeader
1592 (AliMCEvent* mcevent, TString selectHeaderName)
1593 {
1594  Int_t eventN = -1;
1595  if ( (AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler() )
1596  eventN = ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()->GetReadEntry());
1597 
1598  if ( fCurrentEvent != eventN )
1599  {
1600  fCurrentEvent = eventN ;
1601 
1602  //GetJets(mcevent, kTRUE);
1603 
1604  fPyGenHead = GetPythiaEventHeader(mcevent, selectHeaderName,
1607 
1608  if ( fPyVersion == 6 )
1609  {
1610  fMaxPartonicParent = 8 ;
1611  fMinPartonicParent = 5 ;
1612  }
1613  else if ( fPyVersion == 8 )
1614  {
1615  fMaxPartonicParent = 6 ;
1616  fMinPartonicParent = 3 ;
1617  }
1618 
1619  AliDebug(1,Form("For event %d, Pythia v%d name <%s>, process %d <%s>, first generated particle %d,"
1620  " photon partonic parent range [%d,%d]",
1622  fPyProcess, fPyProcessName.Data(),
1624  }
1625 
1626  return fPyGenHead ;
1627 }
1628 
1629 //________________________________________________________
1641 //________________________________________________________
1642 AliGenPythiaEventHeader * AliMCAnalysisUtils::GetPythiaEventHeader
1643 (AliMCEvent* mcevent, TString selectHeaderName,
1644  TString & genName, TString & processName,
1645  Int_t & process, Int_t & firstParticle,
1646  Int_t & pythiaVersion)
1647 {
1648  AliGenPythiaEventHeader * pyGenHead = 0;
1649  firstParticle = 0 ;
1650  genName = "";
1651  processName = "";
1652  pythiaVersion = 0 ;
1653  Int_t ngen = 0 ;
1654 
1655  TList * genlist = mcevent->GetCocktailList();
1656  if ( !genlist )
1657  {
1658  if ( !mcevent->GenEventHeader() ) return pyGenHead;
1659  else ngen = 1 ;
1660  }
1661  else
1662  {
1663  ngen = genlist->GetEntries();
1664  //printf("Cocktail ngen %d\n",ngen);
1665 
1666  if ( ngen < 1 )
1667  {
1668  if ( !mcevent->GenEventHeader() )
1669  return pyGenHead;
1670  else
1671  ngen = 1 ;
1672  }
1673  }
1674 
1675 // printf("GetPythiaEventHeader() - N generators %d, cocktail list %p, gen header %p\n",
1676 // ngen,genlist,mcevent->GenEventHeader());
1677 
1678  if ( ngen == 1 )
1679  {
1680  //printf("GetPythiaEventHeader() - Generator class name %s\n",mcevent->GenEventHeader()->ClassName());
1681  if ( strcmp(mcevent->GenEventHeader()->ClassName(), "AliGenPythiaEventHeader") ) return pyGenHead ;
1682 
1683  pyGenHead = (AliGenPythiaEventHeader*) mcevent->GenEventHeader();
1684 
1685  genName = pyGenHead->GetName();
1686  //printf("GetPythiaEventHeader() - 1 generator: %s\n",genName.Data());
1687  }
1688  else if ( ngen > 1 )
1689  {
1690  // Get the first pythia header
1691  for(Int_t igen = 0; igen < ngen; igen++)
1692  {
1693  AliGenEventHeader * header = (AliGenEventHeader*) genlist->At(igen);
1694 // printf("GetPythiaEventHeader() - igen %d Generator class %s, name %s\n",
1695 // igen, header->ClassName(),header->GetName());
1696 
1697  if ( !header || strcmp(header->ClassName(), "AliGenPythiaEventHeader") )
1698  {
1699  // Count the number of particles before the first header
1700  firstParticle += header->NProduced();
1701 
1702  continue;
1703  }
1704 
1705  genName = header->GetName();
1706 
1707  // In case of cocktail generators, select only the generator with a particular name
1708  if ( selectHeaderName.Length() > 0 && !genName.Contains(selectHeaderName) )
1709  {
1710  firstParticle += header->NProduced();
1711  genName = "";
1712 
1713  continue;
1714  }
1715 
1716  //printf("\t igen %d %s\n",igen,genName.Data());
1717  pyGenHead = (AliGenPythiaEventHeader*) genlist->At(igen);
1718 
1719  break;
1720  } // for igen
1721 
1722 // printf("GetPythiaEventHeader() - N generators %d, pythia pointer %p, %s, first particle %d\n",
1723 // genlist->GetEntries(), pyGenHead, genName.Data(), firstParticle);
1724  } // cocktail
1725 
1726  if ( !pyGenHead ) return pyGenHead;
1727 
1728  // Try to identify what pythia version was used
1729  // and what event type, a bit dangerous, only valid for gamma-jet and jet-jet events
1730  //
1731  process = pyGenHead->ProcessType();
1732 
1733  if ( genName.Contains("Pythia6") ) pythiaVersion = 6;
1734  else if ( genName.Contains("Pythia8") ) pythiaVersion = 8;
1735 
1736  if ( process == 14 || process == 18 || process == 29 )
1737  {
1738  if( !pythiaVersion ) pythiaVersion = 6;
1739  processName = "Gamma-Jet";
1740  }
1741  else if ( process == 11 || process == 12 ||
1742  process == 13 || process == 28 ||
1743  process == 53 || process == 68 || process == 96 )
1744  {
1745  if( !pythiaVersion ) pythiaVersion = 6;
1746  processName = "Jet-Jet";
1747  }
1748  else if (process >= 201 && process <= 205)
1749  {
1750  if( !pythiaVersion ) pythiaVersion = 8;
1751  processName = "Gamma-Jet";
1752  }
1753  else if ( (process >= 111 && process <= 116) || // jet-jet
1754  (process >= 121 && process <= 124) ) // jet-jet (b/c)
1755  {
1756  if( !pythiaVersion ) pythiaVersion = 8;
1757  processName = "Jet-Jet";
1758  }
1759 
1760 // printf("GetPythiaEventHeader() - Pythia version %d, process %d, process name %s\n",
1761 // pythiaVersion, process,processName.Data());
1762 
1763  // TString nameGen = "";
1764  // mcevent->GetCocktailGenerator(firstParticle-1,genName);
1765  // printf("\t -1: %d, %s\n",firstParticle-1,genName.Data());
1766  // mcevent->GetCocktailGenerator(firstParticle,genName);
1767  // printf("\t 0: %d, %s\n",firstParticle,genName.Data());
1768  // mcevent->GetCocktailGenerator(firstParticle+1,genName);
1769  // printf("\t +1: %d, %s\n",firstParticle+1,genName.Data());
1770 
1771  return pyGenHead;
1772 }
1773 
1774 //________________________________________________________
1776 //________________________________________________________
1778 {
1779  if(! opt)
1780  return;
1781 
1782  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1783 
1784  printf("Debug level = %d\n",fDebug);
1785  printf("MC Generator = %s\n",fMCGeneratorString.Data());
1786  printf(" \n");
1787 }
1788 
1789 //________________________________________________________
1797 //________________________________________________________
1798 void AliMCAnalysisUtils::PrintAncestry(AliMCEvent* mcevent, Int_t label, Int_t nGenerMax) const
1799 {
1800  AliVParticle * primary = 0;
1801  Int_t index = label;
1802  Int_t gener = 0;
1803 
1804  AliInfo("*********** Start");
1805  while ( index > 0 )
1806  {
1807  primary = mcevent->GetTrack(index);
1808 
1809  if(!primary)
1810  {
1811  AliWarning("primary pointer not available!!");
1812  return;
1813  }
1814 
1815  Float_t eta = 0;
1816  // Protection against floating point exception
1817  if ( primary->E() == TMath::Abs(primary->Pz()) ||
1818  (primary->E() - primary->Pz()) < 1e-3 ||
1819  (primary->E() + primary->Pz()) < 0 )
1820  eta = -999;
1821  else
1822  eta = primary->Eta();
1823 
1824  Int_t pdg = primary->PdgCode();
1825 
1826  printf("generation %d, label %d, %s, pdg %d, status %d, phys prim %d, E %2.2f, pT %2.2f, p(%2.2f,%2.2f,%2.2f) eta %2.2f, phi %2.2f"
1827  " mother %d, n daughters %d, d1 %d, d2 %d\n",
1828  gener,index,TDatabasePDG::Instance()->GetParticle(pdg)->GetName(),pdg,primary->MCStatusCode(),primary->IsPhysicalPrimary(),
1829  primary->E(),primary->Pt(),primary->Px(),primary->Py(),primary->Pz(),
1830  eta,primary->Phi()*TMath::RadToDeg(),
1831  primary->GetMother(),primary->GetNDaughters(),primary->GetDaughterLabel(0), primary->GetDaughterLabel(1));
1832 
1833  gener++;
1834  index = primary->GetMother();
1835  if ( nGenerMax < gener ) index = -1; // stop digging ancestry
1836  } // while
1837 
1838  AliInfo("*********** End");
1839 }
1840 
1841 
1842 //__________________________________________________
1844 //__________________________________________________
1846 {
1847  AliInfo
1848  (
1849  Form
1850  ("Tag %d: photon %d, conv %d, prompt %d, frag %d, isr %d,\n"
1851  " pi0 decay %d, eta decay %d, other decay %d, lost decay %d, in calo decay %d, pi0 %d, eta %d,\n"
1852  " electron %d, muon %d,pion %d, proton %d, neutron %d,\n"
1853  " kaon %d, a-proton %d, a-neutron %d, unk %d, bad %d",
1854  tag,
1855  CheckTagBit(tag,kMCPhoton),
1857  CheckTagBit(tag,kMCPrompt),
1859  CheckTagBit(tag,kMCISR),
1860  CheckTagBit(tag,kMCPi0Decay),
1861  CheckTagBit(tag,kMCEtaDecay),
1865  CheckTagBit(tag,kMCPi0),
1866  CheckTagBit(tag,kMCEta),
1867  CheckTagBit(tag,kMCElectron),
1868  CheckTagBit(tag,kMCMuon),
1869  CheckTagBit(tag,kMCPion),
1870  CheckTagBit(tag,kMCProton),
1872  CheckTagBit(tag,kMCKaon),
1873  CheckTagBit(tag,kMCAntiProton),
1875  CheckTagBit(tag,kMCUnknown),
1877  )
1878  );
1879 }
1880 
1881 //__________________________________________________
1883 //__________________________________________________
1885 {
1886  fMCGenerator = mcgen ;
1887  if (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
1888  else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
1889  else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
1890  else
1891  {
1892  fMCGeneratorString = "";
1894  }
1895 }
1896 
1897 //____________________________________________________
1899 //____________________________________________________
1901 {
1902  fMCGeneratorString = mcgen ;
1903 
1904  if (mcgen == "PYTHIA") fMCGenerator = kPythia;
1905  else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
1906  else if(mcgen == "HIJING") fMCGenerator = kHijing;
1907  else
1908  {
1910  fMCGeneratorString = "" ;
1911  }
1912 }
1913 
1914 
1915 
Int_t pdg
void SetMCGenerator(Int_t mcgen)
Set the generator type.
Int_t fCurrentEvent
Current Event number - GetJets()
double Double_t
Definition: External.C:58
TLorentzVector GetMother(Int_t label, const AliMCEvent *mcevent, Bool_t &ok)
TLorentzVector GetMotherWithPDG(Int_t label, Int_t pdg, const AliMCEvent *mcevent, Bool_t &ok, Int_t &momLabel)
TString fPyGenName
Pythia header assigned name.
AliGenPythiaEventHeader * CheckAndGetPythiaEventHeader(AliMCEvent *mcevent, TString selecHeaderName)
TLorentzVector fMotherMom
! particle momentum
TList * GetJets(AliMCEvent *mcevent, Bool_t check)
void PrintMCTag(Int_t tag) const
Print the assigned origins to this particle.
static AliGenPythiaEventHeader * GetPythiaEventHeader(AliMCEvent *mcevent, TString selecHeaderName, TString &genName, TString &processName, Int_t &process, Int_t &firstParticle, Int_t &pythiaVersion)
TLorentzVector GetDaughter(Int_t daughter, Int_t label, const AliMCEvent *mcevent, Int_t &pdg, Int_t &status, Bool_t &ok, Int_t &daugLabel, TVector3 &prodVertex)
void CheckOverlapped2GammaDecay(const Int_t *labels, const UShort_t *edepFrac, Int_t nlabels, Int_t mesonIndex, Float_t clusE, const AliMCEvent *mcevent, Int_t &tag)
Int_t fMinPartonicParent
Minimum label of partonic parent of direct photon.
void GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliMCEvent *mcevent, Float_t &asy, Float_t &angle, Bool_t &ok)
In case of an eta or pi0 decay into 2 photons, get the asymmetry in the energy of the photons...
TLorentzVector GetGrandMother(Int_t label, const AliMCEvent *mcevent, Int_t &pdg, Int_t &status, Bool_t &ok, Int_t &grandMomLabel, Int_t &greatMomLabel)
Int_t CheckOrigin(Int_t label, AliMCEvent *mcevent, TString selectHeaderName, Float_t clusE)
virtual ~AliMCAnalysisUtils()
Destructor.
void SetTagBit(Int_t &tag, UInt_t set) const
TLorentzVector fGMotherMom
! particle momentum
Int_t fDebug
Debug level.
TLorentzVector fDaughMom2
! particle momentum
AliGenPythiaEventHeader * fPyGenHead
! pythia event header of current event
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Int_t fPyFirstParticle
First Pythia generated particle in array.
TLorentzVector GetFirstMotherWithPDG(Int_t label, Int_t pdg, const AliMCEvent *mcevent, Bool_t &ok, Int_t &momLabel, Int_t &gparentlabel)
TString fPyProcessName
Pythia process name, Gamma-Jet or Jet-Jet.
void PrintAncestry(AliMCEvent *mcevent, Int_t label, Int_t nGenerMax=1000) const
Int_t fMCGenerator
MC generator used to generate data in simulation.
TLorentzVector fDaughMom
! particle momentum
TList * fJetsList
List of jets - GetJets()
Int_t fPyProcess
Pythia process code.
TString fMCGeneratorString
MC generator used to generate data in simulation.
void CheckLostDecayPair(const TObjArray *arrayCluster, Int_t iMom, Int_t iParent, const AliMCEvent *mcevent, Int_t &tag)
Check on AODs if the current decay photon has the second photon companion lost.
Int_t GetNDaughters(Int_t label, const AliMCEvent *mcevent, Bool_t &ok)
unsigned short UShort_t
Definition: External.C:28
Int_t fPyVersion
Pythia guessed version.
const char Option_t
Definition: External.C:48
Int_t CheckCommonAncestor(Int_t index1, Int_t index2, const AliMCEvent *mcevent, Int_t &ancPDG, Int_t &ancStatus, TLorentzVector &momentum, TVector3 &prodVertex)
AliMCAnalysisUtils()
Constructor.
bool Bool_t
Definition: External.C:53
Class with analysis utils for simulations.
TString meson
Bool_t CheckTagBit(Int_t tag, UInt_t test) const
Int_t GetNOverlaps(const Int_t *label, UInt_t nlabels, Int_t mctag, Int_t mesonLabel, AliMCEvent *mcevent, Int_t *overpdg, Int_t *overlabel)
Int_t fMaxPartonicParent
Minimum label of partonic parent of direct photon.