PWG4/PartCorrBase/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 /* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
  16
  17 //_________________________________________________________________________
  18 // Class for analysis utils for MC data
  19 // stored in stack or event header.
  20 // Contains:
  21 //  - method to check the origin of a given track/cluster
  22 //  - method to obtain the generated jets
  23 //
  24 //*-- Author: Gustavo Conesa (LNF-INFN)
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27
  28 // --- ROOT system ---
  29 #include <TMath.h>
  30 #include <TList.h>
  31 #include "TParticle.h"
  32 #include "TDatabasePDG.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(), fCurrentEvent(-1), fDebug(-1),
  46 fJetsList(new TList), fMCGenerator("PYTHIA")
  47 {
  48   //Ctor
  49 }
  50 /*
  51  //____________________________________________________________________________
  52  AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
  53  TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
  54  fJetsList(new TList), fMCGenerator(mcutils.fMCGenerator)
  55  {
  56  // cpy ctor
  57
  58  }
  59  */
  60
  61 //_________________________________________________________________________
  62 //AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
  63 //{
  64 //  // assignment operator
  65 //
  66 //  if(&mcutils == this) return *this;
  67 //  fCurrentEvent = mcutils.fCurrentEvent ;
  68 //  fDebug        = mcutils.fDebug;
  69 //  fJetsList     = mcutils.fJetsList;
  70 //  fMCGenerator  = mcutils.fMCGenerator;
  71 //
  72 //  return *this;
  73 //}
  74
  75 //____________________________________________________________________________
  76 AliMCAnalysisUtils::~AliMCAnalysisUtils()
  77 {
  78   // Remove all pointers.
  79
  80   if (fJetsList) {
  81     fJetsList->Clear();
  82     delete fJetsList ;
  83   }
  84 }
  85
  86
  87 //_________________________________________________________________________
  88 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, const Int_t nlabels, AliCaloTrackReader* reader, const Int_t input = 0) {
  89         //Play with the montecarlo particles if available
  90         Int_t tag = 0;
  91
  92   if(nlabels<=0) {
  93     printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
  94     return kMCBadLabel;
  95   }
  96
  97         //Select where the information is, ESD-galice stack or AOD mcparticles branch
  98         if(reader->ReadStack()){
  99                 tag = CheckOriginInStack(label, nlabels, reader->GetStack());
 100         }
 101         else if(reader->ReadAODMCParticles()){
 102                 tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(input));
 103         }
 104
 105         return tag ;
 106 }
 107
 108 //_________________________________________________________________________
 109 Int_t AliMCAnalysisUtils::CheckCommonAncestor(const Int_t index1, const Int_t index2, AliCaloTrackReader* reader,
 110                                               Int_t & ancPDG, Int_t & ancStatus, TLorentzVector & momentum) {
 111   //Check the first common ancestor of 2 clusters, given the most likely labels of the primaries generating such clusters.
 112   Int_t label1[100];
 113   Int_t label2[100];
 114   label1[0]= index1;
 115   label2[0]= index2;
 116   Int_t counter1 = 0;
 117   Int_t counter2 = 0;
 118
 119   if(label1[0]==label2[0]) {
 120     //printf("AliMCAnalysisUtils::CheckCommonAncestor() - Already the same label: %d\n",label1[0]);
 121     counter1=1;
 122     counter2=1;
 123   }
 124   else{
 125     if(reader->ReadAODMCParticles()){
 126       TClonesArray * mcparticles = reader->GetAODMCParticles(0);
 127
 128       Int_t label=label1[0];
 129       while(label > -1 && counter1 < 99){
 130         counter1++;
 131         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 132         if(mom){
 133          label  = mom->GetMother() ;
 134          label1[counter1]=label;
 135         }
 136         //printf("\t counter %d, label %d\n", counter1,label);
 137       }
 138       //printf("Org label2=%d,\n",label2[0]);
 139       label=label2[0];
 140       while(label > -1 && counter2 < 99){
 141         counter2++;
 142         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 143         if(mom){
 144           label  = mom->GetMother() ;
 145           label2[counter2]=label;
 146         }
 147         //printf("\t counter %d, label %d\n", counter2,label);
 148       }
 149     }//AOD MC
 150     else { //Kine stack from ESDs
 151       AliStack * stack = reader->GetStack();
 152       Int_t label=label1[0];
 153       while(label > -1 && counter1 < 99){
 154         counter1++;
 155         TParticle * mom = stack->Particle(label);
 156         if(mom){
 157           label  = mom->GetFirstMother() ;
 158           label1[counter1]=label;
 159         }
 160         //printf("\t counter %d, label %d\n", counter1,label);
 161       }
 162       //printf("Org label2=%d,\n",label2[0]);
 163       label=label2[0];
 164       while(label > -1 && counter2 < 99){
 165         counter2++;
 166         TParticle * mom = stack->Particle(label);
 167         if(mom){
 168           label  = mom->GetFirstMother() ;
 169           label2[counter2]=label;
 170         }
 171         //printf("\t counter %d, label %d\n", counter2,label);
 172       }
 173     }// Kine stack from ESDs
 174   }//First labels not the same
 175
 176   if(counter1==99 || counter2==99) printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
 177   //printf("CheckAncestor:\n");
 178   Int_t commonparents = 0;
 179   Int_t ancLabel = -1;
 180   //printf("counters %d %d \n",counter1, counter2);
 181   for (Int_t c1 = 0; c1 < counter1; c1++) {
 182     for (Int_t c2 = 0; c2 < counter2; c2++) {
 183       if(label1[c1]==label2[c2] && label1[c1]>-1) {
 184         ancLabel = label1[c1];
 185         commonparents++;
 186         if(reader->ReadAODMCParticles()){
 187           AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles(0)->At(label1[c1]);
 188           if (mom) {
 189             ancPDG    = mom->GetPdgCode();
 190             ancStatus = mom->GetStatus();
 191             momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
 192           }
 193         }
 194         else {
 195           TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
 196           if (mom) {
 197             ancPDG    = mom->GetPdgCode();
 198             ancStatus = mom->GetStatusCode();
 199             mom->Momentum(momentum);
 200           }
 201         }
 202         //First ancestor found, end the loops
 203         counter1=0;
 204         counter2=0;
 205       }//Ancestor found
 206     }//second cluster loop
 207   }//first cluster loop
 208
 209   return ancLabel;
 210 }
 211
 212 //_________________________________________________________________________
 213 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) {
 214         //Play with the montecarlo particles if available
 215         Int_t tag = 0;
 216
 217   if(label<0) {
 218     printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
 219     return kMCBadLabel;
 220   }
 221
 222         Int_t labels[]={label};
 223
 224         //Select where the information is, ESD-galice stack or AOD mcparticles branch
 225         if(reader->ReadStack()){
 226                 tag = CheckOriginInStack(labels, 1,reader->GetStack());
 227         }
 228         else if(reader->ReadAODMCParticles()){
 229                 tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(input));
 230         }
 231
 232         return tag ;
 233 }
 234
 235 //_________________________________________________________________________
 236 Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, const Int_t nlabels, AliStack* stack) {
 237   // Play with the MC stack if available. Tag particles depending on their origin.
 238   // Do same things as in CheckOriginInAOD but different input.
 239
 240   //generally speaking, label is the MC label of a reconstructed
 241   //entity (track, cluster, etc) for which we want to know something
 242   //about its heritage, but one can also use it directly with stack
 243   //particles not connected to reconstructed entities
 244
 245   if(!stack) {
 246     if (fDebug >=0)
 247       printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
 248     return -1;
 249   }
 250
 251   Int_t tag = 0;
 252   Int_t label=labels[0];//Most significant particle contributing to the cluster
 253
 254   if(label >= 0 && label < stack->GetNtrack()){
 255     //MC particle of interest is the "mom" of the entity
 256     TParticle * mom = stack->Particle(label);
 257     Int_t iMom     = label;
 258     Int_t mPdgSign = mom->GetPdgCode();
 259     Int_t mPdg     = TMath::Abs(mPdgSign);
 260     Int_t mStatus  = mom->GetStatusCode() ;
 261     Int_t iParent  = mom->GetFirstMother() ;
 262     if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 263
 264     //GrandParent of the entity
 265     TParticle * parent = NULL;
 266     Int_t pPdg = -1;
 267     Int_t pStatus =-1;
 268     if(iParent >= 0){
 269       parent = stack->Particle(iParent);
 270       if(parent){
 271         pPdg = TMath::Abs(parent->GetPdgCode());
 272         pStatus = parent->GetStatusCode();
 273       }
 274     }
 275     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
 276
 277     if(fDebug > 2 ) {
 278                   printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
 279                   printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 280                   printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 281     }
 282
 283     //Check if "mother" of entity is converted, if not, get the first non converted mother
 284     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
 285       SetTagBit(tag,kMCConversion);
 286       //Check if the mother is photon or electron with status not stable
 287       while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
 288         //Mother
 289         iMom     = mom->GetFirstMother();
 290         mom      = stack->Particle(iMom);
 291         mPdgSign = mom->GetPdgCode();
 292         mPdg     = TMath::Abs(mPdgSign);
 293         mStatus  = mom->GetStatusCode() ;
 294         iParent  = mom->GetFirstMother() ;
 295         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 296
 297         //GrandParent
 298         if(iParent >= 0){
 299           parent = stack->Particle(iParent);
 300           if(parent){
 301             pPdg = TMath::Abs(parent->GetPdgCode());
 302             pStatus = parent->GetStatusCode();
 303           }
 304         }
 305         else {// in case of gun/box simulations
 306           pPdg    = 0;
 307           pStatus = 0;
 308           break;
 309         }
 310       }//while
 311       if(fDebug > 2 ) {
 312         printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
 313         printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 314         printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 315       }
 316
 317     }//mother and parent are electron or photon and have status 0
 318     else if((mPdg == 22 || mPdg == 11) && mStatus == 0){
 319       //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
 320       if(pPdg == 2112 ||  pPdg == 211  ||  pPdg == 321 ||
 321          pPdg == 2212 ||  pPdg == 130  ||  pPdg == 13 ) {
 322         SetTagBit(tag,kMCConversion);
 323         iMom     = mom->GetFirstMother();
 324         mom      = stack->Particle(iMom);
 325         mPdgSign = mom->GetPdgCode();
 326         mPdg     = TMath::Abs(mPdgSign);
 327
 328         if(fDebug > 2 ) {
 329           printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
 330           printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 331         }
 332       }//hadron converted
 333
 334       //Comment for the next lines, we do not check the parent of the hadron for the moment.
 335       //iParent =  mom->GetFirstMother() ;
 336       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 337
 338       //GrandParent
 339       //if(iParent >= 0){
 340       //        parent = stack->Particle(iParent);
 341       //        pPdg = TMath::Abs(parent->GetPdgCode());
 342       //}
 343     }
 344     // conversion into electrons/photons checked
 345
 346     //first check for typical charged particles
 347     if     (mPdg     ==    13) SetTagBit(tag,kMCMuon);
 348     else if(mPdg     ==   211) SetTagBit(tag,kMCPion);
 349     else if(mPdg     ==   321) SetTagBit(tag,kMCKaon);
 350     else if(mPdgSign ==  2212) SetTagBit(tag,kMCProton);
 351     else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
 352     else if(mPdgSign ==  2112) SetTagBit(tag,kMCNeutron);
 353     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
 354
 355     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 356     else if(mPdg == 111)  {
 357       SetTagBit(tag,kMCPi0Decay);
 358       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
 359       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 360     }
 361     else if(mPdg == 221) {
 362       SetTagBit(tag,kMCEtaDecay);
 363       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
 364       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
 365     }
 366     //Photons
 367     else if(mPdg == 22){
 368       SetTagBit(tag,kMCPhoton);
 369       if(mStatus == 1){ //undecayed particle
 370         if(fMCGenerator == "PYTHIA"){
 371           if(iParent < 8 && iParent > 5) {//outgoing partons
 372             if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 373             else SetTagBit(tag,kMCFragmentation);
 374           }//Outgoing partons
 375           else  if(iParent <= 5) {
 376             SetTagBit(tag, kMCISR); //Initial state radiation
 377           }
 378           else if(pStatus == 11){//Decay
 379             if(pPdg == 111) {
 380               SetTagBit(tag,kMCPi0Decay);
 381               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 382               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 383             }
 384             else if (pPdg == 221) {
 385               SetTagBit(tag, kMCEtaDecay);
 386               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 387               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 388             }
 389             else SetTagBit(tag,kMCOtherDecay);
 390           }//Decay
 391           else {
 392             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",
 393                                    mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
 394             if(pPdg == 111) {
 395               SetTagBit(tag,kMCPi0Decay);
 396               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 397               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 398             }
 399             else if (pPdg == 221) {
 400               SetTagBit(tag, kMCEtaDecay);
 401               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 402               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 403             }
 404             else SetTagBit(tag,kMCOtherDecay);
 405           }
 406         }//PYTHIA
 407
 408         else if(fMCGenerator == "HERWIG"){
 409           if(pStatus < 197){//Not decay
 410             while(1){
 411               if(parent){
 412                 if(parent->GetFirstMother()<=5) break;
 413                 iParent = parent->GetFirstMother();
 414                 parent=stack->Particle(iParent);
 415                 pStatus= parent->GetStatusCode();
 416                 pPdg = TMath::Abs(parent->GetPdgCode());
 417               } else break;
 418             }//Look for the parton
 419
 420             if(iParent < 8 && iParent > 5) {
 421               if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 422               else SetTagBit(tag,kMCFragmentation);
 423             }
 424             else SetTagBit(tag,kMCISR);//Initial state radiation
 425           }//Not decay
 426           else{//Decay
 427             if(pPdg == 111) {
 428               SetTagBit(tag,kMCPi0Decay);
 429               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
 430               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 431             }
 432             else if (pPdg == 221) {
 433               SetTagBit(tag,kMCEtaDecay);
 434               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
 435               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
 436             }
 437             else SetTagBit(tag,kMCOtherDecay);
 438           }//Decay
 439         }//HERWIG
 440
 441         else SetTagBit(tag,kMCUnknown);
 442
 443       }//Status 1 : created by event generator
 444
 445       else if(mStatus == 0){ // geant
 446         if(pPdg == 111) {
 447           SetTagBit(tag,kMCPi0Decay);
 448           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
 449           CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 450         }
 451         else if (pPdg == 221) {
 452           SetTagBit(tag,kMCEtaDecay);
 453           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
 454           CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
 455         }
 456         else  SetTagBit(tag,kMCOtherDecay);
 457       }//status 0 : geant generated
 458
 459     }//Mother Photon
 460
 461     //Electron check.  Where did that electron come from?
 462     else if(mPdg == 11){ //electron
 463       if(pPdg == 11 && parent){
 464         Int_t iGrandma = parent->GetFirstMother();
 465         if(iGrandma >= 0) {
 466           TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
 467           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 468
 469           if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 470           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 471         }
 472       }
 473       SetTagBit(tag,kMCElectron);
 474       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
 475       if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
 476       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
 477       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
 478       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
 479         if(parent){
 480           Int_t iGrandma = parent->GetFirstMother();
 481           if(iGrandma >= 0) {
 482             TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
 483             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 484             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
 485             else SetTagBit(tag,kMCEFromC); //c-->e
 486           } else SetTagBit(tag,kMCEFromC); //c-->e
 487         }//parent
 488       } else {
 489         //if it is not from any of the above, where is it from?
 490         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 491         else SetTagBit(tag,kMCOtherDecay);
 492         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);
 493       }
 494     }//electron check
 495     //Cluster was made by something else
 496     else {
 497       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
 498       SetTagBit(tag,kMCUnknown);
 499     }
 500   }//Good label value
 501   else{// Bad label
 502
 503     if(label < 0 && (fDebug >= 0))
 504       printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***:  label %d \n", label);
 505     if(label >=  stack->GetNtrack() &&  (fDebug >= 0))
 506       printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 507     SetTagBit(tag,kMCUnknown);
 508   }//Bad label
 509
 510   return tag;
 511
 512 }
 513
 514
 515 //_________________________________________________________________________
 516 Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, const Int_t nlabels, TClonesArray *mcparticles) {
 517   // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
 518   // Do same things as in CheckOriginInStack but different input.
 519   if(!mcparticles) {
 520     if(fDebug >= 0)
 521       printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
 522     return -1;
 523   }
 524
 525   Int_t tag = 0;
 526   Int_t label=labels[0];//Most significant particle contributing to the cluster
 527
 528   Int_t nprimaries = mcparticles->GetEntriesFast();
 529   if(label >= 0 && label < nprimaries){
 530     //Mother
 531     AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 532     Int_t iMom     = label;
 533     Int_t mPdgSign = mom->GetPdgCode();
 534     Int_t mPdg     = TMath::Abs(mPdgSign);
 535     Int_t iParent  = mom->GetMother() ;
 536     if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 537
 538     //GrandParent
 539     AliAODMCParticle * parent = NULL ;
 540     Int_t pPdg = -1;
 541     if(iParent >= 0){
 542       parent = (AliAODMCParticle *) mcparticles->At(iParent);
 543       pPdg = TMath::Abs(parent->GetPdgCode());
 544     }
 545     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
 546
 547     if(fDebug > 2 ) {
 548                   printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
 549                   printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 550                   if(parent)
 551         printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 552     }
 553
 554     //Check if mother is converted, if not, get the first non converted mother
 555     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
 556       SetTagBit(tag,kMCConversion);
 557       //Check if the mother is photon or electron with status not stable
 558       while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
 559         //Mother
 560         iMom     = mom->GetMother();
 561         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
 562         mPdgSign = mom->GetPdgCode();
 563         mPdg     = TMath::Abs(mPdgSign);
 564         iParent  = mom->GetMother() ;
 565         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 566
 567         //GrandParent
 568         if(iParent >= 0 && parent){
 569           parent = (AliAODMCParticle *) mcparticles->At(iParent);
 570           pPdg = TMath::Abs(parent->GetPdgCode());
 571         }
 572         // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 573         // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 574
 575       }//while
 576
 577       if(fDebug > 2 ) {
 578         printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
 579         printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 580         if(parent)
 581           printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 582       }
 583
 584     }//mother and parent are electron or photon and have status 0 and parent is photon or electron
 585     else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){
 586       //Still a conversion but only one electron/photon generated. Just from hadrons
 587       if(pPdg == 2112 ||  pPdg == 211 ||  pPdg == 321 ||
 588          pPdg == 2212 ||  pPdg == 130 ||  pPdg == 13 ) {
 589         SetTagBit(tag,kMCConversion);
 590         iMom     = mom->GetMother();
 591         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
 592         mPdgSign = mom->GetPdgCode();
 593         mPdg     = TMath::Abs(mPdgSign);
 594
 595         if(fDebug > 2 ) {
 596           printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
 597           printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 598         }
 599       }//hadron converted
 600
 601       //Comment for next lines, we do not check the parent of the hadron for the moment.
 602       //iParent =  mom->GetMother() ;
 603       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 604
 605       //GrandParent
 606       //if(iParent >= 0){
 607       //        parent = (AliAODMCParticle *) mcparticles->At(iParent);
 608       //        pPdg = TMath::Abs(parent->GetPdgCode());
 609       //}
 610     }
 611
 612     // conversion into electrons/photons checked
 613
 614     //first check for typical charged particles
 615     if     (mPdg     ==    13) SetTagBit(tag,kMCMuon);
 616     else if(mPdg     ==   211) SetTagBit(tag,kMCPion);
 617     else if(mPdg     ==   321) SetTagBit(tag,kMCKaon);
 618     else if(mPdgSign ==  2212) SetTagBit(tag,kMCProton);
 619     else if(mPdgSign ==  2112) SetTagBit(tag,kMCNeutron);
 620     else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
 621     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
 622
 623     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 624     else if(mPdg == 111)  {
 625       SetTagBit(tag,kMCPi0Decay);
 626       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
 627       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 628     }
 629     else if(mPdg == 221)  {
 630       SetTagBit(tag,kMCEtaDecay);
 631       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
 632       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 633     }
 634     //Photons
 635     else if(mPdg == 22){
 636       SetTagBit(tag,kMCPhoton);
 637       if(mom->IsPhysicalPrimary()){ //undecayed particle
 638         if(iParent < 8 && iParent > 5) {//outgoing partons
 639           if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 640           else SetTagBit(tag,kMCFragmentation);
 641         }//Outgoing partons
 642         else if(iParent <= 5) {
 643           SetTagBit(tag, kMCISR); //Initial state radiation
 644         }
 645         else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
 646           if(pPdg == 111){
 647             SetTagBit(tag,kMCPi0Decay);
 648             if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
 649             CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 650           }
 651           else if (pPdg == 221) {
 652             SetTagBit(tag, kMCEtaDecay);
 653             if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
 654             CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 655           }
 656           else SetTagBit(tag,kMCOtherDecay);
 657         }//Decay
 658         else {
 659           if(parent)printf("AliMCAnalysisUtils::ChecOrigingInAOD() - what is it? Mother mPdg %d, is primary? %d, is physical %d \n    Parent  iParent %d, pPdg %d, is primary? %d, is physical? %d\n",
 660                  mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
 661           SetTagBit(tag,kMCOtherDecay);//Check
 662         }
 663       }// Pythia generated
 664       else if(!mom->IsPrimary()){       //Decays
 665         if(pPdg == 111){
 666           SetTagBit(tag,kMCPi0Decay);
 667           if(fDebug > 2 ) 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, const Int_t nlabels, const Int_t mesonIndex,
 741                                                     AliStack *stack, Int_t &tag) {
 742         //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
 743
 744         if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
 745                 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 746                           labels[0],stack->GetNtrack(), nlabels);
 747                 return;
 748         }
 749
 750         TParticle * meson = stack->Particle(mesonIndex);
 751         Int_t mesonPdg    = meson->GetPdgCode();
 752         if(mesonPdg!=111 && mesonPdg!=221){
 753                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 754                 return;
 755         }
 756
 757         if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
 758
 759         //Check if meson decayed into 2 daughters or if both were kept.
 760         if(meson->GetNDaughters() != 2){
 761                 if(fDebug > 2)
 762                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 763                 return;
 764         }
 765
 766         //Get the daughters
 767         Int_t iPhoton0 = meson->GetDaughter(0);
 768         Int_t iPhoton1 = meson->GetDaughter(1);
 769         TParticle *photon0 = stack->Particle(iPhoton0);
 770         TParticle *photon1 = stack->Particle(iPhoton1);
 771
 772         //Check if both daughters are photons
 773         if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
 774                 if(fDebug > 2)
 775                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 776                 return;
 777         }
 778
 779         if(fDebug > 2)
 780                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 781
 782         //Check if both photons contribute to the cluster
 783         Bool_t okPhoton0 = kFALSE;
 784         Bool_t okPhoton1 = kFALSE;
 785
 786         if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
 787
 788         for(Int_t i = 0; i < nlabels; i++){
 789                 if(fDebug > 3) printf("\t  at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
 790
 791                 //If we already found both, break the loop
 792                 if(okPhoton0 && okPhoton1) break;
 793
 794                 Int_t index =   labels[i];
 795                 if      (iPhoton0 == index) {
 796                         okPhoton0 = kTRUE;
 797                         continue;
 798                 }
 799                 else if (iPhoton1 == index) {
 800                         okPhoton1 = kTRUE;
 801                         continue;
 802                 }
 803
 804                 //Trace back the mother in case it was a conversion
 805                 TParticle * daught = stack->Particle(index);
 806                 Int_t tmpindex = daught->GetFirstMother();
 807                 if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
 808                 while(tmpindex>=0){
 809                         //MC particle of interest is the mother
 810                         if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
 811                         daught   = stack->Particle(tmpindex);
 812                         if      (iPhoton0 == tmpindex) {
 813                                 okPhoton0 = kTRUE;
 814                                 break;
 815                         }
 816                         else if (iPhoton1 == tmpindex) {
 817                                 okPhoton1 = kTRUE;
 818                                 break;
 819                         }
 820                         tmpindex = daught->GetFirstMother();
 821                 }//While to check if pi0/eta daughter was one of these contributors to the cluster
 822
 823                 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
 824                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
 825
 826         }//loop on list of labels
 827
 828         //If both photons contribute tag as the corresponding meson.
 829         if(okPhoton0 && okPhoton1){
 830                 if(fDebug > 2)
 831                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
 832
 833                 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 834                 else  SetTagBit(tag,kMCEta);
 835         }
 836
 837 }
 838
 839 //_________________________________________________________________________
 840 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
 841                                                     TClonesArray *mcparticles, Int_t &tag) {
 842         //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
 843
 844         if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1) {
 845                 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 846                           labels[0],mcparticles->GetEntriesFast(), nlabels);
 847                 return;
 848         }
 849
 850
 851   AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
 852         Int_t mesonPdg = meson->GetPdgCode();
 853         if(mesonPdg != 111 && mesonPdg != 221) {
 854                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 855                 return;
 856         }
 857
 858         if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
 859
 860
 861         //Get the daughters
 862         if(meson->GetNDaughters() != 2){
 863                 if(fDebug > 2)
 864                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 865                 return;
 866         }
 867         Int_t iPhoton0 = meson->GetDaughter(0);
 868         Int_t iPhoton1 = meson->GetDaughter(1);
 869         //if((iPhoton0 == -1) || (iPhoton1 == -1)){
 870         //      if(fDebug > 2)
 871         //              printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overalapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
 872         //      return;
 873         //}
 874         AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
 875         AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
 876
 877         //Check if both daughters are photons
 878         if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22){
 879                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 880                 return;
 881         }
 882
 883         if(fDebug > 2)
 884                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 885
 886         //Check if both photons contribute to the cluster
 887         Bool_t okPhoton0 = kFALSE;
 888         Bool_t okPhoton1 = kFALSE;
 889
 890         if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
 891
 892         for(Int_t i = 0; i < nlabels; i++){
 893                 if(fDebug > 3) printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
 894
 895                 //If we already found both, break the loop
 896                 if(okPhoton0 && okPhoton1) break;
 897
 898                 Int_t index =   labels[i];
 899                 if      (iPhoton0 == index) {
 900                         okPhoton0 = kTRUE;
 901                         continue;
 902                 }
 903                 else if (iPhoton1 == index) {
 904                         okPhoton1 = kTRUE;
 905                         continue;
 906                 }
 907
 908                 //Trace back the mother in case it was a conversion
 909                 AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
 910                 Int_t tmpindex = daught->GetMother();
 911                 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
 912
 913                 while(tmpindex>=0){
 914
 915                         //MC particle of interest is the mother
 916                         if(fDebug > 3) printf("\t parent index %d\n",tmpindex);
 917                         daught   = (AliAODMCParticle*) mcparticles->At(tmpindex);
 918                         //printf("tmpindex %d\n",tmpindex);
 919                         if      (iPhoton0 == tmpindex) {
 920                                 okPhoton0 = kTRUE;
 921                                 break;
 922                         }
 923                         else if (iPhoton1 == tmpindex) {
 924                                 okPhoton1 = kTRUE;
 925                                 break;
 926                         }
 927                         tmpindex = daught->GetMother();
 928                 }//While to check if pi0/eta daughter was one of these contributors to the cluster
 929
 930                 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
 931                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
 932
 933         }//loop on list of labels
 934
 935         //If both photons contribute tag as the corresponding meson.
 936         if(okPhoton0 && okPhoton1){
 937                 if(fDebug > 2)
 938                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
 939
 940                 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 941                 else  SetTagBit(tag,kMCEta);
 942         }
 943
 944 }
 945
 946 //_________________________________________________________________________
 947 TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * const reader){
 948   //Return list of jets (TParticles) and index of most likely parton that originated it.
 949   AliStack * stack = reader->GetStack();
 950   Int_t iEvent = reader->GetEventNumber();
 951   AliGenEventHeader * geh = reader->GetGenEventHeader();
 952   if(fCurrentEvent!=iEvent){
 953     fCurrentEvent = iEvent;
 954     fJetsList = new TList;
 955     Int_t nTriggerJets = 0;
 956     Float_t tmpjet[]={0,0,0,0};
 957
 958     //printf("Event %d %d\n",fCurrentEvent,iEvent);
 959     //Get outgoing partons
 960     if(stack->GetNtrack() < 8) return fJetsList;
 961     TParticle * parton1 =  stack->Particle(6);
 962     TParticle * parton2 =  stack->Particle(7);
 963     if(fDebug > 2){
 964       printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 965              parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
 966       printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 967              parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
 968                 }
 969     //          //Trace the jet from the mother parton
 970     //          Float_t pt  = 0;
 971     //          Float_t pt1 = 0;
 972     //          Float_t pt2 = 0;
 973     //          Float_t e   = 0;
 974     //          Float_t e1  = 0;
 975     //          Float_t e2  = 0;
 976     //          TParticle * tmptmp = new TParticle;
 977     //          for(Int_t i = 0; i< stack->GetNprimary(); i++){
 978     //                  tmptmp = stack->Particle(i);
 979
 980     //                  if(tmptmp->GetStatusCode() == 1){
 981     //                          pt = tmptmp->Pt();
 982     //                          e =  tmptmp->Energy();
 983     //                          Int_t imom = tmptmp->GetFirstMother();
 984     //                          Int_t imom1 = 0;
 985     //                          //printf("1st imom %d\n",imom);
 986     //                          while(imom > 5){
 987     //                                  imom1=imom;
 988     //                                  tmptmp = stack->Particle(imom);
 989     //                                  imom = tmptmp->GetFirstMother();
 990     //                                  //printf("imom %d       \n",imom);
 991     //                          }
 992     //                          //printf("Last imom %d %d\n",imom1, imom);
 993     //                          if(imom1 == 6) {
 994     //                                  pt1+=pt;
 995     //                                  e1+=e;
 996     //                          }
 997     //                          else if (imom1 == 7){
 998     //                                  pt2+=pt;
 999     //                                  e2+=e;                                  }
1000     //                  }// status 1
1001
1002     //          }// for
1003
1004     //          printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
1005
1006                 //Get the jet, different way for different generator
1007                 //PYTHIA
1008     if(fMCGenerator == "PYTHIA"){
1009       TParticle * jet =  0x0;
1010       AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
1011       nTriggerJets =  pygeh->NTriggerJets();
1012       if(fDebug > 1)
1013         printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
1014
1015       Int_t iparton = -1;
1016       for(Int_t i = 0; i< nTriggerJets; i++){
1017         iparton=-1;
1018         pygeh->TriggerJet(i, tmpjet);
1019         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1020         //Assign an outgoing parton as mother
1021         Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1022         Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1023         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1024         else  jet->SetFirstMother(6);
1025         //jet->Print();
1026         if(fDebug > 1)
1027           printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1028                  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1029         fJetsList->Add(jet);
1030       }
1031     }//Pythia triggered jets
1032     //HERWIG
1033     else if (fMCGenerator=="HERWIG"){
1034       Int_t pdg = -1;
1035       //Check parton 1
1036       TParticle * tmp = parton1;
1037       if(parton1->GetPdgCode()!=22){
1038         while(pdg != 94){
1039           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1040           tmp = stack->Particle(tmp->GetFirstDaughter());
1041           pdg = tmp->GetPdgCode();
1042         }//while
1043
1044         //Add found jet to list
1045         TParticle *jet1 = new TParticle(*tmp);
1046         jet1->SetFirstMother(6);
1047         fJetsList->Add(jet1);
1048         //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1049         //tmp = stack->Particle(tmp->GetFirstDaughter());
1050         //tmp->Print();
1051         //jet1->Print();
1052         if(fDebug > 1)
1053           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1054                  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
1055       }//not photon
1056
1057       //Check parton 2
1058       pdg = -1;
1059       tmp = parton2;
1060       Int_t i = -1;
1061       if(parton2->GetPdgCode()!=22){
1062         while(pdg != 94){
1063           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1064           i = tmp->GetFirstDaughter();
1065           tmp = stack->Particle(tmp->GetFirstDaughter());
1066           pdg = tmp->GetPdgCode();
1067         }//while
1068         //Add found jet to list
1069         TParticle *jet2 = new TParticle(*tmp);
1070         jet2->SetFirstMother(7);
1071         fJetsList->Add(jet2);
1072         //jet2->Print();
1073         if(fDebug > 1)
1074           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1075                  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1076         //Int_t first =  tmp->GetFirstDaughter();
1077         //Int_t last  =  tmp->GetLastDaughter();
1078         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1079                                 //      for(Int_t d = first ; d < last+1; d++){
1080         //                                              tmp = stack->Particle(d);
1081         //                                              if(i == tmp->GetFirstMother())
1082         //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1083         //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1084         //                         }
1085         //tmp->Print();
1086       }//not photon
1087     }//Herwig generated jets
1088   }
1089
1090   return fJetsList;
1091 }
1092
1093
1094 //________________________________________________________________
1095 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1096 {
1097   //Print some relevant parameters set for the analysis
1098
1099   if(! opt)
1100     return;
1101
1102   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1103
1104   printf("Debug level    = %d\n",fDebug);
1105   printf("MC Generator   = %s\n",fMCGenerator.Data());
1106   printf(" \n");
1107
1108 }
1109
1110