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