AliPhysics  1168478 (1168478)
 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 ( pPdg > 100 )
399  {
401 
402  AliDebug(2,Form("Generator decay photon from parent pdg %d",pPdg));
403  }
404  else if( mom->IsPhysicalPrimary() && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle
405  {
406  if(iParent < 8 && iParent > 5 )
407  {
408  //outgoing partons
409  if(pPdg == 22) SetTagBit(tag,kMCPrompt);
410  else SetTagBit(tag,kMCFragmentation);
411  }//Outgoing partons
412  else if( iParent <= 5 && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) )
413  {
414  SetTagBit(tag, kMCISR); //Initial state radiation
415  }
416  else if ( pPdg < 23 )
417  {
419 
420  AliDebug(2,Form("Generator fragmentation photon from parent pdg %d",pPdg));
421  }
422  else
423  {
424  AliDebug(2,Form("Generator physical primary (pythia/herwig) unknown photon from parent pdg %d",pPdg));
425 
426  SetTagBit(tag,kMCUnknown);
427  }
428  }//Physical primary
429  else
430  {
431  AliDebug(2,Form("Generator unknown photon from parent pdg %d",pPdg));
432 
433  SetTagBit(tag,kMCUnknown);
434  }
435 // //Old Herwig selection for ESDs, maybe should be considered.
436 // if(fMCGenerator == kHerwig)
437 // {
438 // if(pStatus < 197)
439 // {//Not decay
440 // while(1)
441 // {
442 // if(parent)
443 // {
444 // if(parent->GetMother()<=5) break;
445 // iParent = parent->GetMother();
446 // parent = mcevent->GetTrack(iParent);
447 // pStatus = parent->MCStatusCode();
448 // pPdg = TMath::Abs(parent->PdgCode());
449 // } else break;
450 // }//Look for the parton
451 //
452 // if(iParent < 8 && iParent > 5)
453 // {
454 // if(pPdg == 22) SetTagBit(tag,kMCPrompt);
455 // else SetTagBit(tag,kMCFragmentation);
456 // }
457 // else SetTagBit(tag,kMCISR);//Initial state radiation
458 // }//Not decay
459 // else SetTagBit(tag,kMCUnknown);
460 // }//HERWIG
461 
462  } // Mother Photon
463 
464  // Electron check. Where did that electron come from?
465  else if(mPdg == 11)
466  {
467  //electron
468  if(pPdg == 11 && parent)
469  {
470  Int_t iGrandma = parent->GetMother();
471  if(iGrandma >= 0)
472  {
473  AliVParticle * gma = mcevent->GetTrack(iGrandma);
474  Int_t gPdg = TMath::Abs(gma->PdgCode());
475 
476  if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
477  else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
478  }
479  }
480 
481  SetTagBit(tag,kMCElectron);
482 
483  AliDebug(1,"Checking ancestors of electrons");
484 
485  if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay
486  else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay
487  else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
488  else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
489  {
490  //c-hadron decay check
491  if(parent)
492  {
493  Int_t iGrandma = parent->GetMother();
494  if(iGrandma >= 0)
495  {
496  AliVParticle * gma = mcevent->GetTrack(iGrandma); //charm's mother
497  Int_t gPdg = TMath::Abs(gma->PdgCode());
498  if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
499  else SetTagBit(tag,kMCEFromC); //c-hadron decay
500  }
501  else SetTagBit(tag,kMCEFromC); //c-hadron decay
502  }//parent
503  } else
504  { //prompt or other decay
505 
506  TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
507  TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
508 
509  AliDebug(1,Form("Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)",foo->GetName(), pPdg,foo1->GetName(),mPdg));
510 
511  if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
512  else SetTagBit(tag,kMCOtherDecay);
513  }
514  }// electron check
515  // cluster was made by something else
516  else
517  {
518  AliDebug(1,Form("\t Setting kMCUnknown for cluster with pdg = %d, Parent pdg = %d",mPdg,pPdg));
519  SetTagBit(tag,kMCUnknown);
520  }
521 
522  return tag;
523 }
524 
525 //_________________________________________________________________________________________
528 //_________________________________________________________________________________________
530  Int_t mesonIndex, const AliMCEvent* mcevent,
531  Int_t & tag)
532 {
533  if(labels[0] < 0 || labels[0] > mcevent->GetNumberOfTracks() || nlabels <= 1)
534  {
535  AliDebug(2,Form("Exit : label[0] %d, n primaries %d, nlabels %d",labels[0],mcevent->GetNumberOfTracks(), nlabels));
536  return;
537  }
538 
539  AliVParticle * meson = mcevent->GetTrack(mesonIndex);
540  Int_t mesonPdg = meson->PdgCode();
541  if(mesonPdg != 111 && mesonPdg != 221)
542  {
543  AliWarning(Form("Wrong pi0/eta PDG : %d",mesonPdg));
544  return;
545  }
546 
547  AliDebug(2,Form("pdg %d, label %d, ndaughters %d", mesonPdg, mesonIndex, meson->GetNDaughters()));
548 
549  // Get the daughters
550  if(meson->GetNDaughters() != 2)
551  {
552 // if(meson->GetNDaughters()>2)
553 // {
554 // printf("xx More than 2 daughters <%d> for pdg %d with E %2.2f, pT %2.2f:\n",
555 // meson->GetNDaughters(),meson->PdgCode(),meson->E(),meson->Pt());
556 // for(Int_t idaugh = 0; idaugh < meson->GetNDaughters(); idaugh++)
557 // {
558 // if(meson->GetDaughterLabel(idaugh) < 0)
559 // {
560 // printf("\t no daughther with index %d:\n",idaugh);
561 // continue;
562 // }
563 // AliVParticle * daugh = mcevent->GetTrack(meson->GetDaughterLabel(idaugh));
564 // printf("\t daughther pdg %d, E %2.2f, pT %2.2f:\n",daugh->PdgCode(),daugh->E(),daugh->Pt());
565 // }
566 // }
567 
568  AliDebug(2,Form("Not overlapped. Number of daughters is %d, not 2",meson->GetNDaughters()));
569  return;
570  }
571 
572  Int_t iPhoton0 = meson->GetDaughterLabel(0);
573  Int_t iPhoton1 = meson->GetDaughterLabel(1);
574 
575  //if((iPhoton0 == -1) || (iPhoton1 == -1))
576  //{
577  // if(fDebug > 2)
578  // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
579  // return;
580  //}
581 
582  AliVParticle *photon0 = mcevent->GetTrack(iPhoton0);
583  AliVParticle *photon1 = mcevent->GetTrack(iPhoton1);
584 
585  // Check if both daughters are photons
586  if(photon0->PdgCode() != 22 && photon1->PdgCode()!=22)
587  {
588  AliWarning(Form("Not overlapped. PDG: daughter 1 = %d, of daughter 2 = %d",photon0->PdgCode(),photon1->PdgCode()));
589  return;
590  }
591 
592  AliDebug(2,Form("Daughter labels : photon0 = %d, photon1 = %d",iPhoton0,iPhoton1));
593 
594  // Check if both photons contribute to the cluster
595  Bool_t okPhoton0 = kFALSE;
596  Bool_t okPhoton1 = kFALSE;
597 
598  AliDebug(3,"Labels loop:");
599 
600  Bool_t conversion = kFALSE;
601 
602 // TLorentzVector d0, d1;
603 // photon0->Momentum(d0);
604 // photon1->Momentum(d1);
605 // Float_t openAngle = d0.Angle(d1.Vect());
606 
607  for(Int_t i = 0; i < nlabels; i++)
608  {
609  AliDebug(3, Form("\t label %d/%d: %d, ok? %d, %d", i, nlabels, labels[i], okPhoton0, okPhoton1));
610 
611  if ( labels[i] < 0 ) continue;
612 
613  // If we already found both, break the loop
614  if(okPhoton0 && okPhoton1) break;
615 
616  Int_t index = labels[i];
617  if (iPhoton0 == index)
618  {
619  okPhoton0 = kTRUE;
620  continue;
621  }
622  else if (iPhoton1 == index)
623  {
624  okPhoton1 = kTRUE;
625  continue;
626  }
627 
628  // Trace back the mother in case it was a conversion
629  if(index >= mcevent->GetNumberOfTracks())
630  {
631  AliWarning(Form("Particle index %d larger than size of list %d!!",index,mcevent->GetNumberOfTracks()));
632  continue;
633  }
634 
635  AliVParticle * daught = mcevent->GetTrack(index);
636  Int_t tmpindex = daught->GetMother();
637  AliDebug(3,Form("Conversion? : mother %d",tmpindex));
638 
639  while(tmpindex>=0)
640  {
641  // MC particle of interest is the mother
642  AliDebug(3,Form("\t parent index %d",tmpindex));
643  daught = mcevent->GetTrack(tmpindex);
644  //printf("tmpindex %d\n",tmpindex);
645  if (iPhoton0 == tmpindex)
646  {
647  conversion = kTRUE;
648  okPhoton0 = kTRUE;
649  break;
650  }
651  else if (iPhoton1 == tmpindex)
652  {
653  conversion = kTRUE;
654  okPhoton1 = kTRUE;
655  break;
656  }
657 
658  tmpindex = daught->GetMother();
659 
660  }//While to check if pi0/eta daughter was one of these contributors to the cluster
661 
662  //if(i == 0 && (!okPhoton0 && !okPhoton1)) AliDebug(1,"Something happens, first label should be from a photon decay!");
663 
664  }//loop on list of labels
665 
666  //If both photons contribute tag as the corresponding meson.
667  if(okPhoton0 && okPhoton1)
668  {
669  AliDebug(2,Form("%s OVERLAPPED DECAY",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()));
670 
671  if(!CheckTagBit(tag,kMCConversion) && conversion)
672  {
673  AliDebug(2,"Second decay photon produced a conversion");
674  SetTagBit(tag,kMCConversion) ;
675  }
676 
677  if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
678  else SetTagBit(tag,kMCEta);
679 
680 // printf("\t Meson E %2.2f, pT %2.2f, Overlapped open angle %2.2f deg, cells %2.1f, conversion %d\n",
681 // meson->E(), meson->Pt(),openAngle*TMath::RadToDeg(),openAngle/0.015,conversion);
682 //
683 // Float_t phi = photon0->Phi()*TMath::RadToDeg();
684 // if(phi < 0) phi+=360;
685 // printf("\t \t daughther1: label %d, E %2.2f, pT %2.2f, eta%2.2f, phi %2.2f, status %d, phys prim %d:\n",
686 // iPhoton0, photon0->E(),photon0->Pt(),
687 // photon0->Eta(),phi,
688 // photon0->MCStatusCode(),photon0->IsPhysicalPrimary());
689 //
690 // phi = photon1->Phi()*TMath::RadToDeg();
691 // if(phi < 0) phi+=360;
692 // printf("\t \t daughther2: label %d, E %2.2f, pT %2.2f, eta%2.2f, phi %2.2f, status %d, phys prim %d:\n",
693 // iPhoton1, photon1->E(),photon1->Pt(),
694 // photon1->Eta(),phi,
695 // photon1->MCStatusCode(),photon1->IsPhysicalPrimary());
696  }
697 }
698 
699 //________________________________________________________________________________________________________
701 //________________________________________________________________________________________________________
702 void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster, Int_t iMom, Int_t iParent,
703  const AliMCEvent* mcevent, Int_t & tag)
704 {
705  if(!arrayCluster || iMom < 0 || iParent < 0|| !mcevent) return;
706 
707  AliVParticle * parent = mcevent->GetTrack(iParent);
708 
709  //printf("*** Check label %d with parent %d\n",iMom, iParent);
710 
711  if(parent->GetNDaughters()!=2)
712  {
714  //printf("\t ndaugh = %d\n",parent->GetNDaughters());
715  return ;
716  }
717 
718  Int_t pairLabel = -1;
719  if ( iMom != parent->GetDaughterLabel(0) ) pairLabel = parent->GetDaughterLabel(0);
720  else if( iMom != parent->GetDaughterLabel(1) ) pairLabel = parent->GetDaughterLabel(1);
721 
722  if(pairLabel<0)
723  {
724  //printf("\t pair Label not found = %d\n",pairLabel);
726  return ;
727  }
728 
729  //printf("\t *** find pair %d\n",pairLabel);
730 
731  for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
732  {
733  AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
734  //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
735 
736  for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
737  {
738  Int_t label = cluster->GetLabels()[ilab];
739 
740  //printf("\t \t label %d\n",label);
741 
742  if ( label==pairLabel )
743  {
744  //printf("\t \t Pair found\n");
746  return ;
747  }
748  else if ( label== iParent || label== iMom )
749  {
750  //printf("\t \t skip\n");
751  continue;
752  }
753  else // check the ancestry
754  {
755  AliVParticle * mother = mcevent->GetTrack(label);
756 
757  if ( !mother )
758  {
759  AliInfo(Form("MC Mother not available for label %d",label));
760  continue;
761  }
762 
763  Int_t momPDG = TMath::Abs(mother->PdgCode());
764  if ( momPDG!=11 && momPDG!=22 ) continue;
765 
766  // Check if "mother" of entity is converted, if not, get the first non converted mother
767  Int_t iParentClus = mother->GetMother();
768  if(iParentClus < 0) continue;
769 
770  AliVParticle * parentClus = mcevent->GetTrack(iParentClus);
771  if(!parentClus) continue;
772 
773  Int_t parentClusPDG = TMath::Abs(parentClus->PdgCode());
774  Int_t parentClusStatus = parentClus->MCStatusCode();
775 
776  if ( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0 )
777  {
778  //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
779  continue;
780  }
781 
782  //printf("\t \t Conversion\n");
783 
784  // Check if the mother is photon or electron with status not stable
785  while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
786  {
787  //New Mother
788  label = iParentClus;
789  momPDG = parentClusPDG;
790 
791  iParentClus = parentClus->GetMother();
792  if ( iParentClus < 0 ) break;
793 
794  parentClus = mcevent->GetTrack(iParentClus);
795  if ( !parentClus ) break;
796 
797  parentClusPDG = TMath::Abs(parentClus->PdgCode());
798  parentClusStatus = parentClus->MCStatusCode() ;
799  }//while
800 
801  if ( (momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel) )
802  {
804  //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
805  return ;
806  }
807  else
808  {
809  //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
810  continue;
811  }
812 
813  }
814  }
815  } // cluster loop
816 
818 }
819 
820 //_____________________________________________________________________
822 //_____________________________________________________________________
823 TList * AliMCAnalysisUtils::GetJets(AliMCEvent* mcevent, AliGenEventHeader * mcheader, Int_t eventNumber)
824 {
825  if(fCurrentEvent==eventNumber) return fJetsList ;
826 
827  fCurrentEvent = eventNumber;
828 
829  if (fJetsList) fJetsList->Clear();
830  else fJetsList = new TList;
831 
832  Int_t nTriggerJets = 0;
833  Float_t tmpjet[]={0,0,0,0};
834 
835  //printf("Event %d %d\n",fCurrentEvent,iEvent);
836 
837  // Get outgoing partons
838  if(mcevent->GetNumberOfTracks() < 8) return fJetsList;
839 
840  AliVParticle * parton1 = mcevent->GetTrack(6);
841  AliVParticle * parton2 = mcevent->GetTrack(7);
842 
843  AliDebug(2,Form("Parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
844  parton1->GetName(),parton1->Pt(),parton1->E(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()));
845  AliDebug(2,Form("Parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
846  parton2->GetName(),parton2->Pt(),parton2->E(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()));
847 
848  // //Trace the jet from the mother parton
849  // Float_t pt = 0;
850  // Float_t pt1 = 0;
851  // Float_t pt2 = 0;
852  // Float_t e = 0;
853  // Float_t e1 = 0;
854  // Float_t e2 = 0;
855  // TParticle * tmptmp = new TParticle;
856  // for(Int_t i = 0; i< stack->GetNprimary(); i++){
857  // tmptmp = stack->Particle(i);
858 
859  // if(tmptmp->GetStatusCode() == 1){
860  // pt = tmptmp->Pt();
861  // e = tmptmp->Energy();
862  // Int_t imom = tmptmp->GetMother();
863  // Int_t imom1 = 0;
864  // //printf("1st imom %d\n",imom);
865  // while(imom > 5){
866  // imom1=imom;
867  // tmptmp = stack->Particle(imom);
868  // imom = tmptmp->GetMother();
869  // //printf("imom %d \n",imom);
870  // }
871  // //printf("Last imom %d %d\n",imom1, imom);
872  // if(imom1 == 6) {
873  // pt1+=pt;
874  // e1+=e;
875  // }
876  // else if (imom1 == 7){
877  // pt2+=pt;
878  // e2+=e; }
879  // }// status 1
880 
881  // }// for
882 
883  // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
884 
885  //Get the jet, different way for different generator
886  //PYTHIA
887  if(fMCGenerator == kPythia)
888  {
889  TParticle * jet = 0x0;
890  AliGenPythiaEventHeader* pygeh = (AliGenPythiaEventHeader*) mcheader;
891  nTriggerJets = pygeh->NTriggerJets();
892  AliDebug(2,Form("PythiaEventHeader: Njets: %d",nTriggerJets));
893 
894  for(Int_t i = 0; i< nTriggerJets; i++)
895  {
896  pygeh->TriggerJet(i, tmpjet);
897 
898  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
899 
900  // Assign an outgoing parton as mother
901  Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
902  Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
903 
904  if(phidiff1 > phidiff2) jet->SetFirstMother(7);
905  else jet->SetFirstMother(6);
906 
907  //jet->Print();
908  AliDebug(1,Form("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
909  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()));
910  fJetsList->Add(jet);
911  }
912  }//Pythia triggered jets
913  //HERWIG
914  else if (fMCGenerator == kHerwig)
915  {
916  Int_t pdg = -1;
917 
918  // Check parton 1
919  AliVParticle * tmp = parton1;
920  if(parton1->PdgCode()!=22)
921  {
922  while(pdg != 94)
923  {
924  if(tmp->GetFirstDaughter()==-1) return fJetsList;
925 
926  tmp = mcevent->GetTrack(tmp->GetFirstDaughter());
927  pdg = tmp->PdgCode();
928  }//while
929 
930  // Add found jet to list
931  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);
932 
933  jet1->SetFirstMother(6);
934 
935  fJetsList->Add(jet1);
936 
937  //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
938  //tmp = stack->Particle(tmp->GetFirstDaughter());
939  //tmp->Print();
940  //jet1->Print();
941  AliDebug(1,Form("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
942  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()));
943  }//not photon
944 
945  // Check parton 2
946  pdg = -1;
947  tmp = parton2;
948  if(parton2->PdgCode()!=22)
949  {
950  while(pdg != 94)
951  {
952  if(tmp->GetFirstDaughter()==-1) return fJetsList;
953 
954  tmp = mcevent->GetTrack(tmp->GetFirstDaughter());
955  pdg = tmp->PdgCode();
956  }//while
957 
958  // Add found jet to list
959  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);
960 
961  jet2->SetFirstMother(7);
962 
963  fJetsList->Add(jet2);
964 
965  //jet2->Print();
966  AliDebug(2,Form("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
967  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()));
968  //Int_t first = tmp->GetFirstDaughter();
969  //Int_t last = tmp->GetLastDaughter();
970  //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
971  // for(Int_t d = first ; d < last+1; d++){
972  // tmp = stack->Particle(d);
973  // if(i == tmp->GetMother())
974  // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
975  // d,tmp->GetMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
976  // }
977  //tmp->Print();
978  }//not photon
979  }//Herwig generated jets
980 
981  return fJetsList;
982 }
983 
984 //______________________________________________________________________________
986 //______________________________________________________________________________
987 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
988  const AliMCEvent* mcevent,
989  Int_t & pdg, Int_t & status,
990  Bool_t & ok, Int_t & daughlabel, TVector3 & prodVertex)
991 {
992  fDaughMom.SetPxPyPzE(0,0,0,0);
993 
994  if(!mcevent)
995  {
996  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
997 
998  ok = kFALSE;
999  return fDaughMom;
1000  }
1001 
1002  Int_t nprimaries = mcevent->GetNumberOfTracks();
1003  if(label < 0 || label >= nprimaries)
1004  {
1005  ok = kFALSE;
1006  return fDaughMom;
1007  }
1008 
1009  AliVParticle * momP = mcevent->GetTrack(label);
1010  daughlabel = momP->GetDaughterLabel(idaugh);
1011 
1012  if(daughlabel < 0 || daughlabel >= nprimaries)
1013  {
1014  ok = kFALSE;
1015  return fDaughMom;
1016  }
1017 
1018  AliVParticle * daughP = mcevent->GetTrack(daughlabel);
1019  fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1020  pdg = daughP->PdgCode();
1021  status = daughP->MCStatusCode();
1022  prodVertex.SetXYZ(daughP->Xv(),daughP->Yv(),daughP->Zv());
1023 
1024  ok = kTRUE;
1025 
1026  return fDaughMom;
1027 }
1028 
1029 //______________________________________________________________________________________________________
1031 //______________________________________________________________________________________________________
1032 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent, Bool_t & ok)
1033 {
1034  Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1035 
1036  return GetMother(label,mcevent,pdg,status, ok,momlabel);
1037 }
1038 
1039 //_________________________________________________________________________________________
1040 // \return the kinematics of the particle that generated the signal.
1041 //_________________________________________________________________________________________
1042 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent,
1043  Int_t & pdg, Int_t & status, Bool_t & ok)
1044 {
1045  Int_t momlabel = -1;
1046 
1047  return GetMother(label,mcevent,pdg,status, ok,momlabel);
1048 }
1049 
1050 //______________________________________________________________________________________________________
1052 //______________________________________________________________________________________________________
1053 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliMCEvent* mcevent,
1054  Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1055 {
1056  fMotherMom.SetPxPyPzE(0,0,0,0);
1057 
1058  if(!mcevent)
1059  {
1060  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1061 
1062  ok = kFALSE;
1063  return fMotherMom;
1064  }
1065 
1066  Int_t nprimaries = mcevent->GetNumberOfTracks();
1067  if(label < 0 || label >= nprimaries)
1068  {
1069  ok = kFALSE;
1070  return fMotherMom;
1071  }
1072 
1073  AliVParticle * momP = mcevent->GetTrack(label);
1074  fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1075  pdg = momP->PdgCode();
1076  status = momP->MCStatusCode();
1077  momlabel = momP->GetMother();
1078 
1079  ok = kTRUE;
1080 
1081  return fMotherMom;
1082 }
1083 
1084 //___________________________________________________________________________________
1086 //___________________________________________________________________________________
1088  const AliMCEvent* mcevent,
1089  Bool_t & ok, Int_t & momlabel)
1090 {
1091  fGMotherMom.SetPxPyPzE(0,0,0,0);
1092 
1093  if ( !mcevent )
1094  {
1095  AliWarning("MCEvent is not available, check analysis settings in configuration file!!");
1096 
1097  ok = kFALSE;
1098  return fGMotherMom;
1099  }
1100 
1101  Int_t nprimaries = mcevent->GetNumberOfTracks();
1102  if ( label < 0 || label >= nprimaries )
1103  {
1104  ok = kFALSE;
1105  return fGMotherMom;
1106  }
1107 
1108 
1109  AliVParticle * momP = mcevent->GetTrack(label);
1110 
1111  if(momP->PdgCode()==pdg)
1112  {
1113  AliDebug(2,"PDG of mother is already the one requested!");
1114  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1115 
1116  ok=kTRUE;
1117  return fGMotherMom;
1118  }
1119 
1120  Int_t grandmomLabel = momP->GetMother();
1121  Int_t grandmomPDG = -1;
1122  AliVParticle * grandmomP = 0x0;
1123 
1124  while (grandmomLabel >=0 )
1125  {
1126  grandmomP = mcevent->GetTrack(grandmomLabel);
1127  grandmomPDG = grandmomP->PdgCode();
1128  if(grandmomPDG==pdg)
1129  {
1130  //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1131  momlabel = grandmomLabel;
1132  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1133  break;
1134  }
1135 
1136  grandmomLabel = grandmomP->GetMother();
1137  }
1138 
1139  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found!",pdg));
1140 
1141  ok = kTRUE;
1142 
1143  return fGMotherMom;
1144 }
1145 
1146 //______________________________________________________________________________________________
1148 //______________________________________________________________________________________________
1149 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliMCEvent* mcevent,
1150  Int_t & pdg, Int_t & status, Bool_t & ok,
1151  Int_t & grandMomLabel, Int_t & greatMomLabel)
1152 {
1153  fGMotherMom.SetPxPyPzE(0,0,0,0);
1154 
1155  if ( !mcevent )
1156  {
1157  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1158 
1159  ok = kFALSE;
1160  return fGMotherMom;
1161  }
1162 
1163  Int_t nprimaries = mcevent->GetNumberOfTracks();
1164  if ( label < 0 || label >= nprimaries )
1165  {
1166  ok = kFALSE;
1167  return fGMotherMom;
1168  }
1169 
1170  AliVParticle * momP = mcevent->GetTrack(label);
1171 
1172  grandMomLabel = momP->GetMother();
1173 
1174  AliVParticle * grandmomP = 0x0;
1175 
1176  if(grandMomLabel >=0 )
1177  {
1178  grandmomP = mcevent->GetTrack(grandMomLabel);
1179  pdg = grandmomP->PdgCode();
1180  status = grandmomP->MCStatusCode();
1181 
1182  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1183  greatMomLabel = grandmomP->GetMother();
1184 
1185  }
1186 
1187  ok = kTRUE;
1188 
1189  return fGMotherMom;
1190 }
1191 
1192 //_______________________________________________________________________________________________________________
1194 //_______________________________________________________________________________________________________________
1196  Float_t & asym, Float_t & angle, Bool_t & ok)
1197 {
1198  if(!mcevent)
1199  {
1200  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1201 
1202  ok = kFALSE;
1203  return;
1204  }
1205 
1206  Int_t nprimaries = mcevent->GetNumberOfTracks();
1207  if ( label < 0 || label >= nprimaries )
1208  {
1209  ok = kFALSE;
1210  return ;
1211  }
1212 
1213  AliVParticle * momP = mcevent->GetTrack(label);
1214 
1215  Int_t grandmomLabel = momP->GetMother();
1216  Int_t grandmomPDG = -1;
1217  AliVParticle * grandmomP = 0x0;
1218 
1219  while (grandmomLabel >=0 )
1220  {
1221  grandmomP = mcevent->GetTrack(grandmomLabel);
1222  grandmomPDG = grandmomP->PdgCode();
1223 
1224  if(grandmomPDG==pdg) break;
1225 
1226  grandmomLabel = grandmomP->GetMother();
1227  }
1228 
1229  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1230  {
1231  AliVParticle * d1 = mcevent->GetTrack(grandmomP->GetDaughterLabel(0));
1232  AliVParticle * d2 = mcevent->GetTrack(grandmomP->GetDaughterLabel(1));
1233 
1234  if(d1->PdgCode() == 22 && d1->PdgCode() == 22)
1235  {
1236  asym = (d1->E()-d2->E())/grandmomP->E();
1237  fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1238  fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1239  angle = fDaughMom.Angle(fDaughMom2.Vect());
1240  }
1241  }
1242  else
1243  {
1244  ok = kFALSE;
1245  AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
1246  return;
1247  }
1248 
1249  ok = kTRUE;
1250 }
1251 
1252 //_________________________________________________________________________________________________
1254 //_________________________________________________________________________________________________
1255 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliMCEvent* mcevent, Bool_t & ok)
1256 {
1257  if ( !mcevent )
1258  {
1259  AliWarning("MCEvent is not available, check analysis settings in configuration file, STOP!!");
1260 
1261  ok=kFALSE;
1262  return -1;
1263  }
1264 
1265  Int_t nprimaries = mcevent->GetNumberOfTracks();
1266  if ( label < 0 || label >= nprimaries )
1267  {
1268  ok = kFALSE;
1269  return -1;
1270  }
1271 
1272  AliVParticle * momP = mcevent->GetTrack(label);
1273 
1274  ok = kTRUE;
1275 
1276  return momP->GetNDaughters();
1277 }
1278 
1279 //_________________________________________________________________________________
1284 //_________________________________________________________________________________
1286  Int_t mctag, Int_t mesonLabel,
1287  AliMCEvent* mcevent,
1288  Int_t *overpdg, Int_t *overlabel)
1289 {
1290  Int_t ancPDG = 0, ancStatus = -1;
1291  TVector3 prodVertex;
1292  Int_t ancLabel = 0;
1293  Int_t noverlaps = 0;
1294  Bool_t ok = kFALSE;
1295 
1296  for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1297  {
1298  ancLabel = CheckCommonAncestor(label[0],label[ilab],mcevent,ancPDG,ancStatus,fMotherMom,prodVertex);
1299 
1300  //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1301  // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1302 
1303  Bool_t overlap = kFALSE;
1304 
1305  if ( ancLabel < 0 )
1306  {
1307  overlap = kTRUE;
1308  //printf("\t \t \t No Label = %d\n",ancLabel);
1309  }
1310  else if ( ( ancPDG==111 || ancPDG==221 ) &&
1311  ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta) ) &&
1312  ( (mesonLabel != ancLabel) && mesonLabel >=0 ) ) // in case the label is not provided check it is larger than 0
1313  {
1314  //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1315  overlap = kTRUE;
1316  }
1317  else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1318  {
1319  //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1320  overlap = kTRUE ;
1321  }
1322 
1323  if( !overlap ) continue ;
1324 
1325  // We have at least one overlap
1326 
1327  //printf("Overlap!!!!!!!!!!!!!!\n");
1328 
1329  noverlaps++;
1330 
1331  // What is the origin of the overlap?
1332  Bool_t mOK = 0, gOK = 0;
1333  Int_t mpdg = -999999, gpdg = -1;
1334  Int_t mstatus = -1, gstatus = -1;
1335  Int_t gLabel = -1, ggLabel = -1;
1336 
1337  GetMother (label[ilab],mcevent,mpdg,mstatus,mOK);
1338  fGMotherMom =
1339  GetGrandMother(label[ilab],mcevent,gpdg,gstatus,gOK, gLabel,ggLabel);
1340 
1341  //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
1342 
1343  if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
1344  ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
1345  gLabel >=0 )
1346  {
1347  Int_t labeltmp = gLabel;
1348  while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
1349  {
1350  mpdg=gpdg;
1351  fGMotherMom = GetGrandMother(labeltmp,mcevent,gpdg,gstatus,ok, gLabel,ggLabel);
1352  labeltmp=gLabel;
1353  }
1354  }
1355  overpdg [noverlaps-1] = mpdg;
1356  overlabel[noverlaps-1] = label[ilab];
1357  }
1358 
1359  return noverlaps ;
1360 }
1361 
1362 //________________________________________________________
1364 //________________________________________________________
1365 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1366 {
1367  if(! opt)
1368  return;
1369 
1370  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1371 
1372  printf("Debug level = %d\n",fDebug);
1373  printf("MC Generator = %s\n",fMCGeneratorString.Data());
1374  printf(" \n");
1375 }
1376 
1377 //________________________________________________________
1385 //________________________________________________________
1386 void AliMCAnalysisUtils::PrintAncestry(AliMCEvent* mcevent, Int_t label, Int_t nGenerMax) const
1387 {
1388  AliVParticle * primary = 0;
1389  Int_t index = label;
1390  Int_t gener = 0;
1391 
1392  AliInfo("*********** Start");
1393  while ( index > 0 )
1394  {
1395  primary = mcevent->GetTrack(index);
1396 
1397  if(!primary)
1398  {
1399  AliWarning("primary pointer not available!!");
1400  return;
1401  }
1402 
1403  Float_t eta = 0;
1404  // Protection against floating point exception
1405  if ( primary->E() == TMath::Abs(primary->Pz()) ||
1406  (primary->E() - primary->Pz()) < 1e-3 ||
1407  (primary->E() + primary->Pz()) < 0 )
1408  eta = -999;
1409  else
1410  eta = primary->Eta();
1411 
1412  Int_t pdg = primary->PdgCode();
1413 
1414  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"
1415  " mother %d, n daughters %d, d1 %d, d2 %d\n",
1416  gener,index,TDatabasePDG::Instance()->GetParticle(pdg)->GetName(),pdg,primary->MCStatusCode(),primary->IsPhysicalPrimary(),
1417  primary->E(),primary->Pt(),primary->Px(),primary->Py(),primary->Pz(),
1418  eta,primary->Phi()*TMath::RadToDeg(),
1419  primary->GetMother(),primary->GetNDaughters(),primary->GetDaughterLabel(0), primary->GetDaughterLabel(1));
1420 
1421  gener++;
1422  index = primary->GetMother();
1423  if ( nGenerMax < gener ) index = -1; // stop digging ancestry
1424  } // while
1425 
1426  AliInfo("*********** End");
1427 }
1428 
1429 
1430 //__________________________________________________
1432 //__________________________________________________
1434 {
1435  AliInfo
1436  (
1437  Form
1438  ("Tag %d: photon %d, conv %d, prompt %d, frag %d, isr %d,\n"
1439  " pi0 decay %d, eta decay %d, other decay %d, lost decay %d, in calo decay %d, pi0 %d, eta %d,\n"
1440  " electron %d, muon %d,pion %d, proton %d, neutron %d,\n"
1441  " kaon %d, a-proton %d, a-neutron %d, unk %d, bad %d",
1442  tag,
1443  CheckTagBit(tag,kMCPhoton),
1445  CheckTagBit(tag,kMCPrompt),
1447  CheckTagBit(tag,kMCISR),
1448  CheckTagBit(tag,kMCPi0Decay),
1449  CheckTagBit(tag,kMCEtaDecay),
1453  CheckTagBit(tag,kMCPi0),
1454  CheckTagBit(tag,kMCEta),
1455  CheckTagBit(tag,kMCElectron),
1456  CheckTagBit(tag,kMCMuon),
1457  CheckTagBit(tag,kMCPion),
1458  CheckTagBit(tag,kMCProton),
1460  CheckTagBit(tag,kMCKaon),
1461  CheckTagBit(tag,kMCAntiProton),
1463  CheckTagBit(tag,kMCUnknown),
1465  )
1466  );
1467 }
1468 
1469 //__________________________________________________
1471 //__________________________________________________
1473 {
1474  fMCGenerator = mcgen ;
1475  if (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
1476  else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
1477  else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
1478  else
1479  {
1480  fMCGeneratorString = "";
1482  }
1483 }
1484 
1485 //____________________________________________________
1487 //____________________________________________________
1489 {
1490  fMCGeneratorString = mcgen ;
1491 
1492  if (mcgen == "PYTHIA") fMCGenerator = kPythia;
1493  else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
1494  else if(mcgen == "HIJING") fMCGenerator = kHijing;
1495  else
1496  {
1498  fMCGeneratorString = "" ;
1499  }
1500 }
1501 
1502 
1503 
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.
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)