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