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(Int_t index1, 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       {
  96         counter1++;
  97         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
  98         if(mom)
  99         {
 100           label  = mom->GetMother() ;
 101           label1[counter1]=label;
 102         }
 103         //printf("\t counter %d, label %d\n", counter1,label);
 104       }
 105
 106       //printf("Org label2=%d,\n",label2[0]);
 107       label=label2[0];
 108       if(label < 0) return -1;
 109
 110       while(label > -1 && counter2 < 99)
 111       {
 112         counter2++;
 113         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 114         if(mom)
 115         {
 116           label  = mom->GetMother() ;
 117           label2[counter2]=label;
 118         }
 119         //printf("\t counter %d, label %d\n", counter2,label);
 120       }
 121     }//AOD MC
 122     else
 123     { //Kine stack from ESDs
 124       AliStack * stack = reader->GetStack();
 125
 126       Int_t label=label1[0];
 127       while(label > -1 && counter1 < 99)
 128       {
 129         counter1++;
 130         TParticle * mom = stack->Particle(label);
 131         if(mom)
 132         {
 133           label  = mom->GetFirstMother() ;
 134           label1[counter1]=label;
 135         }
 136         //printf("\t counter %d, label %d\n", counter1,label);
 137       }
 138
 139       //printf("Org label2=%d,\n",label2[0]);
 140
 141       label=label2[0];
 142       while(label > -1 && counter2 < 99)
 143       {
 144         counter2++;
 145         TParticle * mom = stack->Particle(label);
 146         if(mom)
 147         {
 148           label  = mom->GetFirstMother() ;
 149           label2[counter2]=label;
 150         }
 151         //printf("\t counter %d, label %d\n", counter2,label);
 152       }
 153     }// Kine stack from ESDs
 154   }//First labels not the same
 155
 156   if((counter1==99 || counter2==99) && fDebug >=0)
 157     printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
 158   //printf("CheckAncestor:\n");
 159
 160   Int_t commonparents = 0;
 161   Int_t ancLabel = -1;
 162   //printf("counters %d %d \n",counter1, counter2);
 163   for (Int_t c1 = 0; c1 < counter1; c1++)
 164   {
 165     for (Int_t c2 = 0; c2 < counter2; c2++)
 166     {
 167       if(label1[c1]==label2[c2] && label1[c1]>-1)
 168       {
 169         ancLabel = label1[c1];
 170         commonparents++;
 171
 172         if(reader->ReadAODMCParticles())
 173         {
 174           AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]);
 175
 176           if (mom)
 177           {
 178             ancPDG    = mom->GetPdgCode();
 179             ancStatus = mom->GetStatus();
 180             momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
 181             prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
 182           }
 183         }
 184         else
 185         {
 186           TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
 187
 188           if (mom)
 189           {
 190             ancPDG    = mom->GetPdgCode();
 191             ancStatus = mom->GetStatusCode();
 192             mom->Momentum(momentum);
 193             prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz());
 194           }
 195         }
 196
 197         //First ancestor found, end the loops
 198         counter1=0;
 199         counter2=0;
 200       }//Ancestor found
 201     }//second cluster loop
 202   }//first cluster loop
 203
 204   if(ancLabel < 0)
 205   {
 206     ancPDG    = -10000;
 207     ancStatus = -10000;
 208     momentum.SetXYZT(0,0,0,0);
 209     prodVertex.SetXYZ(-10,-10,-10);
 210   }
 211
 212   return ancLabel;
 213 }
 214
 215 //________________________________________________________________________________________
 216 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, Int_t nlabels,
 217                                       const AliCaloTrackReader* reader, TString calorimeter)
 218 {
 219   //Play with the montecarlo particles if available
 220   Int_t tag = 0;
 221
 222   if(nlabels<=0) {
 223     printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
 224     return kMCBadLabel;
 225   }
 226
 227   TObjArray* arrayCluster = 0;
 228   if      ( calorimeter == "EMCAL" ) arrayCluster = reader->GetEMCALClusters();
 229   else if ( calorimeter ==  "PHOS" ) arrayCluster= reader->GetPHOSClusters();
 230
 231   //Select where the information is, ESD-galice stack or AOD mcparticles branch
 232   if(reader->ReadStack()){
 233     tag = CheckOriginInStack(label, nlabels, reader->GetStack(), arrayCluster);
 234   }
 235   else if(reader->ReadAODMCParticles()){
 236     tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(),arrayCluster);
 237   }
 238
 239   return tag ;
 240 }
 241
 242 //____________________________________________________________________________________________________
 243 Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, const AliCaloTrackReader* reader, TString calorimeter)
 244 {
 245   //Play with the montecarlo particles if available
 246   Int_t tag = 0;
 247
 248   if(label<0) {
 249     printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
 250     return kMCBadLabel;
 251   }
 252
 253   TObjArray* arrayCluster = 0;
 254   if      ( calorimeter == "EMCAL" ) arrayCluster = reader->GetEMCALClusters();
 255   else if ( calorimeter ==  "PHOS" ) arrayCluster= reader->GetPHOSClusters();
 256
 257   Int_t labels[]={label};
 258
 259   //Select where the information is, ESD-galice stack or AOD mcparticles branch
 260   if(reader->ReadStack()){
 261     tag = CheckOriginInStack(labels, 1,reader->GetStack(),arrayCluster);
 262   }
 263   else if(reader->ReadAODMCParticles()){
 264     tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(),arrayCluster);
 265   }
 266
 267   return tag ;
 268 }
 269
 270 //__________________________________________________________________________________________
 271 Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, Int_t nlabels,
 272                                              AliStack* stack, const TObjArray* arrayCluster)
 273 {
 274   // Play with the MC stack if available. Tag particles depending on their origin.
 275   // Do same things as in CheckOriginInAOD but different input.
 276
 277   //generally speaking, label is the MC label of a reconstructed
 278   //entity (track, cluster, etc) for which we want to know something
 279   //about its heritage, but one can also use it directly with stack
 280   //particles not connected to reconstructed entities
 281
 282   if(!stack)
 283   {
 284     if (fDebug >=0)
 285       printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
 286     return -1;
 287   }
 288
 289   Int_t tag = 0;
 290   Int_t label=labels[0];//Most significant particle contributing to the cluster
 291
 292   if(label >= 0 && label < stack->GetNtrack())
 293   {
 294     //MC particle of interest is the "mom" of the entity
 295     TParticle * mom = stack->Particle(label);
 296     Int_t iMom     = label;
 297     Int_t mPdgSign = mom->GetPdgCode();
 298     Int_t mPdg     = TMath::Abs(mPdgSign);
 299     Int_t mStatus  = mom->GetStatusCode() ;
 300     Int_t iParent  = mom->GetFirstMother() ;
 301
 302     if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 303
 304     //GrandParent of the entity
 305     TParticle * parent = NULL;
 306     Int_t pPdg = -1;
 307     Int_t pStatus =-1;
 308     if(iParent >= 0)
 309     {
 310       parent = stack->Particle(iParent);
 311
 312       if(parent)
 313       {
 314         pPdg = TMath::Abs(parent->GetPdgCode());
 315         pStatus = parent->GetStatusCode();
 316       }
 317     }
 318     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
 319
 320     if(fDebug > 2 )
 321     {
 322       printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
 323       printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 324       printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 325     }
 326
 327     //Check if "mother" of entity is converted, if not, get the first non converted mother
 328     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0)
 329     {
 330       SetTagBit(tag,kMCConversion);
 331
 332       //Check if the mother is photon or electron with status not stable
 333       while ((pPdg == 22 || pPdg == 11) && mStatus != 1)
 334       {
 335         //Mother
 336         iMom     = mom->GetFirstMother();
 337         mom      = stack->Particle(iMom);
 338         mPdgSign = mom->GetPdgCode();
 339         mPdg     = TMath::Abs(mPdgSign);
 340         mStatus  = mom->GetStatusCode() ;
 341         iParent  = mom->GetFirstMother() ;
 342
 343         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 344
 345         //GrandParent
 346         if(iParent >= 0)
 347         {
 348           parent = stack->Particle(iParent);
 349           if(parent)
 350           {
 351             pPdg = TMath::Abs(parent->GetPdgCode());
 352             pStatus = parent->GetStatusCode();
 353           }
 354         }
 355         else {// in case of gun/box simulations
 356           pPdg    = 0;
 357           pStatus = 0;
 358           break;
 359         }
 360       }//while
 361
 362       if(fDebug > 2 )
 363       {
 364         printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
 365         printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 366         printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 367       }
 368
 369     }//mother and parent are electron or photon and have status 0
 370     else if((mPdg == 22 || mPdg == 11) && mStatus == 0)
 371     {
 372       //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
 373       if(pPdg == 2112 ||  pPdg == 211  ||  pPdg == 321 ||
 374          pPdg == 2212 ||  pPdg == 130  ||  pPdg == 13 )
 375       {
 376         SetTagBit(tag,kMCConversion);
 377         iMom     = mom->GetFirstMother();
 378         mom      = stack->Particle(iMom);
 379         mPdgSign = mom->GetPdgCode();
 380         mPdg     = TMath::Abs(mPdgSign);
 381
 382         if(fDebug > 2 )
 383         {
 384           printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
 385           printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 386         }
 387       }//hadron converted
 388
 389       //Comment for the next lines, we do not check the parent of the hadron for the moment.
 390       //iParent =  mom->GetFirstMother() ;
 391       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 392
 393       //GrandParent
 394       //if(iParent >= 0){
 395       //        parent = stack->Particle(iParent);
 396       //        pPdg = TMath::Abs(parent->GetPdgCode());
 397       //}
 398     }
 399     // conversion into electrons/photons checked
 400
 401     //first check for typical charged particles
 402     if     (mPdg     ==    13) SetTagBit(tag,kMCMuon);
 403     else if(mPdg     ==   211) SetTagBit(tag,kMCPion);
 404     else if(mPdg     ==   321) SetTagBit(tag,kMCKaon);
 405     else if(mPdgSign ==  2212) SetTagBit(tag,kMCProton);
 406     else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
 407     else if(mPdgSign ==  2112) SetTagBit(tag,kMCNeutron);
 408     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
 409
 410     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 411     else if(mPdg == 111)
 412     {
 413
 414       SetTagBit(tag,kMCPi0Decay);
 415
 416       if(fDebug > 2 )
 417         printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
 418
 419       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 420
 421     }
 422     else if(mPdg == 221)
 423     {
 424       SetTagBit(tag,kMCEtaDecay);
 425
 426       if(fDebug > 2 )
 427         printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
 428
 429       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
 430     }
 431     //Photons
 432     else if(mPdg == 22)
 433     {
 434       SetTagBit(tag,kMCPhoton);
 435
 436       if(pPdg == 111)
 437       {
 438         SetTagBit(tag,kMCPi0Decay);
 439
 440         if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 441
 442         CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 443         // In case it did not merge, check if the decay companion is lost
 444         if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo))
 445           CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag);
 446
 447         //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
 448       }
 449       else if (pPdg == 221)
 450       {
 451         SetTagBit(tag, kMCEtaDecay);
 452
 453         if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 454
 455         CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 456         // In case it did not merge, check if the decay companion is lost
 457         if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo))
 458           CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag);
 459       }
 460       else if(mStatus == 1)
 461       { //undecayed particle
 462         if(fMCGenerator == "PYTHIA")
 463         {
 464           if(iParent < 8 && iParent > 5)
 465           {//outgoing partons
 466             if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 467             else           SetTagBit(tag,kMCFragmentation);
 468           }//Outgoing partons
 469           else  if(iParent <= 5)
 470           {
 471             SetTagBit(tag, kMCISR); //Initial state radiation
 472           }
 473           else  SetTagBit(tag,kMCUnknown);
 474          }//PYTHIA
 475
 476         else if(fMCGenerator == "HERWIG")
 477         {
 478           if(pStatus < 197)
 479           {//Not decay
 480             while(1)
 481             {
 482               if(parent)
 483               {
 484                 if(parent->GetFirstMother()<=5) break;
 485                 iParent = parent->GetFirstMother();
 486                 parent=stack->Particle(iParent);
 487                 pStatus= parent->GetStatusCode();
 488                 pPdg = TMath::Abs(parent->GetPdgCode());
 489               } else break;
 490             }//Look for the parton
 491
 492             if(iParent < 8 && iParent > 5)
 493             {
 494               if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 495               else           SetTagBit(tag,kMCFragmentation);
 496             }
 497             else SetTagBit(tag,kMCISR);//Initial state radiation
 498           }//Not decay
 499           else  SetTagBit(tag,kMCUnknown);
 500         }//HERWIG
 501       }
 502       else  SetTagBit(tag,kMCOtherDecay);
 503
 504     }//Mother Photon
 505
 506     //Electron check.  Where did that electron come from?
 507     else if(mPdg == 11){ //electron
 508       if(pPdg == 11 && parent)
 509       {
 510         Int_t iGrandma = parent->GetFirstMother();
 511         if(iGrandma >= 0)
 512         {
 513           TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
 514           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 515
 516           if      (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 517           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 518         }
 519       }
 520
 521       SetTagBit(tag,kMCElectron);
 522
 523       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
 524
 525       if      (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag, kMCDecayDalitz) ; } //Pi0 Dalitz decay
 526       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag, kMCDecayDalitz) ; } //Eta Dalitz decay
 527       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
 528       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
 529         if(parent)
 530         {
 531           Int_t iGrandma = parent->GetFirstMother();
 532           if(iGrandma >= 0)
 533           {
 534             TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
 535             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 536             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
 537             else SetTagBit(tag,kMCEFromC); //c-->e
 538           }
 539           else SetTagBit(tag,kMCEFromC); //c-->e
 540         }//parent
 541       }
 542       else
 543       {
 544         //if it is not from any of the above, where is it from?
 545         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 546
 547         else SetTagBit(tag,kMCOtherDecay);
 548
 549         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);
 550       }
 551     }//electron check
 552     //Cluster was made by something else
 553     else
 554     {
 555       if(fDebug > 0)
 556         printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",
 557                mom->GetName(),mPdg,pPdg);
 558
 559       SetTagBit(tag,kMCUnknown);
 560     }
 561   }//Good label value
 562   else
 563   {// Bad label
 564     if(label < 0 && (fDebug >= 0))
 565       printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***:  label %d \n", label);
 566
 567     if(label >=  stack->GetNtrack() &&  (fDebug >= 0))
 568       printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 569
 570     SetTagBit(tag,kMCUnknown);
 571   }//Bad label
 572
 573   return tag;
 574
 575 }
 576
 577
 578 //________________________________________________________________________________________________________
 579 Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, Int_t nlabels,
 580                                            const TClonesArray *mcparticles, const TObjArray* arrayCluster)
 581 {
 582   // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
 583   // Do same things as in CheckOriginInStack but different input.
 584   if(!mcparticles)
 585   {
 586     if(fDebug >= 0)
 587       printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
 588     return -1;
 589   }
 590
 591   Int_t tag = 0;
 592   Int_t label=labels[0];//Most significant particle contributing to the cluster
 593
 594   Int_t nprimaries = mcparticles->GetEntriesFast();
 595   if(label >= 0 && label < nprimaries)
 596   {
 597     //Mother
 598     AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 599     Int_t iMom     = label;
 600     Int_t mPdgSign = mom->GetPdgCode();
 601     Int_t mPdg     = TMath::Abs(mPdgSign);
 602     Int_t iParent  = mom->GetMother() ;
 603
 604     if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 605
 606     //GrandParent
 607     AliAODMCParticle * parent = NULL ;
 608     Int_t pPdg = -1;
 609     if(iParent >= 0)
 610     {
 611       parent = (AliAODMCParticle *) mcparticles->At(iParent);
 612       pPdg = TMath::Abs(parent->GetPdgCode());
 613     }
 614     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
 615
 616     if(fDebug > 2 )
 617     {
 618       printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
 619
 620       printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
 621              iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 622
 623       if(parent)
 624         printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
 625                iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 626     }
 627
 628     //Check if mother is converted, if not, get the first non converted mother
 629     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary())
 630     {
 631       SetTagBit(tag,kMCConversion);
 632
 633       //Check if the mother is photon or electron with status not stable
 634       while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
 635       {
 636         //Mother
 637         iMom     = mom->GetMother();
 638         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
 639         mPdgSign = mom->GetPdgCode();
 640         mPdg     = TMath::Abs(mPdgSign);
 641         iParent  = mom->GetMother() ;
 642         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 643
 644         //GrandParent
 645         if(iParent >= 0 && parent)
 646         {
 647           parent = (AliAODMCParticle *) mcparticles->At(iParent);
 648           pPdg = TMath::Abs(parent->GetPdgCode());
 649         }
 650         // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 651         // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 652
 653       }//while
 654
 655       if(fDebug > 2 )
 656       {
 657         printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
 658
 659         printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
 660                iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 661
 662         if(parent)
 663           printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
 664                  iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 665       }
 666
 667     }//mother and parent are electron or photon and have status 0 and parent is photon or electron
 668     else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary())
 669     {
 670       //Still a conversion but only one electron/photon generated. Just from hadrons
 671       if(pPdg == 2112 ||  pPdg == 211 ||  pPdg == 321 ||
 672          pPdg == 2212 ||  pPdg == 130 ||  pPdg == 13 )
 673       {
 674         SetTagBit(tag,kMCConversion);
 675         iMom     = mom->GetMother();
 676         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
 677         mPdgSign = mom->GetPdgCode();
 678         mPdg     = TMath::Abs(mPdgSign);
 679
 680         if(fDebug > 2 )
 681         {
 682           printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
 683           printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 684         }
 685       }//hadron converted
 686
 687       //Comment for next lines, we do not check the parent of the hadron for the moment.
 688       //iParent =  mom->GetMother() ;
 689       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 690
 691       //GrandParent
 692       //if(iParent >= 0){
 693       //        parent = (AliAODMCParticle *) mcparticles->At(iParent);
 694       //        pPdg = TMath::Abs(parent->GetPdgCode());
 695       //}
 696     }
 697
 698     //printf("Final mother mPDG %d\n",mPdg);
 699
 700     // conversion into electrons/photons checked
 701
 702     //first check for typical charged particles
 703     if     (mPdg     ==    13) SetTagBit(tag,kMCMuon);
 704     else if(mPdg     ==   211) SetTagBit(tag,kMCPion);
 705     else if(mPdg     ==   321) SetTagBit(tag,kMCKaon);
 706     else if(mPdgSign ==  2212) SetTagBit(tag,kMCProton);
 707     else if(mPdgSign ==  2112) SetTagBit(tag,kMCNeutron);
 708     else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
 709     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
 710
 711     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 712     else if(mPdg == 111)
 713     {
 714       SetTagBit(tag,kMCPi0Decay);
 715
 716       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
 717
 718       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 719     }
 720     else if(mPdg == 221)
 721     {
 722       SetTagBit(tag,kMCEtaDecay);
 723
 724       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
 725
 726       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 727     }
 728     //Photons
 729     else if(mPdg == 22)
 730     {
 731       SetTagBit(tag,kMCPhoton);
 732
 733       if(pPdg == 111)
 734       {
 735         SetTagBit(tag,kMCPi0Decay);
 736
 737         if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
 738
 739         CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 740         // In case it did not merge, check if the decay companion is lost
 741         if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo) && !CheckTagBit(tag,kMCDecayPairLost))
 742           CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag);
 743
 744         //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
 745       }
 746       else if (pPdg == 221)
 747       {
 748         SetTagBit(tag, kMCEtaDecay);
 749
 750         if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
 751
 752         CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 753         // In case it did not merge, check if the decay companion is lost
 754         if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo))
 755           CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag);
 756       }
 757       else if(mom->IsPhysicalPrimary() && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) //undecayed particle
 758       {
 759         if(iParent < 8 && iParent > 5 )
 760         {//outgoing partons
 761           if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 762           else SetTagBit(tag,kMCFragmentation);
 763         }//Outgoing partons
 764         else if(iParent <= 5 && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG"))
 765         {
 766           SetTagBit(tag, kMCISR); //Initial state radiation
 767         }
 768         else SetTagBit(tag,kMCUnknown);
 769       }//Physical primary
 770       else SetTagBit(tag,kMCOtherDecay);
 771
 772     }//Mother Photon
 773
 774     //Electron check.  Where did that electron come from?
 775     else if(mPdg == 11)
 776     { //electron
 777       if(pPdg == 11 && parent)
 778       {
 779         Int_t iGrandma = parent->GetMother();
 780         if(iGrandma >= 0)
 781         {
 782           AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
 783           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 784
 785           if      (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 786           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 787         }
 788       }
 789
 790       SetTagBit(tag,kMCElectron);
 791
 792       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
 793
 794       if      (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay
 795       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay
 796       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
 797       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
 798       { //c-hadron decay check
 799         if(parent)
 800         {
 801           Int_t iGrandma = parent->GetMother();
 802           if(iGrandma >= 0)
 803           {
 804             AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
 805             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 806             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
 807             else SetTagBit(tag,kMCEFromC); //c-hadron decay
 808           }
 809           else SetTagBit(tag,kMCEFromC); //c-hadron decay
 810         }//parent
 811       } else
 812       { //prompt or other decay
 813         TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
 814         TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
 815
 816         if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",
 817                               foo->GetName(), pPdg,foo1->GetName(),mPdg);
 818
 819         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 820         else             SetTagBit(tag,kMCOtherDecay);
 821       }
 822     }//electron check
 823     //cluster was made by something else
 824     else
 825     {
 826       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",
 827                             mPdg,pPdg);
 828       SetTagBit(tag,kMCUnknown);
 829     }
 830   }//Good label value
 831   else
 832   {//Bad label
 833     if(label < 0 && (fDebug >= 0) )
 834       printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***:  label %d \n", label);
 835
 836     if(label >=  mcparticles->GetEntriesFast() &&  (fDebug >= 0) )
 837       printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***:  label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
 838
 839     SetTagBit(tag,kMCUnknown);
 840
 841   }//Bad label
 842
 843   return tag;
 844
 845 }
 846
 847 //_________________________________________________________________________________________
 848 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,    Int_t nlabels,
 849                                                     Int_t mesonIndex, AliStack *stack,
 850                                                     Int_t &tag)
 851 {
 852   //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
 853
 854   if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
 855     if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 856                           labels[0],stack->GetNtrack(), nlabels);
 857     return;
 858   }
 859
 860   TParticle * meson = stack->Particle(mesonIndex);
 861   Int_t mesonPdg    = meson->GetPdgCode();
 862   if(mesonPdg!=111 && mesonPdg!=221)
 863   {
 864     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 865     return;
 866   }
 867
 868   if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s, label %d\n",meson->GetName(), mesonIndex);
 869
 870   //Check if meson decayed into 2 daughters or if both were kept.
 871   if(meson->GetNDaughters() != 2)
 872   {
 873     if(fDebug > 2)
 874       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 875     return;
 876   }
 877
 878   //Get the daughters
 879   Int_t iPhoton0 = meson->GetDaughter(0);
 880   Int_t iPhoton1 = meson->GetDaughter(1);
 881   TParticle *photon0 = stack->Particle(iPhoton0);
 882   TParticle *photon1 = stack->Particle(iPhoton1);
 883
 884   //Check if both daughters are photons
 885   if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22)
 886   {
 887     if(fDebug > 2)
 888       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 889     return;
 890   }
 891
 892   if(fDebug > 2)
 893     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 894
 895   //Check if both photons contribute to the cluster
 896   Bool_t okPhoton0 = kFALSE;
 897   Bool_t okPhoton1 = kFALSE;
 898
 899   if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Labels loop:\n");
 900
 901   Bool_t conversion = kFALSE;
 902
 903   for(Int_t i = 0; i < nlabels; i++)
 904   {
 905     if(fDebug > 3) printf("\t  at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
 906
 907     //If we already found both, break the loop
 908     if(okPhoton0 && okPhoton1) break;
 909
 910     Int_t index = labels[i];
 911     if      (iPhoton0 == index)
 912     {
 913       okPhoton0 = kTRUE;
 914       continue;
 915     }
 916     else if (iPhoton1 == index)
 917     {
 918       okPhoton1 = kTRUE;
 919       continue;
 920     }
 921
 922     //Trace back the mother in case it was a conversion
 923
 924     if(index >= stack->GetNtrack())
 925     {
 926       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",
 927              index,stack->GetNtrack());
 928       continue;
 929     }
 930
 931     TParticle * daught = stack->Particle(index);
 932     Int_t tmpindex = daught->GetFirstMother();
 933     if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
 934     while(tmpindex>=0)
 935     {
 936       //MC particle of interest is the mother
 937       if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
 938       daught   = stack->Particle(tmpindex);
 939       if      (iPhoton0 == tmpindex)
 940       {
 941         conversion = kTRUE;
 942         okPhoton0  = kTRUE;
 943         break;
 944       }
 945       else if (iPhoton1 == tmpindex)
 946       {
 947         conversion = kTRUE;
 948         okPhoton1  = kTRUE;
 949         break;
 950       }
 951
 952       tmpindex = daught->GetFirstMother();
 953
 954     }//While to check if pi0/eta daughter was one of these contributors to the cluster
 955
 956     if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
 957       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Something happens, first label should be from a photon decay!\n");
 958
 959   }//loop on list of labels
 960
 961   //If both photons contribute tag as the corresponding meson.
 962   if(okPhoton0 && okPhoton1)
 963   {
 964     if(fDebug > 2)
 965       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s OVERLAPPED DECAY \n", meson->GetName());
 966
 967     if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ;
 968
 969     if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 970     else                SetTagBit(tag,kMCEta);
 971   }
 972
 973 }
 974
 975 //________________________________________________________________________________________________________
 976 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex,
 977                                                     const TClonesArray *mcparticles, Int_t & tag   )
 978 {
 979   //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
 980
 981   if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1)
 982   {
 983     if(fDebug > 2)
 984       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 985              labels[0],mcparticles->GetEntriesFast(), nlabels);
 986     return;
 987   }
 988
 989   AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
 990   Int_t mesonPdg = meson->GetPdgCode();
 991   if(mesonPdg != 111 && mesonPdg != 221)
 992   {
 993     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 994     return;
 995   }
 996
 997   if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
 998
 999
