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