1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // Class for analysis utils for MC data
18 // stored in stack or event header.
20 // - method to check the origin of a given track/cluster
21 // - method to obtain the generated jets
23 //*-- Author: Gustavo Conesa (LNF-INFN)
24 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
30 #include "TParticle.h"
31 #include "TDatabasePDG.h"
34 //---- ANALYSIS system ----
35 #include "AliMCAnalysisUtils.h"
36 #include "AliCaloTrackReader.h"
38 #include "AliGenPythiaEventHeader.h"
39 #include "AliAODMCParticle.h"
41 ClassImp(AliMCAnalysisUtils)
43 //________________________________________
44 AliMCAnalysisUtils::AliMCAnalysisUtils() :
49 fMCGenerator("PYTHIA")
54 //_______________________________________
55 AliMCAnalysisUtils::~AliMCAnalysisUtils()
57 // Remove all pointers.
65 //_____________________________________________________________________________________________
66 Int_t AliMCAnalysisUtils::CheckCommonAncestor(Int_t index1, Int_t index2,
67 const AliCaloTrackReader* reader,
68 Int_t & ancPDG, Int_t & ancStatus,
69 TLorentzVector & momentum, TVector3 & prodVertex)
71 //Check the first common ancestor of 2 clusters, given the most likely labels of the primaries generating such clusters.
79 if(label1[0]==label2[0])
81 //printf("AliMCAnalysisUtils::CheckCommonAncestor() - Already the same label: %d\n",label1[0]);
87 if(reader->ReadAODMCParticles())
89 TClonesArray * mcparticles = reader->GetAODMCParticles();
91 Int_t label=label1[0];
92 if(label < 0) return -1;
94 while(label > -1 && counter1 < 99)
97 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
100 label = mom->GetMother() ;
101 label1[counter1]=label;
103 //printf("\t counter %d, label %d\n", counter1,label);
106 //printf("Org label2=%d,\n",label2[0]);
108 if(label < 0) return -1;
110 while(label > -1 && counter2 < 99)
113 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
116 label = mom->GetMother() ;
117 label2[counter2]=label;
119 //printf("\t counter %d, label %d\n", counter2,label);
123 { //Kine stack from ESDs
124 AliStack * stack = reader->GetStack();
126 Int_t label=label1[0];
127 while(label > -1 && counter1 < 99)
130 TParticle * mom = stack->Particle(label);
133 label = mom->GetFirstMother() ;
134 label1[counter1]=label;
136 //printf("\t counter %d, label %d\n", counter1,label);
139 //printf("Org label2=%d,\n",label2[0]);
142 while(label > -1 && counter2 < 99)
145 TParticle * mom = stack->Particle(label);
148 label = mom->GetFirstMother() ;
149 label2[counter2]=label;
151 //printf("\t counter %d, label %d\n", counter2,label);
153 }// Kine stack from ESDs
154 }//First labels not the same
156 if((counter1==99 || counter2==99) && fDebug >=0)
157 printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
158 //printf("CheckAncestor:\n");
160 Int_t commonparents = 0;
162 //printf("counters %d %d \n",counter1, counter2);
163 for (Int_t c1 = 0; c1 < counter1; c1++)
165 for (Int_t c2 = 0; c2 < counter2; c2++)
167 if(label1[c1]==label2[c2] && label1[c1]>-1)
169 ancLabel = label1[c1];
172 if(reader->ReadAODMCParticles())
174 AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]);
178 ancPDG = mom->GetPdgCode();
179 ancStatus = mom->GetStatus();
180 momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
181 prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
186 TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
190 ancPDG = mom->GetPdgCode();
191 ancStatus = mom->GetStatusCode();
192 mom->Momentum(momentum);
193 prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz());
197 //First ancestor found, end the loops
201 }//second cluster loop
202 }//first cluster loop
208 momentum.SetXYZT(0,0,0,0);
209 prodVertex.SetXYZ(-10,-10,-10);
215 //_______________________________________________________________________
216 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, Int_t nlabels,
217 const AliCaloTrackReader* reader)
219 //Play with the montecarlo particles if available
223 printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
227 //Select where the information is, ESD-galice stack or AOD mcparticles branch
228 if(reader->ReadStack()){
229 tag = CheckOriginInStack(label, nlabels, reader->GetStack());
231 else if(reader->ReadAODMCParticles()){
232 tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles());
238 //__________________________________________________________________________________
239 Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, const AliCaloTrackReader* reader)
241 //Play with the montecarlo particles if available
245 printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
249 Int_t labels[]={label};
251 //Select where the information is, ESD-galice stack or AOD mcparticles branch
252 if(reader->ReadStack()){
253 tag = CheckOriginInStack(labels, 1,reader->GetStack());
255 else if(reader->ReadAODMCParticles()){
256 tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles());
262 //_______________________________________________________________________________________________
263 Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, Int_t nlabels, AliStack* stack)
265 // Play with the MC stack if available. Tag particles depending on their origin.
266 // Do same things as in CheckOriginInAOD but different input.
268 //generally speaking, label is the MC label of a reconstructed
269 //entity (track, cluster, etc) for which we want to know something
270 //about its heritage, but one can also use it directly with stack
271 //particles not connected to reconstructed entities
276 printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
281 Int_t label=labels[0];//Most significant particle contributing to the cluster
283 if(label >= 0 && label < stack->GetNtrack())
285 //MC particle of interest is the "mom" of the entity
286 TParticle * mom = stack->Particle(label);
288 Int_t mPdgSign = mom->GetPdgCode();
289 Int_t mPdg = TMath::Abs(mPdgSign);
290 Int_t mStatus = mom->GetStatusCode() ;
291 Int_t iParent = mom->GetFirstMother() ;
293 if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
295 //GrandParent of the entity
296 TParticle * parent = NULL;
301 parent = stack->Particle(iParent);
305 pPdg = TMath::Abs(parent->GetPdgCode());
306 pStatus = parent->GetStatusCode();
309 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
313 printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
314 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
315 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
318 //Check if "mother" of entity is converted, if not, get the first non converted mother
319 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0)
321 SetTagBit(tag,kMCConversion);
323 //Check if the mother is photon or electron with status not stable
324 while ((pPdg == 22 || pPdg == 11) && mStatus != 1)
327 iMom = mom->GetFirstMother();
328 mom = stack->Particle(iMom);
329 mPdgSign = mom->GetPdgCode();
330 mPdg = TMath::Abs(mPdgSign);
331 mStatus = mom->GetStatusCode() ;
332 iParent = mom->GetFirstMother() ;
334 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
339 parent = stack->Particle(iParent);
342 pPdg = TMath::Abs(parent->GetPdgCode());
343 pStatus = parent->GetStatusCode();
346 else {// in case of gun/box simulations
355 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
356 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
357 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
360 }//mother and parent are electron or photon and have status 0
361 else if((mPdg == 22 || mPdg == 11) && mStatus == 0)
363 //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
364 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
365 pPdg == 2212 || pPdg == 130 || pPdg == 13 )
367 SetTagBit(tag,kMCConversion);
368 iMom = mom->GetFirstMother();
369 mom = stack->Particle(iMom);
370 mPdgSign = mom->GetPdgCode();
371 mPdg = TMath::Abs(mPdgSign);
375 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
376 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
380 //Comment for the next lines, we do not check the parent of the hadron for the moment.
381 //iParent = mom->GetFirstMother() ;
382 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
386 // parent = stack->Particle(iParent);
387 // pPdg = TMath::Abs(parent->GetPdgCode());
390 // conversion into electrons/photons checked
392 //first check for typical charged particles
393 if (mPdg == 13) SetTagBit(tag,kMCMuon);
394 else if(mPdg == 211) SetTagBit(tag,kMCPion);
395 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
396 else if(mPdgSign == 2212) SetTagBit(tag,kMCProton);
397 else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
398 else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron);
399 else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
401 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
405 SetTagBit(tag,kMCPi0Decay);
408 printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
410 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
414 SetTagBit(tag,kMCEtaDecay);
417 printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
419 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
424 SetTagBit(tag,kMCPhoton);
428 SetTagBit(tag,kMCPi0Decay);
430 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
432 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
434 else if (pPdg == 221)
436 SetTagBit(tag, kMCEtaDecay);
438 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
440 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
442 else if(mStatus == 1)
443 { //undecayed particle
444 if(fMCGenerator == "PYTHIA")
446 if(iParent < 8 && iParent > 5)
448 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
449 else SetTagBit(tag,kMCFragmentation);
451 else if(iParent <= 5)
453 SetTagBit(tag, kMCISR); //Initial state radiation
455 else SetTagBit(tag,kMCUnknown);
458 else if(fMCGenerator == "HERWIG")
466 if(parent->GetFirstMother()<=5) break;
467 iParent = parent->GetFirstMother();
468 parent=stack->Particle(iParent);
469 pStatus= parent->GetStatusCode();
470 pPdg = TMath::Abs(parent->GetPdgCode());
472 }//Look for the parton
474 if(iParent < 8 && iParent > 5)
476 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
477 else SetTagBit(tag,kMCFragmentation);
479 else SetTagBit(tag,kMCISR);//Initial state radiation
481 else SetTagBit(tag,kMCUnknown);
484 else SetTagBit(tag,kMCOtherDecay);
488 //Electron check. Where did that electron come from?
489 else if(mPdg == 11){ //electron
490 if(pPdg == 11 && parent)
492 Int_t iGrandma = parent->GetFirstMother();
495 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
496 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
498 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
499 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
503 SetTagBit(tag,kMCElectron);
505 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
507 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
508 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
509 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
510 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
513 Int_t iGrandma = parent->GetFirstMother();
516 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
517 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
518 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
519 else SetTagBit(tag,kMCEFromC); //c-->e
521 else SetTagBit(tag,kMCEFromC); //c-->e
526 //if it is not from any of the above, where is it from?
527 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
529 else SetTagBit(tag,kMCOtherDecay);
531 if(fDebug > 0 && parent) printf("AliMCAnalysisUtils::CheckOriginInStack() - Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg);
534 //Cluster was made by something else
538 printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",
539 mom->GetName(),mPdg,pPdg);
541 SetTagBit(tag,kMCUnknown);
546 if(label < 0 && (fDebug >= 0))
547 printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
549 if(label >= stack->GetNtrack() && (fDebug >= 0))
550 printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
552 SetTagBit(tag,kMCUnknown);
560 //____________________________________________________________________________
561 Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, Int_t nlabels,
562 const TClonesArray *mcparticles)
564 // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
565 // Do same things as in CheckOriginInStack but different input.
569 printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
574 Int_t label=labels[0];//Most significant particle contributing to the cluster
576 Int_t nprimaries = mcparticles->GetEntriesFast();
577 if(label >= 0 && label < nprimaries)
580 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
582 Int_t mPdgSign = mom->GetPdgCode();
583 Int_t mPdg = TMath::Abs(mPdgSign);
584 Int_t iParent = mom->GetMother() ;
586 if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
589 AliAODMCParticle * parent = NULL ;
593 parent = (AliAODMCParticle *) mcparticles->At(iParent);
594 pPdg = TMath::Abs(parent->GetPdgCode());
596 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
600 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
602 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
603 iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
606 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
607 iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
610 //Check if mother is converted, if not, get the first non converted mother
611 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary())
613 SetTagBit(tag,kMCConversion);
615 //Check if the mother is photon or electron with status not stable
616 while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
619 iMom = mom->GetMother();
620 mom = (AliAODMCParticle *) mcparticles->At(iMom);
621 mPdgSign = mom->GetPdgCode();
622 mPdg = TMath::Abs(mPdgSign);
623 iParent = mom->GetMother() ;
624 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
627 if(iParent >= 0 && parent)
629 parent = (AliAODMCParticle *) mcparticles->At(iParent);
630 pPdg = TMath::Abs(parent->GetPdgCode());
632 // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
633 // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
639 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
641 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
642 iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
645 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
646 iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
649 }//mother and parent are electron or photon and have status 0 and parent is photon or electron
650 else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary())
652 //Still a conversion but only one electron/photon generated. Just from hadrons
653 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
654 pPdg == 2212 || pPdg == 130 || pPdg == 13 )
656 SetTagBit(tag,kMCConversion);
657 iMom = mom->GetMother();
658 mom = (AliAODMCParticle *) mcparticles->At(iMom);
659 mPdgSign = mom->GetPdgCode();
660 mPdg = TMath::Abs(mPdgSign);
664 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
665 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
669 //Comment for next lines, we do not check the parent of the hadron for the moment.
670 //iParent = mom->GetMother() ;
671 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
675 // parent = (AliAODMCParticle *) mcparticles->At(iParent);
676 // pPdg = TMath::Abs(parent->GetPdgCode());
680 //printf("Final mother mPDG %d\n",mPdg);
682 // conversion into electrons/photons checked
684 //first check for typical charged particles
685 if (mPdg == 13) SetTagBit(tag,kMCMuon);
686 else if(mPdg == 211) SetTagBit(tag,kMCPion);
687 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
688 else if(mPdgSign == 2212) SetTagBit(tag,kMCProton);
689 else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron);
690 else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
691 else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
693 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
696 SetTagBit(tag,kMCPi0Decay);
698 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
700 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
704 SetTagBit(tag,kMCEtaDecay);
706 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
708 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
713 SetTagBit(tag,kMCPhoton);
717 SetTagBit(tag,kMCPi0Decay);
719 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
721 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
723 else if (pPdg == 221)
725 SetTagBit(tag, kMCEtaDecay);
727 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
729 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
731 else if(mom->IsPhysicalPrimary() && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) //undecayed particle
733 if(iParent < 8 && iParent > 5 )
735 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
736 else SetTagBit(tag,kMCFragmentation);
738 else if(iParent <= 5 && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG"))
740 SetTagBit(tag, kMCISR); //Initial state radiation
742 else SetTagBit(tag,kMCUnknown);
744 else SetTagBit(tag,kMCOtherDecay);
748 //Electron check. Where did that electron come from?
751 if(pPdg == 11 && parent)
753 Int_t iGrandma = parent->GetMother();
756 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
757 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
759 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
760 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
764 SetTagBit(tag,kMCElectron);
766 if (fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
767 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
768 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
769 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
770 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
771 { //c-hadron decay check
774 Int_t iGrandma = parent->GetMother();
777 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
778 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
779 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
780 else SetTagBit(tag,kMCEFromC); //c-hadron decay
782 else SetTagBit(tag,kMCEFromC); //c-hadron decay
785 { //prompt or other decay
786 TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
787 TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
789 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",
790 foo->GetName(), pPdg,foo1->GetName(),mPdg);
792 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
793 else SetTagBit(tag,kMCOtherDecay);
796 //cluster was made by something else
799 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",
801 SetTagBit(tag,kMCUnknown);
806 if(label < 0 && (fDebug >= 0) )
807 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label);
809 if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) )
810 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
812 SetTagBit(tag,kMCUnknown);
820 //_________________________________________________________________________________________
821 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels,
822 Int_t mesonIndex, AliStack *stack,
825 //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
827 if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
828 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
829 labels[0],stack->GetNtrack(), nlabels);
833 TParticle * meson = stack->Particle(mesonIndex);
834 Int_t mesonPdg = meson->GetPdgCode();
835 if(mesonPdg!=111 && mesonPdg!=221)
837 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
841 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s, label %d\n",meson->GetName(), mesonIndex);
843 //Check if meson decayed into 2 daughters or if both were kept.
844 if(meson->GetNDaughters() != 2)
847 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
852 Int_t iPhoton0 = meson->GetDaughter(0);
853 Int_t iPhoton1 = meson->GetDaughter(1);
854 TParticle *photon0 = stack->Particle(iPhoton0);
855 TParticle *photon1 = stack->Particle(iPhoton1);
857 //Check if both daughters are photons
858 if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22)
861 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
866 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
868 //Check if both photons contribute to the cluster
869 Bool_t okPhoton0 = kFALSE;
870 Bool_t okPhoton1 = kFALSE;
872 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Labels loop:\n");
874 Bool_t conversion = kFALSE;
876 for(Int_t i = 0; i < nlabels; i++)
878 if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
880 //If we already found both, break the loop
881 if(okPhoton0 && okPhoton1) break;
883 Int_t index = labels[i];
884 if (iPhoton0 == index)
889 else if (iPhoton1 == index)
895 //Trace back the mother in case it was a conversion
897 if(index >= stack->GetNtrack())
899 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",
900 index,stack->GetNtrack());
904 TParticle * daught = stack->Particle(index);
905 Int_t tmpindex = daught->GetFirstMother();
906 if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
909 //MC particle of interest is the mother
910 if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
911 daught = stack->Particle(tmpindex);
912 if (iPhoton0 == tmpindex)
918 else if (iPhoton1 == tmpindex)
925 tmpindex = daught->GetFirstMother();
927 }//While to check if pi0/eta daughter was one of these contributors to the cluster
929 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
930 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Something happens, first label should be from a photon decay!\n");
932 }//loop on list of labels
934 //If both photons contribute tag as the corresponding meson.
935 if(okPhoton0 && okPhoton1)
938 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s OVERLAPPED DECAY \n", meson->GetName());
940 if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ;
942 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
943 else SetTagBit(tag,kMCEta);
948 //________________________________________________________________________________________________________
949 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex,
950 const TClonesArray *mcparticles, Int_t & tag )
952 //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
954 if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1)
957 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
958 labels[0],mcparticles->GetEntriesFast(), nlabels);
962 AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
963 Int_t mesonPdg = meson->GetPdgCode();
964 if(mesonPdg != 111 && mesonPdg != 221)
966 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
970 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
974 if(meson->GetNDaughters() != 2)
977 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
981 Int_t iPhoton0 = meson->GetDaughter(0);
982 Int_t iPhoton1 = meson->GetDaughter(1);
983 //if((iPhoton0 == -1) || (iPhoton1 == -1)){
985 // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
988 AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
989 AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
991 //Check if both daughters are photons
992 if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22)
994 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overlapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
999 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
1001 //Check if both photons contribute to the cluster
1002 Bool_t okPhoton0 = kFALSE;
1003 Bool_t okPhoton1 = kFALSE;
1006 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
1008 Bool_t conversion = kFALSE;
1010 for(Int_t i = 0; i < nlabels; i++)
1013 printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
1015 if(labels[i]<0) continue;
1017 //If we already found both, break the loop
1018 if(okPhoton0 && okPhoton1) break;
1020 Int_t index = labels[i];
1021 if (iPhoton0 == index)
1026 else if (iPhoton1 == index)
1032 //Trace back the mother in case it was a conversion
1034 if(index >= mcparticles->GetEntriesFast())
1036 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast());
1040 AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
1041 Int_t tmpindex = daught->GetMother();
1043 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
1047 //MC particle of interest is the mother
1049 printf("\t parent index %d\n",tmpindex);
1050 daught = (AliAODMCParticle*) mcparticles->At(tmpindex);
1051 //printf("tmpindex %d\n",tmpindex);
1052 if (iPhoton0 == tmpindex)
1058 else if (iPhoton1 == tmpindex)
1065 tmpindex = daught->GetMother();
1067 }//While to check if pi0/eta daughter was one of these contributors to the cluster
1069 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
1070 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
1072 }//loop on list of labels
1074 //If both photons contribute tag as the corresponding meson.
1075 if(okPhoton0 && okPhoton1)
1078 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",
1079 (TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
1081 if(!CheckTagBit(tag,kMCConversion) && conversion)
1084 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Second decay photon produced a conversion \n");
1085 SetTagBit(tag,kMCConversion) ;
1088 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
1089 else SetTagBit(tag,kMCEta);
1094 //_____________________________________________________________________
1095 TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
1097 //Return list of jets (TParticles) and index of most likely parton that originated it.
1098 AliStack * stack = reader->GetStack();
1099 Int_t iEvent = reader->GetEventNumber();
1100 AliGenEventHeader * geh = reader->GetGenEventHeader();
1101 if(fCurrentEvent!=iEvent){
1102 fCurrentEvent = iEvent;
1103 fJetsList = new TList;
1104 Int_t nTriggerJets = 0;
1105 Float_t tmpjet[]={0,0,0,0};
1107 //printf("Event %d %d\n",fCurrentEvent,iEvent);
1108 //Get outgoing partons
1109 if(stack->GetNtrack() < 8) return fJetsList;
1110 TParticle * parton1 = stack->Particle(6);
1111 TParticle * parton2 = stack->Particle(7);
1113 printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1114 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
1115 printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1116 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
1118 // //Trace the jet from the mother parton
1125 // TParticle * tmptmp = new TParticle;
1126 // for(Int_t i = 0; i< stack->GetNprimary(); i++){
1127 // tmptmp = stack->Particle(i);
1129 // if(tmptmp->GetStatusCode() == 1){
1130 // pt = tmptmp->Pt();
1131 // e = tmptmp->Energy();
1132 // Int_t imom = tmptmp->GetFirstMother();
1134 // //printf("1st imom %d\n",imom);
1137 // tmptmp = stack->Particle(imom);
1138 // imom = tmptmp->GetFirstMother();
1139 // //printf("imom %d \n",imom);
1141 // //printf("Last imom %d %d\n",imom1, imom);
1146 // else if (imom1 == 7){
1153 // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
1155 //Get the jet, different way for different generator
1157 if(fMCGenerator == "PYTHIA")
1159 TParticle * jet = 0x0;
1160 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
1161 nTriggerJets = pygeh->NTriggerJets();
1163 printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
1165 for(Int_t i = 0; i< nTriggerJets; i++)
1167 pygeh->TriggerJet(i, tmpjet);
1168 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1169 //Assign an outgoing parton as mother
1170 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1171 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1172 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1173 else jet->SetFirstMother(6);
1176 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1177 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1178 fJetsList->Add(jet);
1180 }//Pythia triggered jets
1182 else if (fMCGenerator=="HERWIG")
1186 TParticle * tmp = parton1;
1187 if(parton1->GetPdgCode()!=22)
1190 if(tmp->GetFirstDaughter()==-1) return fJetsList;
1191 tmp = stack->Particle(tmp->GetFirstDaughter());
1192 pdg = tmp->GetPdgCode();
1195 //Add found jet to list
1196 TParticle *jet1 = new TParticle(*tmp);
1197 jet1->SetFirstMother(6);
1198 fJetsList->Add(jet1);
1199 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1200 //tmp = stack->Particle(tmp->GetFirstDaughter());
1204 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1205 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
1211 if(parton2->GetPdgCode()!=22)
1215 if(tmp->GetFirstDaughter()==-1) return fJetsList;
1216 tmp = stack->Particle(tmp->GetFirstDaughter());
1217 pdg = tmp->GetPdgCode();
1220 //Add found jet to list
1221 TParticle *jet2 = new TParticle(*tmp);
1222 jet2->SetFirstMother(7);
1223 fJetsList->Add(jet2);
1226 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1227 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1228 //Int_t first = tmp->GetFirstDaughter();
1229 //Int_t last = tmp->GetLastDaughter();
1230 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1231 // for(Int_t d = first ; d < last+1; d++){
1232 // tmp = stack->Particle(d);
1233 // if(i == tmp->GetFirstMother())
1234 // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1235 // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1239 }//Herwig generated jets
1246 //__________________________________________________________________________________________________________
1247 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
1248 const AliCaloTrackReader* reader,
1249 Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel)
1251 //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1253 TLorentzVector daughter(0,0,0,0);
1255 if(reader->ReadStack())
1257 if(!reader->GetStack())
1260 printf("AliMCAnalysisUtils::GetDaughter() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1265 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1267 TParticle * momP = reader->GetStack()->Particle(label);
1268 daughlabel = momP->GetDaughter(idaugh);
1270 TParticle * daughP = reader->GetStack()->Particle(daughlabel);
1271 daughP->Momentum(daughter);
1272 pdg = daughP->GetPdgCode();
1273 status = daughP->GetStatusCode();
1281 else if(reader->ReadAODMCParticles())
1283 TClonesArray* mcparticles = reader->GetAODMCParticles();
1287 printf("AliMCAnalysisUtils::GetDaughter() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1293 Int_t nprimaries = mcparticles->GetEntriesFast();
1294 if(label >= 0 && label < nprimaries)
1296 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1297 daughlabel = momP->GetDaughter(idaugh);
1299 AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
1300 daughter.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1301 pdg = daughP->GetPdgCode();
1302 status = daughP->GetStatus();
1316 //______________________________________________________________________________________________________
1317 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1319 //Return the kinematics of the particle that generated the signal
1321 Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1322 return GetMother(label,reader,pdg,status, ok,momlabel);
1325 //_________________________________________________________________________________________
1326 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1327 Int_t & pdg, Int_t & status, Bool_t & ok)
1329 //Return the kinematics of the particle that generated the signal
1331 Int_t momlabel = -1;
1332 return GetMother(label,reader,pdg,status, ok,momlabel);
1335 //______________________________________________________________________________________________________
1336 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1337 Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1339 //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1341 TLorentzVector mom(0,0,0,0);
1343 if(reader->ReadStack())
1345 if(!reader->GetStack())
1348 printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1353 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1355 TParticle * momP = reader->GetStack()->Particle(label);
1356 momP->Momentum(mom);
1357 pdg = momP->GetPdgCode();
1358 status = momP->GetStatusCode();
1359 momlabel = momP->GetFirstMother();
1367 else if(reader->ReadAODMCParticles())
1369 TClonesArray* mcparticles = reader->GetAODMCParticles();
1373 printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1379 Int_t nprimaries = mcparticles->GetEntriesFast();
1380 if(label >= 0 && label < nprimaries)
1382 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1383 mom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1384 pdg = momP->GetPdgCode();
1385 status = momP->GetStatus();
1386 momlabel = momP->GetMother();
1400 //___________________________________________________________________________________
1401 TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1402 const AliCaloTrackReader* reader,
1403 Bool_t & ok, Int_t & momlabel)
1405 //Return the kinematics of the particle that generated the signal
1407 TLorentzVector grandmom(0,0,0,0);
1410 if(reader->ReadStack())
1412 if(!reader->GetStack())
1415 printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1421 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1423 TParticle * momP = reader->GetStack()->Particle(label);
1425 Int_t grandmomLabel = momP->GetFirstMother();
1426 Int_t grandmomPDG = -1;
1427 TParticle * grandmomP = 0x0;
1428 while (grandmomLabel >=0 )
1430 grandmomP = reader->GetStack()->Particle(grandmomLabel);
1431 grandmomPDG = grandmomP->GetPdgCode();
1432 if(grandmomPDG==pdg)
1434 momlabel = grandmomLabel;
1435 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1439 grandmomLabel = grandmomP->GetFirstMother();
1443 if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1446 else if(reader->ReadAODMCParticles())
1448 TClonesArray* mcparticles = reader->GetAODMCParticles();
1452 printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1458 Int_t nprimaries = mcparticles->GetEntriesFast();
1459 if(label >= 0 && label < nprimaries)
1461 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1463 Int_t grandmomLabel = momP->GetMother();
1464 Int_t grandmomPDG = -1;
1465 AliAODMCParticle * grandmomP = 0x0;
1466 while (grandmomLabel >=0 )
1468 grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1469 grandmomPDG = grandmomP->GetPdgCode();
1470 if(grandmomPDG==pdg)
1472 //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1473 momlabel = grandmomLabel;
1474 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1478 grandmomLabel = grandmomP->GetMother();
1482 if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
1492 //______________________________________________________________________________________________
1493 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
1494 Int_t & pdg, Int_t & status, Bool_t & ok,
1495 Int_t & grandMomLabel, Int_t & greatMomLabel)
1497 //Return the kinematics of the particle that generated the signal
1499 TLorentzVector grandmom(0,0,0,0);
1501 if(reader->ReadStack())
1503 if(!reader->GetStack())
1506 printf("AliMCAnalysisUtils::GetGrandMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1512 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1514 TParticle * momP = reader->GetStack()->Particle(label);
1516 grandMomLabel = momP->GetFirstMother();
1518 TParticle * grandmomP = 0x0;
1520 if (grandMomLabel >=0 )
1522 grandmomP = reader->GetStack()->Particle(grandMomLabel);
1523 pdg = grandmomP->GetPdgCode();
1524 status = grandmomP->GetStatusCode();
1526 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1527 greatMomLabel = grandmomP->GetFirstMother();
1532 else if(reader->ReadAODMCParticles())
1534 TClonesArray* mcparticles = reader->GetAODMCParticles();
1538 printf("AliMCAnalysisUtils::GetGrandMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1544 Int_t nprimaries = mcparticles->GetEntriesFast();
1545 if(label >= 0 && label < nprimaries)
1547 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1549 grandMomLabel = momP->GetMother();
1551 AliAODMCParticle * grandmomP = 0x0;
1553 if(grandMomLabel >=0 )
1555 grandmomP = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1556 pdg = grandmomP->GetPdgCode();
1557 status = grandmomP->GetStatus();
1559 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1560 greatMomLabel = grandmomP->GetMother();
1571 //_______________________________________________________________________________________________________________
1572 void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader,
1573 Float_t & asym, Float_t & angle, Bool_t & ok)
1575 //In case of an eta or pi0 decay into 2 photons, get the asymmetry in the energy of the photons
1578 if(reader->ReadStack())
1580 if(!reader->GetStack())
1583 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1587 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1589 TParticle * momP = reader->GetStack()->Particle(label);
1591 Int_t grandmomLabel = momP->GetFirstMother();
1592 Int_t grandmomPDG = -1;
1593 TParticle * grandmomP = 0x0;
1594 while (grandmomLabel >=0 )
1596 grandmomP = reader->GetStack()->Particle(grandmomLabel);
1597 grandmomPDG = grandmomP->GetPdgCode();
1599 if(grandmomPDG==pdg) break;
1601 grandmomLabel = grandmomP->GetFirstMother();
1605 if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1607 TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1608 TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1609 if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1611 asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
1612 TLorentzVector lvd1, lvd2;
1615 angle = lvd1.Angle(lvd2.Vect());
1621 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1626 else if(reader->ReadAODMCParticles())
1628 TClonesArray* mcparticles = reader->GetAODMCParticles();
1632 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1638 Int_t nprimaries = mcparticles->GetEntriesFast();
1639 if(label >= 0 && label < nprimaries)
1641 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1643 Int_t grandmomLabel = momP->GetMother();
1644 Int_t grandmomPDG = -1;
1645 AliAODMCParticle * grandmomP = 0x0;
1646 while (grandmomLabel >=0 )
1648 grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1649 grandmomPDG = grandmomP->GetPdgCode();
1651 if(grandmomPDG==pdg) break;
1653 grandmomLabel = grandmomP->GetMother();
1657 if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1659 AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1660 AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1661 if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1663 asym = (d1->E()-d2->E())/grandmomP->E();
1664 TLorentzVector lvd1, lvd2;
1665 lvd1.SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1666 lvd2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1667 angle = lvd1.Angle(lvd2.Vect());
1673 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg);
1684 //_________________________________________________________________________________________________
1685 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1687 // Return the the number of daughters of a given MC particle
1690 if(reader->ReadStack())
1692 if(!reader->GetStack())
1695 printf("AliMCAnalysisUtils::GetNDaughters() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1700 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1702 TParticle * momP = reader->GetStack()->Particle(label);
1704 return momP->GetNDaughters();
1712 else if(reader->ReadAODMCParticles())
1714 TClonesArray* mcparticles = reader->GetAODMCParticles();
1718 printf("AliMCAnalysisUtils::GetNDaughters() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1724 Int_t nprimaries = mcparticles->GetEntriesFast();
1725 if(label >= 0 && label < nprimaries)
1727 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1729 return momP->GetNDaughters();
1743 //_________________________________________________________________________________
1744 Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
1745 Int_t mctag, Int_t mesonLabel,
1746 AliCaloTrackReader * reader, Int_t *overpdg)
1748 // Compare the primary depositing more energy with the rest,
1749 // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing
1750 // Give as input the meson label in case it was a pi0 or eta merged cluster
1751 // Init overpdg with nlabels
1753 Int_t ancPDG = 0, ancStatus = -1;
1754 TLorentzVector momentum; TVector3 prodVertex;
1756 Int_t noverlaps = 0;
1759 for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1761 ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,momentum,prodVertex);
1763 //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1764 // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1766 Bool_t overlap = kFALSE;
1771 //printf("\t \t \t No Label = %d\n",ancLabel);
1773 else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
1775 //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1778 else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1780 //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1784 if( !overlap ) continue ;
1786 // We have at least one overlap
1788 //printf("Overlap!!!!!!!!!!!!!!\n");
1792 // What is the origin of the overlap?
1793 Bool_t mOK = 0, gOK = 0;
1794 Int_t mpdg = -999999, gpdg = -1;
1795 Int_t mstatus = -1, gstatus = -1;
1796 Int_t gLabel = -1, ggLabel = -1;
1797 TLorentzVector mother = GetMother (label[ilab],reader,mpdg,mstatus,mOK);
1798 TLorentzVector grandmother = GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
1800 //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
1802 if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
1803 ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
1806 Int_t labeltmp = gLabel;
1807 while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
1810 grandmother = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
1814 overpdg[noverlaps-1] = mpdg;
1821 //________________________________________________________
1822 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1824 //Print some relevant parameters set for the analysis
1829 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1831 printf("Debug level = %d\n",fDebug);
1832 printf("MC Generator = %s\n",fMCGenerator.Data());
1837 //__________________________________________________
1838 void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
1840 // print the assigned origins to this particle
1842 printf("AliMCAnalysisUtils::PrintMCTag() - tag %d \n photon %d, conv %d, prompt %d, frag %d, isr %d, \n pi0 decay %d, eta decay %d, other decay %d pi0 %d, eta %d \n electron %d, muon %d,pion %d, proton %d, neutron %d, \n kaon %d, a-proton %d, a-neutron %d, other %d, unk %d, bad %d\n",
1844 CheckTagBit(tag,kMCPhoton),
1845 CheckTagBit(tag,kMCConversion),
1846 CheckTagBit(tag,kMCPrompt),
1847 CheckTagBit(tag,kMCFragmentation),
1848 CheckTagBit(tag,kMCISR),
1849 CheckTagBit(tag,kMCPi0Decay),
1850 CheckTagBit(tag,kMCEtaDecay),
1851 CheckTagBit(tag,kMCOtherDecay),
1852 CheckTagBit(tag,kMCPi0),
1853 CheckTagBit(tag,kMCEta),
1854 CheckTagBit(tag,kMCElectron),
1855 CheckTagBit(tag,kMCMuon),
1856 CheckTagBit(tag,kMCPion),
1857 CheckTagBit(tag,kMCProton),
1858 CheckTagBit(tag,kMCAntiNeutron),
1859 CheckTagBit(tag,kMCKaon),
1860 CheckTagBit(tag,kMCAntiProton),
1861 CheckTagBit(tag,kMCAntiNeutron),
1862 CheckTagBit(tag,kMCOther),
1863 CheckTagBit(tag,kMCUnknown),
1864 CheckTagBit(tag,kMCBadLabel)