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 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1121 {
1122   //Print some relevant parameters set for the analysis
1123
1124   if(! opt)
1125     return;
1126
1127   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1128
1129   printf("Debug level    = %d\n",fDebug);
1130   printf("MC Generator   = %s\n",fMCGenerator.Data());
1131   printf(" \n");
1132
1133 }
1134
1135 //________________________________________________________
1136 void AliMCAnalysisUtils::PrintMCTag(const Int_t tag) const
1137 {
1138   // print the assigned origins to this particle
1139
1140   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",
1141          tag,
1142          CheckTagBit(tag,kMCPhoton),
1143          CheckTagBit(tag,kMCConversion),
1144          CheckTagBit(tag,kMCPrompt),
1145          CheckTagBit(tag,kMCFragmentation),
1146          CheckTagBit(tag,kMCISR),
1147          CheckTagBit(tag,kMCPi0Decay),
1148          CheckTagBit(tag,kMCEtaDecay),
1149          CheckTagBit(tag,kMCOtherDecay),
1150          CheckTagBit(tag,kMCPi0),
1151          CheckTagBit(tag,kMCEta),
1152          CheckTagBit(tag,kMCElectron),
1153          CheckTagBit(tag,kMCMuon),
1154          CheckTagBit(tag,kMCPion),
1155          CheckTagBit(tag,kMCProton),
1156          CheckTagBit(tag,kMCAntiNeutron),
1157          CheckTagBit(tag,kMCKaon),
1158          CheckTagBit(tag,kMCAntiProton),
1159          CheckTagBit(tag,kMCAntiNeutron),
1160          CheckTagBit(tag,kMCOther),
1161          CheckTagBit(tag,kMCUnknown),
1162          CheckTagBit(tag,kMCBadLabel)
1163          );
1164
1165 }
1166
1167