1000   //Get the daughters
1001   if(meson->GetNDaughters() != 2)
1002   {
1003     if(fDebug > 2)
1004       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
1005     return;
1006   }
1007
1008   Int_t iPhoton0 = meson->GetDaughter(0);
1009   Int_t iPhoton1 = meson->GetDaughter(1);
1010   //if((iPhoton0 == -1) || (iPhoton1 == -1)){
1011   //    if(fDebug > 2)
1012   //            printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
1013   //    return;
1014   //}
1015   AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
1016   AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
1017
1018   //Check if both daughters are photons
1019   if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22)
1020   {
1021     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overlapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
1022     return;
1023   }
1024
1025   if(fDebug > 2)
1026     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
1027
1028   //Check if both photons contribute to the cluster
1029   Bool_t okPhoton0 = kFALSE;
1030   Bool_t okPhoton1 = kFALSE;
1031
1032   if(fDebug > 3)
1033     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
1034
1035   Bool_t conversion = kFALSE;
1036
1037   for(Int_t i = 0; i < nlabels; i++)
1038   {
1039     if(fDebug > 3)
1040       printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
1041
1042     if(labels[i]<0) continue;
1043
1044     //If we already found both, break the loop
1045     if(okPhoton0 && okPhoton1) break;
1046
1047     Int_t index =       labels[i];
1048     if      (iPhoton0 == index)
1049     {
1050       okPhoton0 = kTRUE;
1051       continue;
1052     }
1053     else if (iPhoton1 == index)
1054     {
1055       okPhoton1 = kTRUE;
1056       continue;
1057     }
1058
1059     //Trace back the mother in case it was a conversion
1060
1061     if(index >= mcparticles->GetEntriesFast())
1062     {
1063       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast());
1064       continue;
1065     }
1066
1067     AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
1068     Int_t tmpindex = daught->GetMother();
1069     if(fDebug > 3)
1070       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
1071
1072     while(tmpindex>=0){
1073
1074       //MC particle of interest is the mother
1075       if(fDebug > 3)
1076         printf("\t parent index %d\n",tmpindex);
1077       daught   = (AliAODMCParticle*) mcparticles->At(tmpindex);
1078       //printf("tmpindex %d\n",tmpindex);
1079       if      (iPhoton0 == tmpindex)
1080       {
1081         conversion = kTRUE;
1082         okPhoton0  = kTRUE;
1083         break;
1084       }
1085       else if (iPhoton1 == tmpindex)
1086       {
1087         conversion = kTRUE;
1088         okPhoton1  = kTRUE;
1089         break;
1090       }
1091
1092       tmpindex = daught->GetMother();
1093
1094     }//While to check if pi0/eta daughter was one of these contributors to the cluster
1095
1096     if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
1097       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
1098
1099   }//loop on list of labels
1100
1101   //If both photons contribute tag as the corresponding meson.
1102   if(okPhoton0 && okPhoton1)
1103   {
1104     if(fDebug > 2)
1105       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",
1106              (TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
1107
1108     if(!CheckTagBit(tag,kMCConversion) && conversion)
1109     {
1110       if(fDebug > 2)
1111         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Second decay photon produced a conversion \n");
1112       SetTagBit(tag,kMCConversion) ;
1113     }
1114
1115     if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
1116     else                SetTagBit(tag,kMCEta);
1117   }
1118
1119 }
1120
1121 //______________________________________________________________________________________________________
1122 void    AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster,   Int_t iMom, Int_t iParent,
1123                                                AliStack * stack,                Int_t & tag)
1124 {
1125   // Check on ESDs if the current decay photon has the second photon companion lost
1126
1127   if(!arrayCluster || iMom < 0 || iParent < 0|| !stack) return;
1128
1129   TParticle * parent= stack->Particle(iParent);
1130
1131   if(parent->GetNDaughters()!=2)
1132   {
1133     SetTagBit(tag, kMCDecayPairLost);
1134     return ;
1135   }
1136
1137   Int_t pairLabel = -1;
1138   if     ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
1139   else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
1140
1141   if(pairLabel<0)
1142   {
1143     SetTagBit(tag, kMCDecayPairLost);
1144     return ;
1145   }
1146
1147   for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
1148   {
1149     AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
1150     for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
1151     {
1152       Int_t label = cluster->GetLabels()[ilab];
1153       if(label==pairLabel)
1154       {
1155         SetTagBit(tag, kMCDecayPairInCalo);
1156         return ;
1157       }
1158       else if(label== iParent || label== iMom)
1159       {
1160         continue;
1161       }
1162       else // check the ancestry
1163       {
1164         TParticle * mother = stack->Particle(label);
1165         Int_t momPDG = TMath::Abs(mother->GetPdgCode());
1166         if(momPDG!=11 && momPDG!=22) continue;
1167
1168         //Check if "mother" of entity is converted, if not, get the first non converted mother
1169         Int_t iParentClus = mother->GetFirstMother();
1170         if(iParentClus < 0) continue;
1171
1172         TParticle * parentClus = stack->Particle(iParentClus);
1173         if(!parentClus) continue;
1174
1175         Int_t parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
1176         Int_t parentClusStatus = parentClus->GetStatusCode();
1177
1178         if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0) continue;
1179
1180         //printf("Conversion\n");
1181
1182         //Check if the mother is photon or electron with status not stable
1183         while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
1184         {
1185           //New Mother
1186           label            = iParentClus;
1187           momPDG           = parentClusPDG;
1188
1189           iParentClus      = parentClus->GetFirstMother();
1190           if(iParentClus < 0) break;
1191
1192           parentClus       = stack->Particle(iParentClus);
1193           if(!parentClus) break;
1194
1195           parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
1196           parentClusStatus = parentClus->GetStatusCode() ;
1197         }//while
1198
1199         if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel))
1200         {
1201           SetTagBit(tag, kMCDecayPairInCalo);
1202           //printf("Conversion is paired\n");
1203           return ;
1204         }
1205         else continue;
1206
1207       }
1208     }
1209   } // cluster loop
1210
1211   SetTagBit(tag, kMCDecayPairLost);
1212
1213 }
1214
1215 //______________________________________________________________________________________________________
1216 void    AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray   * arrayCluster,Int_t iMom, Int_t iParent,
1217                                                const TClonesArray* mcparticles, Int_t & tag)
1218 {
1219   // Check on AODs if the current decay photon has the second photon companion lost
1220
1221   if(!arrayCluster || iMom < 0 || iParent < 0|| !mcparticles) return;
1222
1223   AliAODMCParticle * parent = (AliAODMCParticle*) mcparticles->At(iParent);
1224
1225   //printf("*** Check label %d with parent %d\n",iMom, iParent);
1226
1227   if(parent->GetNDaughters()!=2)
1228   {
1229     SetTagBit(tag, kMCDecayPairLost);
1230     //printf("\t ndaugh = %d\n",parent->GetNDaughters());
1231     return ;
1232   }
1233
1234   Int_t pairLabel = -1;
1235   if     ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
1236   else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
1237
1238   if(pairLabel<0)
1239   {
1240     //printf("\t pair Label not found = %d\n",pairLabel);
1241     SetTagBit(tag, kMCDecayPairLost);
1242     return ;
1243   }
1244
1245   //printf("\t *** find pair %d\n",pairLabel);
1246
1247   for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
1248   {
1249     AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
1250     //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
1251     for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
1252     {
1253       Int_t label = cluster->GetLabels()[ilab];
1254
1255       //printf("\t \t label %d\n",label);
1256
1257       if(label==pairLabel)
1258       {
1259         //printf("\t \t Pair found\n");
1260         SetTagBit(tag, kMCDecayPairInCalo);
1261         return ;
1262       }
1263       else if(label== iParent || label== iMom)
1264       {
1265         //printf("\t \t skip\n");
1266         continue;
1267       }
1268       else // check the ancestry
1269       {
1270         AliAODMCParticle * mother = (AliAODMCParticle*) mcparticles->At(label);
1271         Int_t momPDG = TMath::Abs(mother->GetPdgCode());
1272         if(momPDG!=11 && momPDG!=22)
1273         {
1274           //printf("\t \t skip, pdg %d\n",momPDG);
1275           continue;
1276         }
1277
1278         //Check if "mother" of entity is converted, if not, get the first non converted mother
1279         Int_t iParentClus = mother->GetMother();
1280         if(iParentClus < 0) continue;
1281
1282         AliAODMCParticle * parentClus =  (AliAODMCParticle*) mcparticles->At(iParentClus);
1283         if(!parentClus) continue;
1284
1285         Int_t parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
1286         Int_t parentClusStatus = parentClus->GetStatus();
1287
1288         if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0)
1289         {
1290           //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
1291           continue;
1292         }
1293
1294         //printf("\t \t Conversion\n");
1295
1296         //Check if the mother is photon or electron with status not stable
1297         while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
1298         {
1299           //New Mother
1300           label            = iParentClus;
1301           momPDG           = parentClusPDG;
1302
1303           iParentClus      = parentClus->GetMother();
1304           if(iParentClus < 0) break;
1305
1306           parentClus       =  (AliAODMCParticle*) mcparticles->At(iParentClus);
1307           if(!parentClus) break;
1308
1309           parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
1310           parentClusStatus = parentClus->GetStatus() ;
1311         }//while
1312
1313         if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel))
1314         {
1315           SetTagBit(tag, kMCDecayPairInCalo);
1316           //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
1317           return ;
1318         }
1319         else
1320         {
1321           //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
1322           continue;
1323         }
1324
1325       }
1326     }
1327   } // cluster loop
1328
1329
1330   SetTagBit(tag, kMCDecayPairLost);
1331
1332
1333 }
1334
1335 //_____________________________________________________________________
1336 TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
1337 {
1338   //Return list of jets (TParticles) and index of most likely parton that originated it.
1339   AliStack * stack = reader->GetStack();
1340   Int_t iEvent = reader->GetEventNumber();
1341   AliGenEventHeader * geh = reader->GetGenEventHeader();
1342   if(fCurrentEvent!=iEvent){
1343     fCurrentEvent = iEvent;
1344     fJetsList = new TList;
1345     Int_t nTriggerJets = 0;
1346     Float_t tmpjet[]={0,0,0,0};
1347
1348     //printf("Event %d %d\n",fCurrentEvent,iEvent);
1349     //Get outgoing partons
1350     if(stack->GetNtrack() < 8) return fJetsList;
1351     TParticle * parton1 =  stack->Particle(6);
1352     TParticle * parton2 =  stack->Particle(7);
1353     if(fDebug > 2){
1354       printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1355              parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
1356       printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1357              parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
1358                 }
1359     //          //Trace the jet from the mother parton
1360     //          Float_t pt  = 0;
1361     //          Float_t pt1 = 0;
1362     //          Float_t pt2 = 0;
1363     //          Float_t e   = 0;
1364     //          Float_t e1  = 0;
1365     //          Float_t e2  = 0;
1366     //          TParticle * tmptmp = new TParticle;
1367     //          for(Int_t i = 0; i< stack->GetNprimary(); i++){
1368     //                  tmptmp = stack->Particle(i);
1369
1370     //                  if(tmptmp->GetStatusCode() == 1){
1371     //                          pt = tmptmp->Pt();
1372     //                          e =  tmptmp->Energy();
1373     //                          Int_t imom = tmptmp->GetFirstMother();
1374     //                          Int_t imom1 = 0;
1375     //                          //printf("1st imom %d\n",imom);
1376     //                          while(imom > 5){
1377     //                                  imom1=imom;
1378     //                                  tmptmp = stack->Particle(imom);
1379     //                                  imom = tmptmp->GetFirstMother();
1380     //                                  //printf("imom %d       \n",imom);
1381     //                          }
1382     //                          //printf("Last imom %d %d\n",imom1, imom);
1383     //                          if(imom1 == 6) {
1384     //                                  pt1+=pt;
1385     //                                  e1+=e;
1386     //                          }
1387     //                          else if (imom1 == 7){
1388     //                                  pt2+=pt;
1389     //                                  e2+=e;                                  }
1390     //                  }// status 1
1391
1392     //          }// for
1393
1394     //          printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
1395
1396                 //Get the jet, different way for different generator
1397                 //PYTHIA
1398     if(fMCGenerator == "PYTHIA")
1399     {
1400       TParticle * jet =  0x0;
1401       AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
1402       nTriggerJets =  pygeh->NTriggerJets();
1403       if(fDebug > 1)
1404         printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
1405
1406       for(Int_t i = 0; i< nTriggerJets; i++)
1407       {
1408         pygeh->TriggerJet(i, tmpjet);
1409         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1410         //Assign an outgoing parton as mother
1411         Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1412         Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1413         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1414         else  jet->SetFirstMother(6);
1415         //jet->Print();
1416         if(fDebug > 1)
1417           printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1418                  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1419         fJetsList->Add(jet);
1420       }
1421     }//Pythia triggered jets
1422     //HERWIG
1423     else if (fMCGenerator=="HERWIG")
1424     {
1425       Int_t pdg = -1;
1426       //Check parton 1
1427       TParticle * tmp = parton1;
1428       if(parton1->GetPdgCode()!=22)
1429       {
1430         while(pdg != 94){
1431           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1432           tmp = stack->Particle(tmp->GetFirstDaughter());
1433           pdg = tmp->GetPdgCode();
1434         }//while
1435
1436         //Add found jet to list
1437         TParticle *jet1 = new TParticle(*tmp);
1438         jet1->SetFirstMother(6);
1439         fJetsList->Add(jet1);
1440         //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1441         //tmp = stack->Particle(tmp->GetFirstDaughter());
1442         //tmp->Print();
1443         //jet1->Print();
1444         if(fDebug > 1)
1445           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1446                  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
1447       }//not photon
1448
1449       //Check parton 2
1450       pdg = -1;
1451       tmp = parton2;
1452       if(parton2->GetPdgCode()!=22)
1453       {
1454         while(pdg != 94)
1455         {
1456           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1457           tmp = stack->Particle(tmp->GetFirstDaughter());
1458           pdg = tmp->GetPdgCode();
1459         }//while
1460
1461         //Add found jet to list
1462         TParticle *jet2 = new TParticle(*tmp);
1463         jet2->SetFirstMother(7);
1464         fJetsList->Add(jet2);
1465         //jet2->Print();
1466         if(fDebug > 1)
1467           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1468                  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1469         //Int_t first =  tmp->GetFirstDaughter();
1470         //Int_t last  =  tmp->GetLastDaughter();
1471         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1472                                 //      for(Int_t d = first ; d < last+1; d++){
1473         //                                              tmp = stack->Particle(d);
1474         //                                              if(i == tmp->GetFirstMother())
1475         //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1476         //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1477         //                         }
1478         //tmp->Print();
1479       }//not photon
1480     }//Herwig generated jets
1481   }
1482
1483   return fJetsList;
1484 }
1485
1486
1487 //__________________________________________________________________________________________________________
1488 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
1489                                                const AliCaloTrackReader* reader,
1490                                                Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel)
1491 {
1492   //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1493
1494   TLorentzVector daughter(0,0,0,0);
1495
1496   if(reader->ReadStack())
1497   {
1498     if(!reader->GetStack())
1499     {
1500       if (fDebug >=0)
1501         printf("AliMCAnalysisUtils::GetDaughter() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1502
1503       ok=kFALSE;
1504       return daughter;
1505     }
1506
1507     Int_t nprimaries = reader->GetStack()->GetNtrack();
1508     if(label >= 0 && label < nprimaries)
1509     {
1510       TParticle * momP = reader->GetStack()->Particle(label);
1511       daughlabel       = momP->GetDaughter(idaugh);
1512
1513       if(daughlabel < 0 || daughlabel >= nprimaries)
1514       {
1515         ok = kFALSE;
1516         return daughter;
1517       }
1518
1519       TParticle * daughP = reader->GetStack()->Particle(daughlabel);
1520       daughP->Momentum(daughter);
1521       pdg    = daughP->GetPdgCode();
1522       status = daughP->GetStatusCode();
1523     }
1524     else
1525     {
1526       ok = kFALSE;
1527       return daughter;
1528     }
1529   }
1530   else if(reader->ReadAODMCParticles())
1531   {
1532     TClonesArray* mcparticles = reader->GetAODMCParticles();
1533     if(!mcparticles)
1534     {
1535       if(fDebug >= 0)
1536         printf("AliMCAnalysisUtils::GetDaughter() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1537
1538       ok=kFALSE;
1539       return daughter;
1540     }
1541
1542     Int_t nprimaries = mcparticles->GetEntriesFast();
1543     if(label >= 0 && label < nprimaries)
1544     {
1545       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1546       daughlabel              = momP->GetDaughter(idaugh);
1547
1548       if(daughlabel < 0 || daughlabel >= nprimaries)
1549       {
1550         ok = kFALSE;
1551         return daughter;
1552       }
1553
1554       AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
1555       daughter.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1556       pdg    = daughP->GetPdgCode();
1557       status = daughP->GetStatus();
1558     }
1559     else
1560     {
1561       ok = kFALSE;
1562       return daughter;
1563     }
1564   }
1565
1566   ok = kTRUE;
1567
1568   return daughter;
1569 }
1570
1571 //______________________________________________________________________________________________________
1572 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1573 {
1574   //Return the kinematics of the particle that generated the signal
1575
1576   Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1577   return GetMother(label,reader,pdg,status, ok,momlabel);
1578 }
1579
1580 //_________________________________________________________________________________________
1581 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1582                                              Int_t & pdg, Int_t & status, Bool_t & ok)
1583 {
1584   //Return the kinematics of the particle that generated the signal
1585
1586   Int_t momlabel = -1;
1587   return GetMother(label,reader,pdg,status, ok,momlabel);
1588 }
1589
1590 //______________________________________________________________________________________________________
1591 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1592                                              Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1593 {
1594   //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1595
1596   TLorentzVector mom(0,0,0,0);
1597
1598   if(reader->ReadStack())
1599   {
1600     if(!reader->GetStack())
1601     {
1602       if (fDebug >=0)
1603         printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1604
1605       ok=kFALSE;
1606       return mom;
1607     }
1608     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1609     {
1610       TParticle * momP = reader->GetStack()->Particle(label);
1611       momP->Momentum(mom);
1612       pdg    = momP->GetPdgCode();
1613       status = momP->GetStatusCode();
1614       momlabel = momP->GetFirstMother();
1615     }
1616     else
1617     {
1618       ok = kFALSE;
1619       return mom;
1620     }
1621   }
1622   else if(reader->ReadAODMCParticles())
1623   {
1624     TClonesArray* mcparticles = reader->GetAODMCParticles();
1625     if(!mcparticles)
1626     {
1627       if(fDebug >= 0)
1628         printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1629
1630       ok=kFALSE;
1631       return mom;
1632     }
1633
1634     Int_t nprimaries = mcparticles->GetEntriesFast();
1635     if(label >= 0 && label < nprimaries)
1636     {
1637       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1638       mom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1639       pdg    = momP->GetPdgCode();
1640       status = momP->GetStatus();
1641       momlabel = momP->GetMother();
1642     }
1643     else
1644     {
1645       ok = kFALSE;
1646       return mom;
1647     }
1648   }
1649
1650   ok = kTRUE;
1651
1652   return mom;
1653 }
1654
1655 //___________________________________________________________________________________
1656 TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1657                                                     const AliCaloTrackReader* reader,
1658                                                     Bool_t & ok, Int_t & momlabel)
1659 {
1660   //Return the kinematics of the particle that generated the signal
1661
1662   TLorentzVector grandmom(0,0,0,0);
1663
1664
1665   if(reader->ReadStack())
1666   {
1667     if(!reader->GetStack())
1668     {
1669       if (fDebug >=0)
1670         printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1671
1672       ok = kFALSE;
1673       return grandmom;
1674     }
1675
1676     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1677     {
1678       TParticle * momP = reader->GetStack()->Particle(label);
1679
1680       Int_t grandmomLabel = momP->GetFirstMother();
1681       Int_t grandmomPDG   = -1;
1682       TParticle * grandmomP = 0x0;
1683       while (grandmomLabel >=0 )
1684       {
1685         grandmomP   = reader->GetStack()->Particle(grandmomLabel);
1686         grandmomPDG = grandmomP->GetPdgCode();
1687         if(grandmomPDG==pdg)
1688         {
1689           momlabel = grandmomLabel;
1690           grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1691           break;
1692         }
1693
1694         grandmomLabel =  grandmomP->GetFirstMother();
1695
1696       }
1697
1698       if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1699     }
1700   }
1701   else if(reader->ReadAODMCParticles())
1702   {
1703     TClonesArray* mcparticles = reader->GetAODMCParticles();
1704     if(!mcparticles)
1705     {
1706       if(fDebug >= 0)
1707         printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1708
1709       ok=kFALSE;
1710       return grandmom;
1711     }
1712
1713     Int_t nprimaries = mcparticles->GetEntriesFast();
1714     if(label >= 0 && label < nprimaries)
1715     {
1716       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1717
1718       Int_t grandmomLabel = momP->GetMother();
1719       Int_t grandmomPDG   = -1;
1720       AliAODMCParticle * grandmomP = 0x0;
1721       while (grandmomLabel >=0 )
1722       {
1723         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1724         grandmomPDG = grandmomP->GetPdgCode();
1725         if(grandmomPDG==pdg)
1726         {
1727           //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1728           momlabel = grandmomLabel;
1729           grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1730           break;
1731         }
1732
1733         grandmomLabel =  grandmomP->GetMother();
1734
1735       }
1736
1737       if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
1738
1739     }
1740   }
1741
1742   ok = kTRUE;
1743
1744   return grandmom;
1745 }
1746
1747 //______________________________________________________________________________________________
1748 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
1749                                                   Int_t & pdg, Int_t & status, Bool_t & ok,
1750                                                   Int_t & grandMomLabel, Int_t & greatMomLabel)
1751 {
1752   //Return the kinematics of the particle that generated the signal
1753
1754   TLorentzVector grandmom(0,0,0,0);
1755
1756   if(reader->ReadStack())
1757   {
1758     if(!reader->GetStack())
1759     {
1760       if (fDebug >=0)
1761         printf("AliMCAnalysisUtils::GetGrandMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1762
1763       ok = kFALSE;
1764       return grandmom;
1765     }
1766
1767     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1768     {
1769       TParticle * momP = reader->GetStack()->Particle(label);
1770
1771       grandMomLabel = momP->GetFirstMother();
1772
1773       TParticle * grandmomP = 0x0;
1774
1775       if (grandMomLabel >=0 )
1776       {
1777         grandmomP   = reader->GetStack()->Particle(grandMomLabel);
1778         pdg    = grandmomP->GetPdgCode();
1779         status = grandmomP->GetStatusCode();
1780
1781         grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1782         greatMomLabel =  grandmomP->GetFirstMother();
1783
1784       }
1785     }
1786   }
1787   else if(reader->ReadAODMCParticles())
1788   {
1789     TClonesArray* mcparticles = reader->GetAODMCParticles();
1790     if(!mcparticles)
1791     {
1792       if(fDebug >= 0)
1793         printf("AliMCAnalysisUtils::GetGrandMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1794
1795       ok=kFALSE;
1796       return grandmom;
1797     }
1798
1799     Int_t nprimaries = mcparticles->GetEntriesFast();
1800     if(label >= 0 && label < nprimaries)
1801     {
1802       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1803
1804       grandMomLabel = momP->GetMother();
1805
1806       AliAODMCParticle * grandmomP = 0x0;
1807
1808       if(grandMomLabel >=0 )
1809       {
1810         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1811         pdg    = grandmomP->GetPdgCode();
1812         status = grandmomP->GetStatus();
1813
1814         grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1815         greatMomLabel =  grandmomP->GetMother();
1816
1817       }
1818     }
1819   }
1820
1821   ok = kTRUE;
1822
1823   return grandmom;
1824 }
1825
1826 //_______________________________________________________________________________________________________________
1827 void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader,
1828                                                         Float_t & asym, Float_t & angle, Bool_t & ok)
1829 {
1830   //In case of an eta or pi0 decay into 2 photons, get the asymmetry  in the energy of the photons
1831
1832
1833   if(reader->ReadStack())
1834   {
1835     if(!reader->GetStack())
1836     {
1837       if (fDebug >=0)
1838         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1839
1840       ok = kFALSE;
1841     }
1842     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1843     {
1844       TParticle * momP = reader->GetStack()->Particle(label);
1845
1846       Int_t grandmomLabel = momP->GetFirstMother();
1847       Int_t grandmomPDG   = -1;
1848       TParticle * grandmomP = 0x0;
1849       while (grandmomLabel >=0 )
1850       {
1851         grandmomP   = reader->GetStack()->Particle(grandmomLabel);
1852         grandmomPDG = grandmomP->GetPdgCode();
1853
1854         if(grandmomPDG==pdg) break;
1855
1856         grandmomLabel =  grandmomP->GetFirstMother();
1857
1858       }
1859
1860       if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1861       {
1862         TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1863         TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1864         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1865         {
1866           asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
1867           TLorentzVector lvd1, lvd2;
1868           d1->Momentum(lvd1);
1869           d2->Momentum(lvd2);
1870           angle = lvd1.Angle(lvd2.Vect());
1871         }
1872       }
1873       else
1874       {
1875         ok=kFALSE;
1876         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1877       }
1878
1879       } // good label
1880   }
1881   else if(reader->ReadAODMCParticles())
1882   {
1883     TClonesArray* mcparticles = reader->GetAODMCParticles();
1884     if(!mcparticles)
1885     {
1886       if(fDebug >= 0)
1887         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1888
1889       ok=kFALSE;
1890       return;
1891     }
1892
1893     Int_t nprimaries = mcparticles->GetEntriesFast();
1894     if(label >= 0 && label < nprimaries)
1895     {
1896       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1897
1898       Int_t grandmomLabel = momP->GetMother();
1899       Int_t grandmomPDG   = -1;
1900       AliAODMCParticle * grandmomP = 0x0;
1901       while (grandmomLabel >=0 )
1902       {
1903         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1904         grandmomPDG = grandmomP->GetPdgCode();
1905
1906         if(grandmomPDG==pdg) break;
1907
1908         grandmomLabel =  grandmomP->GetMother();
1909
1910       }
1911
1912       if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1913       {
1914         AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1915         AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1916         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1917         {
1918           asym = (d1->E()-d2->E())/grandmomP->E();
1919           TLorentzVector lvd1, lvd2;
1920           lvd1.SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1921           lvd2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1922           angle = lvd1.Angle(lvd2.Vect());
1923         }
1924       }
1925       else
1926       {
1927         ok=kFALSE;
1928         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg);
1929       }
1930
1931     } // good label
1932   }
1933
1934   ok = kTRUE;
1935
1936 }
1937
1938
1939 //_________________________________________________________________________________________________
1940 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1941 {
1942   // Return the the number of daughters of a given MC particle
1943
1944
1945   if(reader->ReadStack())
1946   {
1947     if(!reader->GetStack())
1948     {
1949       if (fDebug >=0)
1950         printf("AliMCAnalysisUtils::GetNDaughters() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1951
1952       ok=kFALSE;
1953       return -1;
1954     }
1955     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1956     {
1957       TParticle * momP = reader->GetStack()->Particle(label);
1958       ok=kTRUE;
1959       return momP->GetNDaughters();
1960     }
1961     else
1962     {
1963       ok = kFALSE;
1964       return -1;
1965     }
1966   }
1967   else if(reader->ReadAODMCParticles())
1968   {
1969     TClonesArray* mcparticles = reader->GetAODMCParticles();
1970     if(!mcparticles)
1971     {
1972       if(fDebug >= 0)
1973         printf("AliMCAnalysisUtils::GetNDaughters() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1974
1975       ok=kFALSE;
1976       return -1;
1977     }
1978
1979     Int_t nprimaries = mcparticles->GetEntriesFast();
1980     if(label >= 0 && label < nprimaries)
1981     {
1982       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1983       ok = kTRUE;
1984       return momP->GetNDaughters();
1985     }
1986     else
1987     {
1988       ok = kFALSE;
1989       return -1;
1990     }
1991   }
1992
1993   ok = kFALSE;
1994
1995   return -1;
1996 }
1997
1998 //_________________________________________________________________________________
1999 Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
2000                                        Int_t mctag, Int_t mesonLabel,
2001                                        AliCaloTrackReader * reader, Int_t *overpdg)
2002 {
2003   // Compare the primary depositing more energy with the rest,
2004   // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing
2005   // Give as input the meson label in case it was a pi0 or eta merged cluster
2006   // Init overpdg with nlabels
2007
2008   Int_t ancPDG = 0, ancStatus = -1;
2009   TLorentzVector momentum; TVector3 prodVertex;
2010   Int_t ancLabel = 0;
2011   Int_t noverlaps = 0;
2012   Bool_t ok = kFALSE;
2013
2014   for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
2015   {
2016     ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,momentum,prodVertex);
2017
2018     //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
2019     //       ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
2020
2021     Bool_t overlap = kFALSE;
2022
2023     if     ( ancLabel < 0 )
2024     {
2025       overlap = kTRUE;
2026       //printf("\t \t \t No Label = %d\n",ancLabel);
2027     }
2028     else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) ||  CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
2029     {
2030       //printf("\t \t  meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
2031       overlap = kTRUE;
2032     }
2033     else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
2034     {
2035       //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
2036       overlap = kTRUE ;
2037     }
2038
2039     if( !overlap ) continue ;
2040
2041     // We have at least one overlap
2042
2043     //printf("Overlap!!!!!!!!!!!!!!\n");
2044
2045     noverlaps++;
2046
2047     // What is the origin of the overlap?
2048     Bool_t  mOK = 0,      gOK = 0;
2049     Int_t   mpdg = -999999,  gpdg = -1;
2050     Int_t   mstatus = -1, gstatus = -1;
2051     Int_t   gLabel = -1, ggLabel = -1;
2052     TLorentzVector mother      = GetMother     (label[ilab],reader,mpdg,mstatus,mOK);
2053     TLorentzVector grandmother = GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
2054
2055     //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
2056
2057     if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
2058         ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
2059        gLabel >=0 )
2060     {
2061       Int_t labeltmp = gLabel;
2062       while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
2063       {
2064         mpdg=gpdg;
2065         grandmother = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
2066         labeltmp=gLabel;
2067       }
2068     }
2069     overpdg[noverlaps-1] = mpdg;
2070   }
2071
2072   return noverlaps ;
2073
2074 }
2075
2076 //________________________________________________________
2077 void AliMCAnalysisUtils::Print(const Option_t * opt) const
2078 {
2079   //Print some relevant parameters set for the analysis
2080
2081   if(! opt)
2082     return;
2083
2084   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2085
2086   printf("Debug level    = %d\n",fDebug);
2087   printf("MC Generator   = %s\n",fMCGenerator.Data());
2088   printf(" \n");
2089
2090 }
2091
2092 //__________________________________________________
2093 void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
2094 {
2095   // print the assigned origins to this particle
2096
2097   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, unk %d, bad %d\n",
2098          tag,
2099          CheckTagBit(tag,kMCPhoton),
2100          CheckTagBit(tag,kMCConversion),
2101          CheckTagBit(tag,kMCPrompt),
2102          CheckTagBit(tag,kMCFragmentation),
2103          CheckTagBit(tag,kMCISR),
2104          CheckTagBit(tag,kMCPi0Decay),
2105          CheckTagBit(tag,kMCEtaDecay),
2106          CheckTagBit(tag,kMCOtherDecay),
2107          CheckTagBit(tag,kMCPi0),
2108          CheckTagBit(tag,kMCEta),
2109          CheckTagBit(tag,kMCElectron),
2110          CheckTagBit(tag,kMCMuon),
2111          CheckTagBit(tag,kMCPion),
2112          CheckTagBit(tag,kMCProton),
2113          CheckTagBit(tag,kMCAntiNeutron),
2114          CheckTagBit(tag,kMCKaon),
2115          CheckTagBit(tag,kMCAntiProton),
2116          CheckTagBit(tag,kMCAntiNeutron),
2117          CheckTagBit(tag,kMCUnknown),
2118          CheckTagBit(tag,kMCBadLabel)
2119          );
2120 }
2121
2122