PWG/CaloTrackCorrBase/AliMCAnalysisUtils.cxx

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