AliPhysics  ff1d528 (ff1d528)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "AliMCAnalysisUtils.h"
25 #include "AliMCEvent.h"
26 #include "AliGenPythiaEventHeader.h"
27 #include "AliVParticle.h"
28 #include "AliLog.h"
29 
33 
34 //________________________________________
36 //________________________________________
38 TObject(),
39 fCurrentEvent(-1),
40 fDebug(0),
41 fJetsList(new TList),
42 fMCGenerator(kPythia),
43 fMCGeneratorString("PYTHIA"),
44 fDaughMom(), fDaughMom2(),
45 fMotherMom(), fGMotherMom()
46 {}
47 
48 //_______________________________________
50 //_______________________________________
52 {
53  if (fJetsList)
54  {
55  fJetsList->Clear();
56  delete fJetsList ;
57  }
58 }
59 
60 //_____________________________________________________________________________________________
63 //_____________________________________________________________________________________________
65  const AliMCEvent* mcevent,
66  Int_t & ancPDG, Int_t & ancStatus,
67  TLorentzVector & momentum, TVector3 & prodVertex)
68 {
69  Int_t label1[100];
70  Int_t label2[100];
71  label1[0]= index1;
72  label2[0]= index2;
73  Int_t counter1 = 0;
74  Int_t counter2 = 0;
75 
76  if(label1[0]==label2[0])
77  {
78  //printf("AliMCAnalysisUtils::CheckCommonAncestor() - Already the same label: %d\n",label1[0]);
79  counter1=1;
80  counter2=1;
81  }
82  else
83  {
84  Int_t label=label1[0];
85  if(label < 0) return -1;
86 
87  while(label > -1 && counter1 < 99)
88  {
89  counter1++;
90  AliVParticle * mom = mcevent->GetTrack(label);
91  if(mom)
92  {
93  label = mom->GetMother() ;
94  label1[counter1]=label;
95  }
96  //printf("\t counter %d, label %d\n", counter1,label);
97  }
98 
99  //printf("Org label2=%d,\n",label2[0]);
100  label=label2[0];
101  if(label < 0) return -1;
102 
103  while(label > -1 && counter2 < 99)
104  {
105  counter2++;
106  AliVParticle * mom = mcevent->GetTrack(label);
107  if(mom)
108  {
109  label = mom->GetMother() ;
110  label2[counter2]=label;
111  }
112  //printf("\t counter %d, label %d\n", counter2,label);
113  }
114  }//First labels not the same
115 
116 // if((counter1==99 || counter2==99) && fDebug >=0)
117 // printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
118  //printf("CheckAncestor:\n");
119 
120  Int_t commonparents = 0;
121  Int_t ancLabel = -1;
122  //printf("counters %d %d \n",counter1, counter2);
123  for (Int_t c1 = 0; c1 < counter1; c1++)
124  {
125  for (Int_t c2 = 0; c2 < counter2; c2++)
126  {
127  if(label1[c1]==label2[c2] && label1[c1]>-1)
128  {
129  ancLabel = label1[c1];
130  commonparents++;
131 
132  AliVParticle * mom = mcevent->GetTrack(label1[c1]);
133 
134  if (mom)
135  {
136  ancPDG = mom->PdgCode();
137  ancStatus = mom->MCStatusCode();
138  momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
139  prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
140  }
141 
142  //First ancestor found, end the loops
143  counter1=0;
144  counter2=0;
145  }//Ancestor found
146  }//second cluster loop
147  }//first cluster loop
148 
149  if(ancLabel < 0)
150  {
151  ancPDG = -10000;
152  ancStatus = -10000;
153  momentum.SetXYZT(0,0,0,0);
154  prodVertex.SetXYZ(-10,-10,-10);
155  }
156 
157  return ancLabel;
158 }
159 
160 //____________________________________________________________________________________________________
166 //_____________________________________________________________________________________________________
167 Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, const AliMCEvent* mcevent)
168 {
169  Int_t labels[] = { label };
170 
171  return CheckOrigin(labels, 1, mcevent);
172 }
173 
174 //__________________________________________________________________________________________
189 //__________________________________________________________________________________________
191  const AliMCEvent* mcevent, const TObjArray* arrayCluster)
192 {
193  if( nlabels <= 0 )
194  {
195  AliWarning("No MC labels available, please check!!!");
196  return kMCBadLabel;
197  }
198 
199  if ( !mcevent )
200  {
201  AliDebug(1,"MCEvent is not available, check analysis settings in configuration file, do nothing with MC!!");
202  return -1;
203  }
204 
205  Int_t tag = 0;
206  Int_t nprimaries = mcevent->GetNumberOfTracks();
207 
208  // Bad label
209  if ( labels[0] < 0 || labels[0] >= nprimaries )
210  {
211  if(labels[0] < 0)
212  AliWarning(Form("*** bad label ***: label %d", labels[0]));
213 
214  if(labels[0] >= nprimaries)
215  AliWarning(Form("*** large label ***: label %d, n tracks %d", labels[0], nprimaries));
216 
217  SetTagBit(tag,kMCBadLabel);
218 
219  return tag;
220  } // Bad label
221 
222  // Most significant particle contributing to the cluster
223  Int_t label=labels[0];
224 
225  // Mother
226  AliVParticle * mom = mcevent->GetTrack(label);
227  Int_t iMom = label;
228  Int_t mPdgSign = mom->PdgCode();
229  Int_t mPdg = TMath::Abs(mPdgSign);
230  Int_t mStatus = mom->MCStatusCode() ;
231  Int_t iParent = mom->GetMother() ;
232 
233  //if(label < 8 && fMCGenerator != kBoxLike) AliDebug(1,Form("Mother is parton %d\n",iParent));
234 
235  //GrandParent
236  AliVParticle * parent = NULL ;
237  Int_t pPdg =-1;
238  Int_t pStatus =-1;
239  if(iParent >= 0)
240  {
241  parent = mcevent->GetTrack(iParent);
242  pPdg = TMath::Abs(parent->PdgCode());
243  pStatus = parent->MCStatusCode();
244  }
245  else AliDebug(1,Form("Parent with label %d",iParent));
246 
247  AliDebug(2,"Cluster most contributing mother and its parent:");
248  AliDebug(2,Form("\t Mother label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d",
249  iMom , mPdg, mStatus, mom->IsPrimary() , mom->IsPhysicalPrimary()));
250  AliDebug(2,Form("\t Parent label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d",
251  iParent, pPdg, pStatus, parent?parent->IsPrimary():-1, parent?parent->IsPhysicalPrimary():-1));
252 
253  //Check if mother is converted, if not, get the first non converted mother
254  if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus==0)
255  {
257 
258  // Check if the mother is photon or electron with status not stable
259  while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
260  {
261  // Mother
262  iMom = mom->GetMother();
263 
264  if(iMom < 0)
265  {
266  AliInfo(Form("pdg = %d, mother = %d, skip",pPdg,iMom));
267  break;
268  }
269 
270  mom = mcevent->GetTrack(iMom);
271  mPdgSign = mom->PdgCode();
272  mPdg = TMath::Abs(mPdgSign);
273  mStatus = mom->MCStatusCode() ;
274  iParent = mom->GetMother() ;
275  //if(label < 8 ) AliDebug(1, Form("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent));
276 
277  // GrandParent
278  if(iParent >= 0 && parent)
279  {
280  parent = mcevent->GetTrack(iParent);
281  pPdg = TMath::Abs(parent->PdgCode());
282  pStatus = parent->MCStatusCode();
283  }
284  // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
285  // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
286 
287  }//while
288 
289  AliDebug(2,"Converted photon/electron:");
290  AliDebug(2,Form("\t Mother label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d"
291  ,iMom , mPdg, mStatus, mom->IsPrimary() , mom->IsPhysicalPrimary()));
292  AliDebug(2,Form("\t Parent label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d"
293  ,iParent, pPdg, pStatus, parent?parent->IsPrimary():-1, parent?parent->IsPhysicalPrimary():-1));
294 
295  } // mother and parent are electron or photon and have status 0 and parent is photon or electron
296  else if((mPdg == 22 || mPdg == 11) && mStatus==0)
297  {
298  // Still a conversion but only one electron/photon generated. Just from hadrons
299  if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
300  pPdg == 2212 || pPdg == 130 || pPdg == 13 )
301  {
303  iMom = mom->GetMother();
304 
305  if(iMom < 0)
306  {
307  AliInfo(Form("pdg = %d, mother = %d, skip",pPdg,iMom));
308  }
309  else
310  {
311  mom = mcevent->GetTrack(iMom);
312  mPdgSign = mom->PdgCode();
313  mPdg = TMath::Abs(mPdgSign);
314  mStatus = mom->MCStatusCode() ;
315 
316  AliDebug(2,"Converted hadron:");
317  AliDebug(2,Form("\t Mother label %d, pdg %d, status %d, Primary? %d, Physical Primary? %d",
318  iMom, mPdg, mStatus, mom->IsPrimary(), mom->IsPhysicalPrimary()));
319  }
320  } // hadron converted
321 
322  //Comment for next lines, we do not check the parent of the hadron for the moment.
323  //iParent = mom->GetMother() ;
324  //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
325 
326  //GrandParent
327  //if(iParent >= 0){
328  // parent = mcevent->GetTrack(iParent);
329  // pPdg = TMath::Abs(parent->PdgCode());
330  //}
331  }
332 
333  //printf("Final mother mPDG %d\n",mPdg);
334 
335  // conversion into electrons/photons checked
336 
337  //first check for typical charged particles
338  if (mPdg == 13) SetTagBit(tag,kMCMuon);
339  else if(mPdg == 211) SetTagBit(tag,kMCPion);
340  else if(mPdg == 321) SetTagBit(tag,kMCKaon);
341  else if(mPdgSign == 2212) SetTagBit(tag,kMCProton);
342  else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron);
343  else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
344  else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
345 
346  //check for pi0 and eta (shouldn't happen unless their decays were turned off)
347  else if(mPdg == 111)
348  {
349  SetTagBit(tag,kMCPi0Decay);
350 
351  AliDebug(2,"First mother is directly pi0, not decayed by generator");
352 
353  CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcevent, tag); //set to kMCPi0 if 2 gammas in same cluster
354  }
355  else if(mPdg == 221)
356  {
357  SetTagBit(tag,kMCEtaDecay);
358 
359  AliDebug(2,"First mother is directly eta, not decayed by generator");
360 
361  CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcevent, tag); //set to kMCEta if 2 gammas in same cluster
362  }
363  //Photons
364  else if(mPdg == 22)
365  {
366  SetTagBit(tag,kMCPhoton);
367 
368  if(pPdg == 111)
369  {
370  SetTagBit(tag,kMCPi0Decay);
371 
372  AliDebug(2,"Generator pi0 decay photon");
373 
374  // Set to kMCPi0 if 2 gammas in same cluster
375  CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcevent, tag);
376 
377  // In case it did not merge, check if the decay companion is lost
379  {
380  CheckLostDecayPair(arrayCluster,iMom, iParent, mcevent, tag);
381  }
382 
383  //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
384  }
385  else if (pPdg == 221)
386  {
387  SetTagBit(tag, kMCEtaDecay);
388 
389  AliDebug(2,"Generator eta decay photon");
390 
391  // Set to kMCEta if 2 gammas in same cluster
392  CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcevent, tag);
393 
394  // In case it did not merge, check if the decay companion is lost
396  CheckLostDecayPair(arrayCluster,iMom, iParent, mcevent, tag);
397  }
398  else if( mom->IsPhysicalPrimary() && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle
399  {
400  if(iParent < 8 && iParent > 5 )
401  {
402  //outgoing partons
403  if(pPdg == 22) SetTagBit(tag,kMCPrompt);
404  else SetTagBit(tag,kMCFragmentation);
405  }//Outgoing partons
406  else if( iParent <= 5 && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) )
407  {
408  SetTagBit(tag, kMCISR); //Initial state radiation
409  }
410  else SetTagBit(tag,kMCUnknown);
411  }//Physical primary
412  else SetTagBit(tag,kMCOtherDecay);
413 
414 // //Old Herwig selection for ESDs, maybe should be considered.
415 // if(fMCGenerator == kHerwig)
416 // {
417 // if(pStatus < 197)
418 // {//Not decay
419 // while(1)
420 // {
421 // if(parent)
422 // {
423 // if(parent->GetMother()<=5) break;
424 // iParent = parent->GetMother();
425 // parent = mcevent->GetTrack(iParent);
426 // pStatus = parent->MCStatusCode();
427 // pPdg = TMath::Abs(parent->PdgCode());
428 // } else break;
429 // }//Look for the parton
430 //
431 // if(iParent < 8 && iParent > 5)
432 // {
433 // if(pPdg == 22) SetTagBit(tag,kMCPrompt);
434 // else SetTagBit(tag,kMCFragmentation);
435 // }
436 // else SetTagBit(tag,kMCISR);//Initial state radiation
437 // }//Not decay
438 // else SetTagBit(tag,kMCUnknown);
439 // }//HERWIG
440 
441  } // Mother Photon
442 
443  // Electron check. Where did that electron come from?
444  else if(mPdg == 11)
445  {
446  //electron
447  if(pPdg == 11 && parent)
448  {
449  Int_t iGrandma = parent->GetMother();
450  if(iGrandma >= 0)
451  {
452  AliVParticle * gma = mcevent->GetTrack(iGrandma);
453  Int_t gPdg = TMath::Abs(gma->PdgCode());
454 
455  if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
456  else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
457  }
458  }
459 
460  SetTagBit(tag,kMCElectron);
461 
462  AliDebug(1,"Checking ancestors of electrons");
463 
464  if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay
465  else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay
466  else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
467  else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
468  {
469  //c-hadron decay check
470  if(parent)
471  {
472  Int_t iGrandma = parent->GetMother();
473  if(iGrandma >= 0)
474  {
475  AliVParticle * gma = mcevent->GetTrack(iGrandma); //charm's mother
476  Int_t gPdg = TMath::Abs(gma->PdgCode());
477  if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
478  else SetTagBit(tag,kMCEFromC); //c-hadron decay
479  }
480  else SetTagBit(tag,kMCEFromC); //c-hadron decay
481  }//parent
482  } else
483  { //prompt or other decay
484 
485  TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
486  TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
487 
488  AliDebug(1,Form("Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)",foo->GetName(), pPdg,foo1->GetName(),mPdg));
489 
490  if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
491  else SetTagBit(tag,kMCOtherDecay);
492  }
493  }// electron check
494  // cluster was made by something else
495  else
496  {
497  AliDebug(1,Form("\t Setting kMCUnknown for cluster with pdg = %d, Parent pdg = %d",mPdg,pPdg));
498  SetTagBit(tag,kMCUnknown);
499  }
500 
501  return tag;
502 }
503 
504 //_________________________________________________________________________________________
507 //_________________________________________________________________________________________
509  Int_t mesonIndex, const AliMCEvent* mcevent,
510  Int_t & tag)
511 {
512  if(labels[0] < 0 || labels[0] > mcevent->GetNumberOfTracks() || nlabels <= 1)
513  {
514  AliDebug(2,Form("Exit : label[0] %d, n primaries %d, nlabels %d",labels[0],mcevent->GetNumberOfTracks(), nlabels));
515  return;
516  }
517 
518  AliVParticle * meson = mcevent->GetTrack(mesonIndex);
519  Int_t mesonPdg = meson->PdgCode();
520  if(mesonPdg != 111 && mesonPdg != 221)
521  {
522  AliWarning(Form("Wrong pi0/eta PDG : %d",mesonPdg));
523  return;
524  }
525 
526  AliDebug(2,Form("pdg %d, label %d, ndaughters %d", mesonPdg, mesonIndex, meson->GetNDaughters()));
527 
528  // Get the daughters
529  if(meson->GetNDaughters() != 2)
530  {
531 // if(meson->GetNDaughters()>2)
532 // {
533 // printf("xx More than 2 daughters <%d> for pdg %d with E %2.2f, pT %2.2f:\n",
534 // meson->GetNDaughters(),meson->PdgCode(),meson->E(),meson->Pt());
535 // for(Int_t idaugh = 0; idaugh < meson->GetNDaughters(); idaugh++)
536 // {
537 // if(meson->GetDaughterLabel(idaugh) < 0)
538 // {
539 // printf("\t no daughther with index %d:\n",idaugh);
540 // continue;
541 // }
542 // AliVParticle * daugh = mcevent->GetTrack(meson->GetDaughterLabel(idaugh));
543 // printf("\t daughther pdg %d, E %2.2f, pT %2.2f:\n",daugh->PdgCode(),daugh->E(),daugh->Pt());
544 // }
545 // }
546 
547  AliDebug(2,Form("Not overlapped. Number of daughters is %d, not 2",meson->GetNDaughters()));
548  return;
549  }
550 
551  Int_t iPhoton0 = meson->GetDaughterLabel(0);
552  Int_t iPhoton1 = meson->GetDaughterLabel(1);
553 
554  //if((iPhoton0 == -1) || (iPhoton1 == -1))
555  //{
556  // if(fDebug > 2)
557  // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
558  // return;
559  //}
560 
561  AliVParticle *photon0 = mcevent->GetTrack(iPhoton0);
562  AliVParticle *photon1 = mcevent->GetTrack(iPhoton1);
563 
564  // Check if both daughters are photons
565  if(photon0->PdgCode() != 22 && photon1->PdgCode()!=22)
566  {
567  AliWarning(Form("Not overlapped. PDG: daughter 1 = %d, of daughter 2 = %d",photon0->PdgCode(),photon1->PdgCode()));
568  return;
569  }
570 
571  AliDebug(2,Form("Daughter labels : photon0 = %d, photon1 = %d",iPhoton0,iPhoton1));
572 
573  // Check if both photons contribute to the cluster
574  Bool_t okPhoton0 = kFALSE;
575  Bool_t okPhoton1 = kFALSE;
576 
577  AliDebug(3,"Labels loop:");
578 
579  Bool_t conversion = kFALSE;
580 
581 // TLorentzVector d0, d1;
582 // photon0->Momentum(d0);
583 // photon1->Momentum(d1);
584 // Float_t openAngle = d0.Angle(d1.Vect());
585 
586  for(Int_t i = 0; i < nlabels; i++)
587  {
588  AliDebug(3, Form("\t label %d/%d: %d, ok? %d, %d", i, nlabels, labels[i], okPhoton0, okPhoton1));
589 
590  if ( labels[i] < 0 ) continue;
591 
592  // If we already found both, break the loop
593  if(okPhoton0 && okPhoton1) break;
594 
595  Int_t index = labels[i];
596  if (iPhoton0 == index)
597  {
598  okPhoton0 = kTRUE;
599  continue;
600  }
601  else if (iPhoton1 == index)
602  {
603  okPhoton1 = kTRUE;
604  continue;
605  }
606 
607  // Trace back the mother in case it was a conversion
608  if(index >= mcevent->GetNumberOfTracks())
609  {
610  AliWarning(Form("Particle index %d larger than size of list %d!!",index,mcevent->GetNumberOfTracks()));
611  continue;
612  }
613 
614  AliVParticle * daught = mcevent->GetTrack(index);
615  Int_t tmpindex = daught->GetMother();
616  AliDebug(3,Form("Conversion? : mother %d",tmpindex));
617 
618  while(tmpindex>=0)
619  {
620  // MC particle of interest is the mother
621  AliDebug(3,Form("\t parent index %d",tmpindex));
622  daught = mcevent->GetTrack(tmpindex);
623  //printf("tmpindex %d\n",tmpindex);
624  if (iPhoton0 == tmpindex)
625  {
626  conversion = kTRUE;
627  okPhoton0 = kTRUE;
628  break;
629  }
630  else if (iPhoton1 == tmpindex)
631  {
632  conversion = kTRUE;
633  okPhoton1 = kTRUE;
634  break;
635  }
636 
637  tmpindex = daught->GetMother();
638 
639  }//While to check if pi0/eta daughter was one of these contributors to the cluster
640 
641  //if(i == 0 && (!okPhoton0 && !okPhoton1)) AliDebug(1,"Something happens, first label should be from a photon decay!");
642 
643  }//loop on list of labels
644 
645  //If both photons contribute tag as the corresponding meson.
646  if(okPhoton0 && okPhoton1)
647  {
648  AliDebug(2,Form("%s OVERLAPPED DECAY",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()));
649 
650  if(!CheckTagBit(tag,kMCConversion) && conversion)
651  {
652  AliDebug(2,"Second decay photon produced a conversion");
653  SetTagBit(tag,kMCConversion) ;
654  }
655 
656  if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
657  else SetTagBit(tag,kMCEta);
658 
659 // printf("\t Meson E %2.2f, pT %2.2f, Overlapped open angle %2.2f deg, cells %2.1f, conversion %d\n",
660 // meson->E(), meson->Pt(),openAngle*TMath::RadToDeg(),openAngle/0.015,conversion);
661 //
662 // Float_t phi = photon0->Phi()*TMath::RadToDeg();
663 // if(phi < 0) phi+=360;
664 // printf("\t \t daughther1: label %d, E %2.2f, pT %2.2f, eta%2.2f, phi %2.2f, status %d, phys prim %d:\n",
665 // iPhoton0, photon0->E(),photon0->Pt(),
666 // photon0->Eta(),phi,
667 // photon0->MCStatusCode(),photon0->IsPhysicalPrimary());
668 //
669 // phi = photon1->Phi()*TMath::RadToDeg();
670 // if(phi < 0) phi+=360;
671 // printf("\t \t daughther2: label %d, E %2.2f, pT %2.2f, eta%2.2f, phi %2.2f, status %d, phys prim %d:\n",
672 // iPhoton1, photon1->E(),photon1->Pt(),
673 // photon1->Eta(),phi,
674 // photon1->MCStatusCode(),photon1->IsPhysicalPrimary());
675  }
676 }
677 
678 //________________________________________________________________________________________________________
680 //________________________________________________________________________________________________________
681 void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster, Int_t iMom, Int_t iParent,
682  const AliMCEvent* mcevent, Int_t & tag)
683 {
684  if(!arrayCluster || iMom < 0 || iParent < 0|| !mcevent) return;
685 
686  AliVParticle * parent = mcevent->GetTrack(iParent);
687 
688  //printf("*** Check label %d with parent %d\n",iMom, iParent);
689 
690  if(parent->GetNDaughters()!=2)
691  {
693  //printf("\t ndaugh = %d\n",parent->GetNDaughters());
694  return ;
695  }
696 
697  Int_t pairLabel = -1;
698  if ( iMom != parent->GetDaughterLabel(0) ) pairLabel = parent->GetDaughterLabel(0);
699  else if( iMom != parent->GetDaughterLabel(1) ) pairLabel = parent->GetDaughterLabel(1);
700 
701  if(pairLabel<0)
702  {
703  //printf("\t pair Label not found = %d\n",pairLabel);
705  return ;
706  }
707 
708  //printf("\t *** find pair %d\n",pairLabel);
709 
710  for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
711  {
712  AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
713  //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
714 
715  for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
716  {
717  Int_t label = cluster->GetLabels()[ilab];
718 
719  //printf("\t \t label %d\n",label);
720 
721  if ( label==pairLabel )
722  {
723  //printf("\t \t Pair found\n");
725  return ;
726  }
727  else if ( label== iParent || label== iMom )
728  {
729  //printf("\t \t skip\n");
730  continue;
731  }
732  else // check the ancestry
733  {
734  AliVParticle * mother = mcevent->GetTrack(label);
735 
736  if ( !mother )
737  {
738  AliInfo(Form("MC Mother not available for label %d",label));
739  continue;
740  }
741 
742  Int_t momPDG = TMath::Abs(mother->PdgCode());
743  if ( momPDG!=11 && momPDG!=22 ) continue;
744 
745  // Check if "mother" of entity is converted, if not, get the first non converted mother
746  Int_t iParentClus = mother->GetMother();
747  if(iParentClus < 0) continue;
748 
749  AliVParticle * parentClus = mcevent->GetTrack(iParentClus);
750  if(!parentClus) continue;
751 
752  Int_t parentClusPDG = TMath::Abs(parentClus->PdgCode());
753  Int_t parentClusStatus = parentClus->MCStatusCode();
754 
755  if ( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0 )
756  {
757  //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
758  continue;
759  }
760 
761  //printf("\t \t Conversion\n");
762 
763  // Check if the mother is photon or electron with status not stable
764  while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
765  {
766  //New Mother
767  label = iParentClus;
768  momPDG = parentClusPDG;
769 
770  iParentClus = parentClus->GetMother();
771  if ( iParentClus < 0 ) break;
772 
773  parentClus = mcevent->GetTrack(iParentClus);
774  if ( !parentClus ) break;
775 
776  parentClusPDG = TMath::Abs(parentClus->PdgCode());
777  parentClusStatus = parentClus->MCStatusCode() ;
778  }//while
779 
780  if ( (momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel) )
781  {
783  //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
784  return ;
785  }
786  else
787  {
788  //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
789  continue;
790  }
791 
792  }
793  }
794  } // cluster loop
795 
797 }
798 
799 //_____________________________________________________________________
801 //_____________________________________________________________________
802 TList * AliMCAnalysisUtils::GetJets(AliMCEvent* mcevent, AliGenEventHeader * mcheader, Int_t eventNumber)
803 {
804  if(fCurrentEvent==eventNumber) return fJetsList ;
805 
806  fCurrentEvent = eventNumber;
807 
808  if (fJetsList) fJetsList->Clear();
809  else fJetsList = new TList;
810 
811  Int_t nTriggerJets = 0;
812  Float_t tmpjet[]={0,0,0,0};
813 
814  //printf("Event %d %d\n",fCurrentEvent,iEvent);
815 
816  // Get outgoing partons
817  if(mcevent->GetNumberOfTracks() < 8) return fJetsList;
818 
819  AliVParticle * parton1 = mcevent->GetTrack(6);
820  AliVParticle * parton2 = mcevent->GetTrack(7);
821 
822  AliDebug(2,Form("Parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
823  parton1->GetName(),parton1->Pt(),parton1->E(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()));
824  AliDebug(2,Form("Parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
825  parton2->GetName(),parton2->Pt(),parton2->E(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()));
826 
827  // //Trace the jet from the mother parton
828  // Float_t pt = 0;
829  // Float_t pt1 = 0;
830  // Float_t pt2 = 0;
831  // Float_t e = 0;
832  // Float_t e1 = 0;
833  // Float_t e2 = 0;
834  // TParticle * tmptmp = new TParticle;
835  // for(Int_t i = 0; i< stack->GetNprimary(); i++){
836  // tmptmp = stack->Particle(i);
837 
838  // if(tmptmp->GetStatusCode() == 1){
839  // pt = tmptmp->Pt();
840  // e = tmptmp->Energy();
841  // Int_t imom = tmptmp->GetMother();
842  // Int_t imom1 = 0;
843  // //printf("1st imom %d\n",imom);
844  // while(imom > 5){
845  // imom1=imom;
846  // tmptmp = stack->Particle(imom);
847  // imom = tmptmp->GetMother();
848  // //printf("imom %d \n",imom);
849  // }
850  // //printf("Last imom %d %d\n",imom1, imom);
851  // if(imom1 == 6) {
852  // pt1+=pt;
853  // e1+=e;
854  // }
855  // else if (imom1 == 7){
856  // pt2+=pt;
857  // e2+=e; }
858  // }// status 1
859 
860  // }// for
861 
862  // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
863 
864  //Get the jet, different way for different generator
865  //PYTHIA
866  if(fMCGenerator == kPythia)
867  {
868  TParticle * jet = 0x0;
869  AliGenPythiaEventHeader* pygeh = (AliGenPythiaEventHeader*) mcheader;
870  nTriggerJets = pygeh->NTriggerJets();
871  AliDebug(2,Form("PythiaEventHeader: Njets: %d",nTriggerJets));
872 
873  for(Int_t i = 0; i< nTriggerJets; i++)
874  {
875  pygeh->TriggerJet(i, tmpjet);
876 
877  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
878 
879  // Assign an outgoing parton as mother
880  Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
881  Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
882 
883  if(phidiff1 > phidiff2) jet->SetFirstMother(7);
884  else jet->SetFirstMother(6);
885 
886  //jet->Print();
887  AliDebug(1,Form("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
888  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()));
889  fJetsList->Add(jet);
890  }
891  }//Pythia triggered jets
892  //HERWIG
893  else if (fMCGenerator == kHerwig)
894  {
895  Int_t pdg = -1;
896 
897  // Check parton 1
898  AliVParticle * tmp = parton1;
899  if(parton1->PdgCode()!=22)
900  {
901  while(pdg != 94)
902  {
903  if(tmp->GetFirstDaughter()==-1) return fJetsList;
904 
905  tmp = mcevent->GetTrack(tmp->GetFirstDaughter());
906  pdg = tmp->PdgCode();
907  }//while
908 
909  // Add found jet to list
910  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);
911 
912  jet1->SetFirstMother(6);
913 
914  fJetsList->Add(jet1);
915 
916  //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
917  //tmp = stack->Particle(tmp->GetFirstDaughter());
918  //tmp->Print();
919  //jet1->Print();
920  AliDebug(1,Form("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
921  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()));
922  }//not photon
923 
924  // Check parton 2
925  pdg = -1;
926  tmp = parton2;
927  if(parton2->PdgCode()!=22)
928  {
929  while(pdg != 94)
930  {
931  if(tmp->GetFirstDaughter()==-1) return fJetsList;
932 
933  tmp = mcevent->GetTrack(tmp->GetFirstDaughter());
934  pdg = tmp->PdgCode();
935  }//while
936 
937  // Add found jet to list
938  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);
939 
940  jet2->SetFirstMother(7);
941 
942  fJetsList->Add(jet2);
943 
944  //jet2->Print();
945  AliDebug(2,Form("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
946  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()));
947  //Int_t first = tmp->GetFirstDaughter();
948  //Int_t last = tmp->GetLastDaughter();
949  //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
950  // for(Int_t d = first ; d < last+1; d++){
951  // tmp = stack->Particle(d);
952  // if(i == tmp->GetMother())
953  // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
954  // d,tmp->GetMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
955  // }
956  //tmp->Print();
957  }//not photon
958  }//Herwig generated jets
959 
960  return fJetsList;
961 }
962 
963 //______________________________________________________________________________
965 //______________________________________________________________________________
966 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
967  const AliMCEvent* mcevent,
968  Int_t & pdg, Int_t & status,
969  Bool_t & ok, Int_t & daughlabel, TVector3 & prodVertex)
970 {
971  fDaughMom.SetPxPyPzE(0,0,0,0);
972 
973  if(!mcevent)
974  {
975  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
976 
977  ok = kFALSE;
978  return fDaughMom;
979  }
980 
981  Int_t nprimaries = mcevent->GetNumberOfTracks();
982  if(label < 0 || label >= nprimaries)
983  {
984  ok = kFALSE;
985  return fDaughMom;
986  }
987 
988  AliVParticle * momP = mcevent->GetTrack(label);
989  daughlabel = momP->GetDaughterLabel(idaugh);
990 
991  if(daughlabel < 0 || daughlabel >= nprimaries)
992  {
993  ok = kFALSE;
994  return fDaughMom;
995  }
996 
997  AliVParticle * daughP = mcevent->GetTrack(daughlabel);
998  fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
999  pdg = daughP->PdgCode();
1000  status = daughP->MCStatusCode();
1001  prodVertex.SetXYZ(daughP->Xv(),daughP->Yv(),daughP->Zv());
1002 
1003  ok = kTRUE;
1004 
1005  return fDaughMom;
1006 }
1007 
1008 //______________________________________________________________________________________________________
1010 //______________________________________________________________________________________________________
1011 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent, Bool_t & ok)
1012 {
1013  Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1014 
1015  return GetMother(label,mcevent,pdg,status, ok,momlabel);
1016 }
1017 
1018 //_________________________________________________________________________________________
1019 // \return the kinematics of the particle that generated the signal.
1020 //_________________________________________________________________________________________
1021 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent,
1022  Int_t & pdg, Int_t & status, Bool_t & ok)
1023 {
1024  Int_t momlabel = -1;
1025 
1026  return GetMother(label,mcevent,pdg,status, ok,momlabel);
1027 }
1028 
1029 //______________________________________________________________________________________________________
1031 //______________________________________________________________________________________________________
1032 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent,
1033  Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1034 {
1035  fMotherMom.SetPxPyPzE(0,0,0,0);
1036 
1037  if(!mcevent)
1038  {
1039  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1040 
1041  ok = kFALSE;
1042  return fMotherMom;
1043  }
1044 
1045  Int_t nprimaries = mcevent->GetNumberOfTracks();
1046  if(label < 0 || label >= nprimaries)
1047  {
1048  ok = kFALSE;
1049  return fMotherMom;
1050  }
1051 
1052  AliVParticle * momP = mcevent->GetTrack(label);
1053  fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1054  pdg = momP->PdgCode();
1055  status = momP->MCStatusCode();
1056  momlabel = momP->GetMother();
1057 
1058  ok = kTRUE;
1059 
1060  return fMotherMom;
1061 }
1062 
1063 //___________________________________________________________________________________
1065 //___________________________________________________________________________________
1067  const AliMCEvent* mcevent,
1068  Bool_t & ok, Int_t & momlabel)
1069 {
1070  fGMotherMom.SetPxPyPzE(0,0,0,0);
1071 
1072  if ( !mcevent )
1073  {
1074  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
1075 
1076  ok = kFALSE;
1077  return fGMotherMom;
1078  }
1079 
1080  Int_t nprimaries = mcevent->GetNumberOfTracks();
1081  if ( label < 0 || label >= nprimaries )
1082  {
1083  ok = kFALSE;
1084  return fGMotherMom;
1085  }
1086 
1087 
1088  AliVParticle * momP = mcevent->GetTrack(label);
1089 
1090  if(momP->PdgCode()==pdg)
1091  {
1092  AliDebug(2,"PDG of mother is already the one requested!");
1093  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1094 
1095  ok=kTRUE;
1096  return fGMotherMom;
1097  }
1098 
1099  Int_t grandmomLabel = momP->GetMother();
1100  Int_t grandmomPDG = -1;
1101  AliVParticle * grandmomP = 0x0;
1102 
1103  while (grandmomLabel >=0 )
1104  {
1105  grandmomP = mcevent->GetTrack(grandmomLabel);
1106  grandmomPDG = grandmomP->PdgCode();
1107  if(grandmomPDG==pdg)
1108  {
1109  //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1110  momlabel = grandmomLabel;
1111  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1112  break;
1113  }
1114 
1115  grandmomLabel = grandmomP->GetMother();
1116  }
1117 
1118  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found!",pdg));
1119 
1120  ok = kTRUE;
1121 
1122  return fGMotherMom;
1123 }
1124 
1125 //______________________________________________________________________________________________
1127 //______________________________________________________________________________________________
1128 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliMCEvent* mcevent,
1129  Int_t & pdg, Int_t & status, Bool_t & ok,
1130  Int_t & grandMomLabel, Int_t & greatMomLabel)
1131 {
1132  fGMotherMom.SetPxPyPzE(0,0,0,0);
1133 
1134  if ( !mcevent )
1135  {
1136  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1137 
1138  ok = kFALSE;
1139  return fGMotherMom;
1140  }
1141 
1142  Int_t nprimaries = mcevent->GetNumberOfTracks();
1143  if ( label < 0 || label >= nprimaries )
1144  {
1145  ok = kFALSE;
1146  return fGMotherMom;
1147  }
1148 
1149  AliVParticle * momP = mcevent->GetTrack(label);
1150 
1151  grandMomLabel = momP->GetMother();
1152 
1153  AliVParticle * grandmomP = 0x0;
1154 
1155  if(grandMomLabel >=0 )
1156  {
1157  grandmomP = mcevent->GetTrack(grandMomLabel);
1158  pdg = grandmomP->PdgCode();
1159  status = grandmomP->MCStatusCode();
1160 
1161  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1162  greatMomLabel = grandmomP->GetMother();
1163 
1164  }
1165 
1166  ok = kTRUE;
1167 
1168  return fGMotherMom;
1169 }
1170 
1171 //_______________________________________________________________________________________________________________
1173 //_______________________________________________________________________________________________________________
1175  Float_t & asym, Float_t & angle, Bool_t & ok)
1176 {
1177  if(!mcevent)
1178  {
1179  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1180 
1181  ok = kFALSE;
1182  return;
1183  }
1184 
1185  Int_t nprimaries = mcevent->GetNumberOfTracks();
1186  if ( label < 0 || label >= nprimaries )
1187  {
1188  ok = kFALSE;
1189  return ;
1190  }
1191 
1192  AliVParticle * momP = mcevent->GetTrack(label);
1193 
1194  Int_t grandmomLabel = momP->GetMother();
1195  Int_t grandmomPDG = -1;
1196  AliVParticle * grandmomP = 0x0;
1197 
1198  while (grandmomLabel >=0 )
1199  {
1200  grandmomP = mcevent->GetTrack(grandmomLabel);
1201  grandmomPDG = grandmomP->PdgCode();
1202 
1203  if(grandmomPDG==pdg) break;
1204 
1205  grandmomLabel = grandmomP->GetMother();
1206  }
1207 
1208  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1209  {
1210  AliVParticle * d1 = mcevent->GetTrack(grandmomP->GetDaughterLabel(0));
1211  AliVParticle * d2 = mcevent->GetTrack(grandmomP->GetDaughterLabel(1));
1212 
1213  if(d1->PdgCode() == 22 && d1->PdgCode() == 22)
1214  {
1215  asym = (d1->E()-d2->E())/grandmomP->E();
1216  fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1217  fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1218  angle = fDaughMom.Angle(fDaughMom2.Vect());
1219  }
1220  }
1221  else
1222  {
1223  ok = kFALSE;
1224  AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
1225  return;
1226  }
1227 
1228  ok = kTRUE;
1229 }
1230 
1231 //_________________________________________________________________________________________________
1233 //_________________________________________________________________________________________________
1234 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliMCEvent* mcevent, Bool_t & ok)
1235 {
1236  if ( !mcevent )
1237  {
1238  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1239 
1240  ok=kFALSE;
1241  return -1;
1242  }
1243 
1244  Int_t nprimaries = mcevent->GetNumberOfTracks();
1245  if ( label < 0 || label >= nprimaries )
1246  {
1247  ok = kFALSE;
1248  return -1;
1249  }
1250 
1251  AliVParticle * momP = mcevent->GetTrack(label);
1252 
1253  ok = kTRUE;
1254 
1255  return momP->GetNDaughters();
1256 }
1257 
1258 //_________________________________________________________________________________
1263 //_________________________________________________________________________________
1265  Int_t mctag, Int_t mesonLabel,
1266  AliMCEvent* mcevent,
1267  Int_t *overpdg, Int_t *overlabel)
1268 {
1269  Int_t ancPDG = 0, ancStatus = -1;
1270  TVector3 prodVertex;
1271  Int_t ancLabel = 0;
1272  Int_t noverlaps = 0;
1273  Bool_t ok = kFALSE;
1274 
1275  for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1276  {
1277  ancLabel = CheckCommonAncestor(label[0],label[ilab],mcevent,ancPDG,ancStatus,fMotherMom,prodVertex);
1278 
1279  //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1280  // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1281 
1282  Bool_t overlap = kFALSE;
1283 
1284  if ( ancLabel < 0 )
1285  {
1286  overlap = kTRUE;
1287  //printf("\t \t \t No Label = %d\n",ancLabel);
1288  }
1289  else if ( ( ancPDG==111 || ancPDG==221 ) &&
1290  ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta) ) &&
1291  ( (mesonLabel != ancLabel) && mesonLabel >=0 ) ) // in case the label is not provided check it is larger than 0
1292  {
1293  //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1294  overlap = kTRUE;
1295  }
1296  else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1297  {
1298  //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1299  overlap = kTRUE ;
1300  }
1301 
1302  if( !overlap ) continue ;
1303 
1304  // We have at least one overlap
1305 
1306  //printf("Overlap!!!!!!!!!!!!!!\n");
1307 
1308  noverlaps++;
1309 
1310  // What is the origin of the overlap?
1311  Bool_t mOK = 0, gOK = 0;
1312  Int_t mpdg = -999999, gpdg = -1;
1313  Int_t mstatus = -1, gstatus = -1;
1314  Int_t gLabel = -1, ggLabel = -1;
1315 
1316  GetMother (label[ilab],mcevent,mpdg,mstatus,mOK);
1317  fGMotherMom =
1318  GetGrandMother(label[ilab],mcevent,gpdg,gstatus,gOK, gLabel,ggLabel);
1319 
1320  //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
1321 
1322  if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
1323  ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
1324  gLabel >=0 )
1325  {
1326  Int_t labeltmp = gLabel;
1327  while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
1328  {
1329  mpdg=gpdg;
1330  fGMotherMom = GetGrandMother(labeltmp,mcevent,gpdg,gstatus,ok, gLabel,ggLabel);
1331  labeltmp=gLabel;
1332  }
1333  }
1334  overpdg [noverlaps-1] = mpdg;
1335  overlabel[noverlaps-1] = label[ilab];
1336  }
1337 
1338  return noverlaps ;
1339 }
1340 
1341 //________________________________________________________
1343 //________________________________________________________
1344 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1345 {
1346  if(! opt)
1347  return;
1348 
1349  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1350 
1351  printf("Debug level = %d\n",fDebug);
1352  printf("MC Generator = %s\n",fMCGeneratorString.Data());
1353  printf(" \n");
1354 }
1355 
1356 //__________________________________________________
1358 //__________________________________________________
1360 {
1361  printf("AliMCAnalysisUtils::PrintMCTag() - tag %d \n photon %d, conv %d, prompt %d, frag %d, isr %d, \n pi0 decay %d, eta decay %d, other decay %d pi0 %d, eta %d \n electron %d, muon %d,pion %d, proton %d, neutron %d, \n kaon %d, a-proton %d, a-neutron %d, unk %d, bad %d\n",
1362  tag,
1363  CheckTagBit(tag,kMCPhoton),
1365  CheckTagBit(tag,kMCPrompt),
1367  CheckTagBit(tag,kMCISR),
1368  CheckTagBit(tag,kMCPi0Decay),
1369  CheckTagBit(tag,kMCEtaDecay),
1371  CheckTagBit(tag,kMCPi0),
1372  CheckTagBit(tag,kMCEta),
1373  CheckTagBit(tag,kMCElectron),
1374  CheckTagBit(tag,kMCMuon),
1375  CheckTagBit(tag,kMCPion),
1376  CheckTagBit(tag,kMCProton),
1378  CheckTagBit(tag,kMCKaon),
1379  CheckTagBit(tag,kMCAntiProton),
1381  CheckTagBit(tag,kMCUnknown),
1383  );
1384 }
1385 
1386 //__________________________________________________
1388 //__________________________________________________
1390 {
1391  fMCGenerator = mcgen ;
1392  if (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
1393  else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
1394  else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
1395  else
1396  {
1397  fMCGeneratorString = "";
1399  }
1400 }
1401 
1402 //____________________________________________________
1404 //____________________________________________________
1406 {
1407  fMCGeneratorString = mcgen ;
1408 
1409  if (mcgen == "PYTHIA") fMCGenerator = kPythia;
1410  else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
1411  else if(mcgen == "HIJING") fMCGenerator = kHijing;
1412  else
1413  {
1415  fMCGeneratorString = "" ;
1416  }
1417 }
1418 
1419 
1420 
Int_t pdg
void SetMCGenerator(Int_t mcgen)
Set the generator type.
Int_t fCurrentEvent
Current Event number - GetJets()
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)
TLorentzVector fMotherMom
! particle momentum
void PrintMCTag(Int_t tag) const
Print the assigned origins to this particle.
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 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)
virtual ~AliMCAnalysisUtils()
Destructor.
void SetTagBit(Int_t &tag, UInt_t set) const
Int_t CheckOrigin(Int_t label, const AliMCEvent *mcevent)
TLorentzVector fGMotherMom
! particle momentum
Int_t fDebug
Debug level.
TLorentzVector fDaughMom2
! particle momentum
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 fMCGenerator
MC generator used to generate data in simulation.
TLorentzVector fDaughMom
! particle momentum
void CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex, const AliMCEvent *mcevent, Int_t &tag)
TList * fJetsList
List of jets - GetJets()
TList * GetJets(AliMCEvent *mcevent, AliGenEventHeader *mcheader, Int_t eventNumber)
TString fMCGeneratorString
MC generator used to generate data in simulation.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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)
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)