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