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