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