PWG/CaloTrackCorrBase/AliMCAnalysisUtils.cxx

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