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 **************************************************************************/
15 /* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
17 //_________________________________________________________________________
18 // Class for analysis utils for MC data
19 // stored in stack or event header.
21 // - method to check the origin of a given track/cluster
22 // - method to obtain the generated jets
24 //*-- Author: Gustavo Conesa (LNF-INFN)
25 //////////////////////////////////////////////////////////////////////////////
28 // --- ROOT system ---
31 #include "TParticle.h"
32 #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() :
45 TObject(), fCurrentEvent(-1), fDebug(-1),
46 fJetsList(new TList), fMCGenerator("PYTHIA")
51 //____________________________________________________________________________
52 AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
53 TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
54 fJetsList(new TList), fMCGenerator(mcutils.fMCGenerator)
61 //_________________________________________________________________________
62 //AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
64 // // assignment operator
66 // if(&mcutils == this) return *this;
67 // fCurrentEvent = mcutils.fCurrentEvent ;
68 // fDebug = mcutils.fDebug;
69 // fJetsList = mcutils.fJetsList;
70 // fMCGenerator = mcutils.fMCGenerator;
75 //____________________________________________________________________________
76 AliMCAnalysisUtils::~AliMCAnalysisUtils()
78 // Remove all pointers.
87 //_________________________________________________________________________
88 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, const Int_t nlabels, AliCaloTrackReader* reader, const Int_t input = 0) {
89 //Play with the montecarlo particles if available
93 printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
97 //Select where the information is, ESD-galice stack or AOD mcparticles branch
98 if(reader->ReadStack()){
99 tag = CheckOriginInStack(label, nlabels, reader->GetStack());
101 else if(reader->ReadAODMCParticles()){
102 tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(input));
108 //_________________________________________________________________________
109 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) {
110 //Play with the montecarlo particles if available
114 printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
118 Int_t labels[]={label};
120 //Select where the information is, ESD-galice stack or AOD mcparticles branch
121 if(reader->ReadStack()){
122 tag = CheckOriginInStack(labels, 1,reader->GetStack());
124 else if(reader->ReadAODMCParticles()){
125 tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(input));
131 //_________________________________________________________________________
132 Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, const Int_t nlabels, AliStack* stack) {
133 // Play with the MC stack if available. Tag particles depending on their origin.
134 // Do same things as in CheckOriginInAOD but different input.
136 //generally speaking, label is the MC label of a reconstructed
137 //entity (track, cluster, etc) for which we want to know something
138 //about its heritage, but one can also use it directly with stack
139 //particles not connected to reconstructed entities
143 printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
148 Int_t label=labels[0];//Most significant particle contributing to the cluster
150 if(label >= 0 && label < stack->GetNtrack()){
151 //MC particle of interest is the "mom" of the entity
152 TParticle * mom = stack->Particle(label);
154 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
155 Int_t mStatus = mom->GetStatusCode() ;
156 Int_t iParent = mom->GetFirstMother() ;
157 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
159 //GrandParent of the entity
160 TParticle * parent = NULL;
164 parent = stack->Particle(iParent);
166 pPdg = TMath::Abs(parent->GetPdgCode());
167 pStatus = parent->GetStatusCode();
170 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
173 printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
174 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
175 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
178 //Check if "mother" of entity is converted, if not, get the first non converted mother
179 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
180 SetTagBit(tag,kMCConversion);
181 //Check if the mother is photon or electron with status not stable
182 while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
184 iMom = mom->GetFirstMother();
185 mom = stack->Particle(iMom);
186 mPdg = TMath::Abs(mom->GetPdgCode());
187 mStatus = mom->GetStatusCode() ;
188 iParent = mom->GetFirstMother() ;
189 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
193 parent = stack->Particle(iParent);
195 pPdg = TMath::Abs(parent->GetPdgCode());
196 pStatus = parent->GetStatusCode();
199 else {// in case of gun/box simulations
206 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
207 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
208 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
211 }//mother and parent are electron or photon and have status 0
212 else if((mPdg == 22 || mPdg == 11) && mStatus == 0){
213 //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
214 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
215 pPdg == 2212 || pPdg == 130 || pPdg == 13 ) {
216 SetTagBit(tag,kMCConversion);
217 iMom = mom->GetFirstMother();
218 mom = stack->Particle(iMom);
219 mPdg = TMath::Abs(mom->GetPdgCode());
222 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
223 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
227 //Comment for the next lines, we do not check the parent of the hadron for the moment.
228 //iParent = mom->GetFirstMother() ;
229 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
233 // parent = stack->Particle(iParent);
234 // pPdg = TMath::Abs(parent->GetPdgCode());
237 // conversion into electrons/photons checked
239 //first check for typical charged particles
240 if(mPdg == 13) SetTagBit(tag,kMCMuon);
241 else if(mPdg == 211) SetTagBit(tag,kMCPion);
242 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
243 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
244 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
245 else if(mPdg == 111) {
246 SetTagBit(tag,kMCPi0Decay);
247 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
248 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
250 else if(mPdg == 221) {
251 SetTagBit(tag,kMCEtaDecay);
252 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
253 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
257 SetTagBit(tag,kMCPhoton);
258 if(mStatus == 1){ //undecayed particle
259 if(fMCGenerator == "PYTHIA"){
260 if(iParent < 8 && iParent > 5) {//outgoing partons
261 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
262 else SetTagBit(tag,kMCFragmentation);
264 else if(iParent <= 5) {
265 SetTagBit(tag, kMCISR); //Initial state radiation
267 else if(pStatus == 11){//Decay
269 SetTagBit(tag,kMCPi0Decay);
270 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
271 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
273 else if (pPdg == 221) {
274 SetTagBit(tag, kMCEtaDecay);
275 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
276 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
278 else SetTagBit(tag,kMCOtherDecay);
281 if(fDebug > 1 && parent) printf("AliMCAnalysisUtils::CheckOrigingInStack() - what is it in PYTHIA? Wrong generator setting? Mother mPdg %d, status %d \n Parent iParent %d, pPdg %d %s, status %d\n",
282 mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
284 SetTagBit(tag,kMCPi0Decay);
285 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
286 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
288 else if (pPdg == 221) {
289 SetTagBit(tag, kMCEtaDecay);
290 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
291 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
293 else SetTagBit(tag,kMCOtherDecay);
297 else if(fMCGenerator == "HERWIG"){
298 if(pStatus < 197){//Not decay
301 if(parent->GetFirstMother()<=5) break;
302 iParent = parent->GetFirstMother();
303 parent=stack->Particle(iParent);
304 pStatus= parent->GetStatusCode();
305 pPdg = TMath::Abs(parent->GetPdgCode());
307 }//Look for the parton
309 if(iParent < 8 && iParent > 5) {
310 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
311 else SetTagBit(tag,kMCFragmentation);
313 else SetTagBit(tag,kMCISR);//Initial state radiation
317 SetTagBit(tag,kMCPi0Decay);
318 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
319 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
321 else if (pPdg == 221) {
322 SetTagBit(tag,kMCEtaDecay);
323 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
324 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
326 else SetTagBit(tag,kMCOtherDecay);
330 else SetTagBit(tag,kMCUnknown);
332 }//Status 1 : created by event generator
334 else if(mStatus == 0){ // geant
336 SetTagBit(tag,kMCPi0Decay);
337 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
338 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
340 else if (pPdg == 221) {
341 SetTagBit(tag,kMCEtaDecay);
342 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
343 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
345 else SetTagBit(tag,kMCOtherDecay);
346 }//status 0 : geant generated
350 //Electron check. Where did that electron come from?
351 else if(mPdg == 11){ //electron
352 if(pPdg == 11 && parent){
353 Int_t iGrandma = parent->GetFirstMother();
355 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
356 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
358 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
359 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
362 SetTagBit(tag,kMCElectron);
363 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
364 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
365 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
366 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
367 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
369 Int_t iGrandma = parent->GetFirstMother();
371 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
372 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
373 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
374 else SetTagBit(tag,kMCEFromC); //c-->e
375 } else SetTagBit(tag,kMCEFromC); //c-->e
378 //if it is not from any of the above, where is it from?
379 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
380 else SetTagBit(tag,kMCOtherDecay);
381 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);
384 //Cluster was made by something else
386 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
387 SetTagBit(tag,kMCUnknown);
392 if(label < 0 && (fDebug >= 0))
393 printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
394 if(label >= stack->GetNtrack() && (fDebug >= 0))
395 printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
396 SetTagBit(tag,kMCUnknown);
404 //_________________________________________________________________________
405 Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, const Int_t nlabels, TClonesArray *mcparticles) {
406 // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
407 // Do same things as in CheckOriginInStack but different input.
410 printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
415 Int_t label=labels[0];//Most significant particle contributing to the cluster
417 Int_t nprimaries = mcparticles->GetEntriesFast();
418 if(label >= 0 && label < nprimaries){
420 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
422 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
423 Int_t iParent = mom->GetMother() ;
424 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
427 AliAODMCParticle * parent = NULL ;
430 parent = (AliAODMCParticle *) mcparticles->At(iParent);
431 pPdg = TMath::Abs(parent->GetPdgCode());
433 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
436 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
437 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
439 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
442 //Check if mother is converted, if not, get the first non converted mother
443 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
444 SetTagBit(tag,kMCConversion);
445 //Check if the mother is photon or electron with status not stable
446 while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
448 iMom = mom->GetMother();
449 mom = (AliAODMCParticle *) mcparticles->At(iMom);
450 mPdg = TMath::Abs(mom->GetPdgCode());
451 iParent = mom->GetMother() ;
452 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
455 if(iParent >= 0 && parent){
456 parent = (AliAODMCParticle *) mcparticles->At(iParent);
457 pPdg = TMath::Abs(parent->GetPdgCode());
459 // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
460 // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
465 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
466 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
468 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
471 }//mother and parent are electron or photon and have status 0 and parent is photon or electron
472 else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){
473 //Still a conversion but only one electron/photon generated. Just from hadrons
474 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
475 pPdg == 2212 || pPdg == 130 || pPdg == 13 ) {
476 SetTagBit(tag,kMCConversion);
477 iMom = mom->GetMother();
478 mom = (AliAODMCParticle *) mcparticles->At(iMom);
479 mPdg = TMath::Abs(mom->GetPdgCode());
482 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
483 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
487 //Comment for next lines, we do not check the parent of the hadron for the moment.
488 //iParent = mom->GetMother() ;
489 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
493 // parent = (AliAODMCParticle *) mcparticles->At(iParent);
494 // pPdg = TMath::Abs(parent->GetPdgCode());
498 // conversion into electrons/photons checked
500 //first check for typical charged particles
501 if(mPdg == 13) SetTagBit(tag,kMCMuon);
502 else if(mPdg == 211) SetTagBit(tag,kMCPion);
503 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
504 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
505 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
506 else if(mPdg == 111) {
507 SetTagBit(tag,kMCPi0Decay);
508 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
509 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
511 else if(mPdg == 221) {
512 SetTagBit(tag,kMCEtaDecay);
513 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
514 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
518 SetTagBit(tag,kMCPhoton);
519 if(mom->IsPhysicalPrimary()){ //undecayed particle
520 if(iParent < 8 && iParent > 5) {//outgoing partons
521 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
522 else SetTagBit(tag,kMCFragmentation);
524 else if(iParent <= 5) {
525 SetTagBit(tag, kMCISR); //Initial state radiation
527 else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
529 SetTagBit(tag,kMCPi0Decay);
530 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
531 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
533 else if (pPdg == 221) {
534 SetTagBit(tag, kMCEtaDecay);
535 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
536 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
538 else SetTagBit(tag,kMCOtherDecay);
541 if(parent)printf("AliMCAnalysisUtils::ChecOrigingInAOD() - what is it? Mother mPdg %d, is primary? %d, is physical %d \n Parent iParent %d, pPdg %d, is primary? %d, is physical? %d\n",
542 mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
543 SetTagBit(tag,kMCOtherDecay);//Check
546 else if(!mom->IsPrimary()){ //Decays
548 SetTagBit(tag,kMCPi0Decay);
549 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
550 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
552 else if (pPdg == 221) {
553 SetTagBit(tag,kMCEtaDecay);
554 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n");
555 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
557 else SetTagBit(tag,kMCOtherDecay);
558 }//not primary : geant generated, decays
560 //printf("UNKNOWN 1, mom pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
561 //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
562 SetTagBit(tag,kMCUnknown);
566 //Electron check. Where did that electron come from?
567 else if(mPdg == 11){ //electron
568 if(pPdg == 11 && parent){
569 Int_t iGrandma = parent->GetMother();
571 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
572 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
574 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
575 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
578 SetTagBit(tag,kMCElectron);
579 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
580 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
581 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
582 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
583 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
585 Int_t iGrandma = parent->GetMother();
587 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
588 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
589 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
590 else SetTagBit(tag,kMCEFromC); //c-hadron decay
591 } else SetTagBit(tag,kMCEFromC); //c-hadron decay
593 } else { //prompt or other decay
594 TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
595 TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
596 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg);
597 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
598 else SetTagBit(tag,kMCOtherDecay);
601 //cluster was made by something else
603 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg);
604 SetTagBit(tag,kMCUnknown);
609 if(label < 0 && (fDebug >= 0) )
610 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label);
611 if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) )
612 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
613 SetTagBit(tag,kMCUnknown);
621 //_________________________________________________________________________
622 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
623 AliStack *stack, Int_t &tag) {
624 //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
626 if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
627 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
628 labels[0],stack->GetNtrack(), nlabels);
632 TParticle * meson = stack->Particle(mesonIndex);
633 Int_t mesonPdg = meson->GetPdgCode();
634 if(mesonPdg!=111 && mesonPdg!=221){
635 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
639 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
641 //Check if meson decayed into 2 daughters or if both were kept.
642 if(meson->GetNDaughters() != 2){
644 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
649 Int_t iPhoton0 = meson->GetDaughter(0);
650 Int_t iPhoton1 = meson->GetDaughter(1);
651 TParticle *photon0 = stack->Particle(iPhoton0);
652 TParticle *photon1 = stack->Particle(iPhoton1);
654 //Check if both daughters are photons
655 if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
657 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
662 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
664 //Check if both photons contribute to the cluster
665 Bool_t okPhoton0 = kFALSE;
666 Bool_t okPhoton1 = kFALSE;
668 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
670 for(Int_t i = 0; i < nlabels; i++){
671 if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
673 //If we already found both, break the loop
674 if(okPhoton0 && okPhoton1) break;
676 Int_t index = labels[i];
677 if (iPhoton0 == index) {
681 else if (iPhoton1 == index) {
686 //Trace back the mother in case it was a conversion
687 TParticle * daught = stack->Particle(index);
688 Int_t tmpindex = daught->GetFirstMother();
689 if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
691 //MC particle of interest is the mother
692 if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
693 daught = stack->Particle(tmpindex);
694 if (iPhoton0 == tmpindex) {
698 else if (iPhoton1 == tmpindex) {
702 tmpindex = daught->GetFirstMother();
703 }//While to check if pi0/eta daughter was one of these contributors to the cluster
705 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
706 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
708 }//loop on list of labels
710 //If both photons contribute tag as the corresponding meson.
711 if(okPhoton0 && okPhoton1){
713 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
715 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
716 else SetTagBit(tag,kMCEta);
721 //_________________________________________________________________________
722 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
723 TClonesArray *mcparticles, Int_t &tag) {
724 //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
726 if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1) {
727 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
728 labels[0],mcparticles->GetEntriesFast(), nlabels);
733 AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
734 Int_t mesonPdg = meson->GetPdgCode();
735 if(mesonPdg != 111 && mesonPdg != 221) {
736 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
740 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
744 if(meson->GetNDaughters() != 2){
746 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
749 Int_t iPhoton0 = meson->GetDaughter(0);
750 Int_t iPhoton1 = meson->GetDaughter(1);
751 //if((iPhoton0 == -1) || (iPhoton1 == -1)){
753 // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overalapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
756 AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
757 AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
759 //Check if both daughters are photons
760 if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22){
761 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
766 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
768 //Check if both photons contribute to the cluster
769 Bool_t okPhoton0 = kFALSE;
770 Bool_t okPhoton1 = kFALSE;
772 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
774 for(Int_t i = 0; i < nlabels; i++){
775 if(fDebug > 3) printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
777 //If we already found both, break the loop
778 if(okPhoton0 && okPhoton1) break;
780 Int_t index = labels[i];
781 if (iPhoton0 == index) {
785 else if (iPhoton1 == index) {
790 //Trace back the mother in case it was a conversion
791 AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
792 Int_t tmpindex = daught->GetMother();
793 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
797 //MC particle of interest is the mother
798 if(fDebug > 3) printf("\t parent index %d\n",tmpindex);
799 daught = (AliAODMCParticle*) mcparticles->At(tmpindex);
800 //printf("tmpindex %d\n",tmpindex);
801 if (iPhoton0 == tmpindex) {
805 else if (iPhoton1 == tmpindex) {
809 tmpindex = daught->GetMother();
810 }//While to check if pi0/eta daughter was one of these contributors to the cluster
812 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
813 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
815 }//loop on list of labels
817 //If both photons contribute tag as the corresponding meson.
818 if(okPhoton0 && okPhoton1){
820 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
822 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
823 else SetTagBit(tag,kMCEta);
828 //_________________________________________________________________________
829 TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * const reader){
830 //Return list of jets (TParticles) and index of most likely parton that originated it.
831 AliStack * stack = reader->GetStack();
832 Int_t iEvent = reader->GetEventNumber();
833 AliGenEventHeader * geh = reader->GetGenEventHeader();
834 if(fCurrentEvent!=iEvent){
835 fCurrentEvent = iEvent;
836 fJetsList = new TList;
837 Int_t nTriggerJets = 0;
838 Float_t tmpjet[]={0,0,0,0};
840 //printf("Event %d %d\n",fCurrentEvent,iEvent);
841 //Get outgoing partons
842 if(stack->GetNtrack() < 8) return fJetsList;
843 TParticle * parton1 = stack->Particle(6);
844 TParticle * parton2 = stack->Particle(7);
846 printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
847 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
848 printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
849 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
851 // //Trace the jet from the mother parton
858 // TParticle * tmptmp = new TParticle;
859 // for(Int_t i = 0; i< stack->GetNprimary(); i++){
860 // tmptmp = stack->Particle(i);
862 // if(tmptmp->GetStatusCode() == 1){
863 // pt = tmptmp->Pt();
864 // e = tmptmp->Energy();
865 // Int_t imom = tmptmp->GetFirstMother();
867 // //printf("1st imom %d\n",imom);
870 // tmptmp = stack->Particle(imom);
871 // imom = tmptmp->GetFirstMother();
872 // //printf("imom %d \n",imom);
874 // //printf("Last imom %d %d\n",imom1, imom);
879 // else if (imom1 == 7){
886 // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
888 //Get the jet, different way for different generator
890 if(fMCGenerator == "PYTHIA"){
891 TParticle * jet = 0x0;
892 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
893 nTriggerJets = pygeh->NTriggerJets();
895 printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
898 for(Int_t i = 0; i< nTriggerJets; i++){
900 pygeh->TriggerJet(i, tmpjet);
901 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
902 //Assign an outgoing parton as mother
903 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
904 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
905 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
906 else jet->SetFirstMother(6);
909 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
910 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
913 }//Pythia triggered jets
915 else if (fMCGenerator=="HERWIG"){
918 TParticle * tmp = parton1;
919 if(parton1->GetPdgCode()!=22){
921 if(tmp->GetFirstDaughter()==-1) return fJetsList;
922 tmp = stack->Particle(tmp->GetFirstDaughter());
923 pdg = tmp->GetPdgCode();
926 //Add found jet to list
927 TParticle *jet1 = new TParticle(*tmp);
928 jet1->SetFirstMother(6);
929 fJetsList->Add(jet1);
930 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
931 //tmp = stack->Particle(tmp->GetFirstDaughter());
935 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
936 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
943 if(parton2->GetPdgCode()!=22){
945 if(tmp->GetFirstDaughter()==-1) return fJetsList;
946 i = tmp->GetFirstDaughter();
947 tmp = stack->Particle(tmp->GetFirstDaughter());
948 pdg = tmp->GetPdgCode();
950 //Add found jet to list
951 TParticle *jet2 = new TParticle(*tmp);
952 jet2->SetFirstMother(7);
953 fJetsList->Add(jet2);
956 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
957 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
958 //Int_t first = tmp->GetFirstDaughter();
959 //Int_t last = tmp->GetLastDaughter();
960 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
961 // for(Int_t d = first ; d < last+1; d++){
962 // tmp = stack->Particle(d);
963 // if(i == tmp->GetFirstMother())
964 // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
965 // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
969 }//Herwig generated jets
976 //________________________________________________________________
977 void AliMCAnalysisUtils::Print(const Option_t * opt) const
979 //Print some relevant parameters set for the analysis
984 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
986 printf("Debug level = %d\n",fDebug);
987 printf("MC Generator = %s\n",fMCGenerator.Data());