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       Int_t iparton = -1;
1166       for(Int_t i = 0; i< nTriggerJets; i++){
1167         iparton=-1;
1168         pygeh->TriggerJet(i, tmpjet);
1169         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1170         //Assign an outgoing parton as mother
1171         Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1172         Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1173         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1174         else  jet->SetFirstMother(6);
1175         //jet->Print();
1176         if(fDebug > 1)
1177           printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1178                  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1179         fJetsList->Add(jet);
1180       }
1181     }//Pythia triggered jets
1182     //HERWIG
1183     else if (fMCGenerator=="HERWIG"){
1184       Int_t pdg = -1;
1185       //Check parton 1
1186       TParticle * tmp = parton1;
1187       if(parton1->GetPdgCode()!=22){
1188         while(pdg != 94){
1189           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1190           tmp = stack->Particle(tmp->GetFirstDaughter());
1191           pdg = tmp->GetPdgCode();
1192         }//while
1193
1194         //Add found jet to list
1195         TParticle *jet1 = new TParticle(*tmp);
1196         jet1->SetFirstMother(6);
1197         fJetsList->Add(jet1);
1198         //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1199         //tmp = stack->Particle(tmp->GetFirstDaughter());
1200         //tmp->Print();
1201         //jet1->Print();
1202         if(fDebug > 1)
1203           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1204                  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
1205       }//not photon
1206
1207       //Check parton 2
1208       pdg = -1;
1209       tmp = parton2;
1210       Int_t i = -1;
1211       if(parton2->GetPdgCode()!=22){
1212         while(pdg != 94){
1213           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1214           i = tmp->GetFirstDaughter();
1215           tmp = stack->Particle(tmp->GetFirstDaughter());
1216           pdg = tmp->GetPdgCode();
1217         }//while
1218         //Add found jet to list
1219         TParticle *jet2 = new TParticle(*tmp);
1220         jet2->SetFirstMother(7);
1221         fJetsList->Add(jet2);
1222         //jet2->Print();
1223         if(fDebug > 1)
1224           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1225                  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1226         //Int_t first =  tmp->GetFirstDaughter();
1227         //Int_t last  =  tmp->GetLastDaughter();
1228         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1229                                 //      for(Int_t d = first ; d < last+1; d++){
1230         //                                              tmp = stack->Particle(d);
1231         //                                              if(i == tmp->GetFirstMother())
1232         //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1233         //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1234         //                         }
1235         //tmp->Print();
1236       }//not photon
1237     }//Herwig generated jets
1238   }
1239
1240   return fJetsList;
1241 }
1242
1243
1244 //__________________________________________________________________________________________________________
1245 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
1246                                                const AliCaloTrackReader* reader,
1247                                                Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel)
1248 {
1249   //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1250
1251   TLorentzVector daughter(0,0,0,0);
1252
1253   if(reader->ReadStack())
1254   {
1255     if(!reader->GetStack())
1256     {
1257       if (fDebug >=0)
1258         printf("AliMCAnalysisUtils::GetDaughter() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1259
1260       ok=kFALSE;
1261       return daughter;
1262     }
1263     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1264     {
1265       TParticle * momP = reader->GetStack()->Particle(label);
1266       daughlabel       = momP->GetDaughter(idaugh);
1267
1268       TParticle * daughP = reader->GetStack()->Particle(daughlabel);
1269       daughP->Momentum(daughter);
1270       pdg    = daughP->GetPdgCode();
1271       status = daughP->GetStatusCode();
1272     }
1273     else
1274     {
1275       ok = kFALSE;
1276       return daughter;
1277     }
1278   }
1279   else if(reader->ReadAODMCParticles())
1280   {
1281     TClonesArray* mcparticles = reader->GetAODMCParticles();
1282     if(!mcparticles)
1283     {
1284       if(fDebug >= 0)
1285         printf("AliMCAnalysisUtils::GetDaughter() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1286
1287       ok=kFALSE;
1288       return daughter;
1289     }
1290
1291     Int_t nprimaries = mcparticles->GetEntriesFast();
1292     if(label >= 0 && label < nprimaries)
1293     {
1294       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1295       daughlabel              = momP->GetDaughter(idaugh);
1296
1297       AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
1298       daughter.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1299       pdg    = daughP->GetPdgCode();
1300       status = daughP->GetStatus();
1301     }
1302     else
1303     {
1304       ok = kFALSE;
1305       return daughter;
1306     }
1307   }
1308
1309   ok = kTRUE;
1310
1311   return daughter;
1312 }
1313
1314 //______________________________________________________________________________________________________
1315 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1316 {
1317   //Return the kinematics of the particle that generated the signal
1318
1319   Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1320   return GetMother(label,reader,pdg,status, ok,momlabel);
1321 }
1322
1323 //_________________________________________________________________________________________
1324 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1325                                              Int_t & pdg, Int_t & status, Bool_t & ok)
1326 {
1327   //Return the kinematics of the particle that generated the signal
1328
1329   Int_t momlabel = -1;
1330   return GetMother(label,reader,pdg,status, ok,momlabel);
1331 }
1332
1333 //______________________________________________________________________________________________________
1334 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1335                                              Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1336 {
1337   //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1338
1339   TLorentzVector mom(0,0,0,0);
1340
1341   if(reader->ReadStack())
1342   {
1343     if(!reader->GetStack())
1344     {
1345       if (fDebug >=0)
1346         printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1347
1348       ok=kFALSE;
1349       return mom;
1350     }
1351     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1352     {
1353       TParticle * momP = reader->GetStack()->Particle(label);
1354       momP->Momentum(mom);
1355       pdg    = momP->GetPdgCode();
1356       status = momP->GetStatusCode();
1357       momlabel = momP->GetFirstMother();
1358     }
1359     else
1360     {
1361       ok = kFALSE;
1362       return mom;
1363     }
1364   }
1365   else if(reader->ReadAODMCParticles())
1366   {
1367     TClonesArray* mcparticles = reader->GetAODMCParticles();
1368     if(!mcparticles)
1369     {
1370       if(fDebug >= 0)
1371         printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1372
1373       ok=kFALSE;
1374       return mom;
1375     }
1376
1377     Int_t nprimaries = mcparticles->GetEntriesFast();
1378     if(label >= 0 && label < nprimaries)
1379     {
1380       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1381       mom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1382       pdg    = momP->GetPdgCode();
1383       status = momP->GetStatus();
1384       momlabel = momP->GetMother();
1385     }
1386     else
1387     {
1388       ok = kFALSE;
1389       return mom;
1390     }
1391   }
1392
1393   ok = kTRUE;
1394
1395   return mom;
1396 }
1397
1398 //___________________________________________________________________________________
1399 TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1400                                                     const AliCaloTrackReader* reader,
1401                                                     Bool_t & ok, Int_t & momlabel)
1402 {
1403   //Return the kinematics of the particle that generated the signal
1404
1405   TLorentzVector grandmom(0,0,0,0);
1406
1407
1408   if(reader->ReadStack())
1409   {
1410     if(!reader->GetStack())
1411     {
1412       if (fDebug >=0)
1413         printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1414
1415       ok = kFALSE;
1416       return grandmom;
1417     }
1418
1419     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1420     {
1421       TParticle * momP = reader->GetStack()->Particle(label);
1422
1423       Int_t grandmomLabel = momP->GetFirstMother();
1424       Int_t grandmomPDG   = -1;
1425       TParticle * grandmomP = 0x0;
1426       while (grandmomLabel >=0 )
1427       {
1428         grandmomP   = reader->GetStack()->Particle(grandmomLabel);
1429         grandmomPDG = grandmomP->GetPdgCode();
1430         if(grandmomPDG==pdg)
1431         {
1432           momlabel = grandmomLabel;
1433           grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1434           break;
1435         }
1436
1437         grandmomLabel =  grandmomP->GetFirstMother();
1438
1439       }
1440
1441       if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1442     }
1443   }
1444   else if(reader->ReadAODMCParticles())
1445   {
1446     TClonesArray* mcparticles = reader->GetAODMCParticles();
1447     if(!mcparticles)
1448     {
1449       if(fDebug >= 0)
1450         printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1451
1452       ok=kFALSE;
1453       return grandmom;
1454     }
1455
1456     Int_t nprimaries = mcparticles->GetEntriesFast();
1457     if(label >= 0 && label < nprimaries)
1458     {
1459       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1460
1461       Int_t grandmomLabel = momP->GetMother();
1462       Int_t grandmomPDG   = -1;
1463       AliAODMCParticle * grandmomP = 0x0;
1464       while (grandmomLabel >=0 )
1465       {
1466         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1467         grandmomPDG = grandmomP->GetPdgCode();
1468         if(grandmomPDG==pdg)
1469         {
1470           //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1471           momlabel = grandmomLabel;
1472           grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1473           break;
1474         }
1475
1476         grandmomLabel =  grandmomP->GetMother();
1477
1478       }
1479
1480       if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
1481
1482     }
1483   }
1484
1485   ok = kTRUE;
1486
1487   return grandmom;
1488 }
1489
1490 //______________________________________________________________________________________________
1491 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
1492                                                   Int_t & pdg, Int_t & status, Bool_t & ok,
1493                                                   Int_t & grandMomLabel, Int_t & greatMomLabel)
1494 {
1495   //Return the kinematics of the particle that generated the signal
1496
1497   TLorentzVector grandmom(0,0,0,0);
1498
1499   if(reader->ReadStack())
1500   {
1501     if(!reader->GetStack())
1502     {
1503       if (fDebug >=0)
1504         printf("AliMCAnalysisUtils::GetGrandMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1505
1506       ok = kFALSE;
1507       return grandmom;
1508     }
1509
1510     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1511     {
1512       TParticle * momP = reader->GetStack()->Particle(label);
1513
1514       grandMomLabel = momP->GetFirstMother();
1515
1516       TParticle * grandmomP = 0x0;
1517
1518       if (grandMomLabel >=0 )
1519       {
1520         grandmomP   = reader->GetStack()->Particle(grandMomLabel);
1521         pdg    = grandmomP->GetPdgCode();
1522         status = grandmomP->GetStatusCode();
1523
1524         grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1525         greatMomLabel =  grandmomP->GetFirstMother();
1526
1527       }
1528     }
1529   }
1530   else if(reader->ReadAODMCParticles())
1531   {
1532     TClonesArray* mcparticles = reader->GetAODMCParticles();
1533     if(!mcparticles)
1534     {
1535       if(fDebug >= 0)
1536         printf("AliMCAnalysisUtils::GetGrandMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1537
1538       ok=kFALSE;
1539       return grandmom;
1540     }
1541
1542     Int_t nprimaries = mcparticles->GetEntriesFast();
1543     if(label >= 0 && label < nprimaries)
1544     {
1545       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1546
1547       grandMomLabel = momP->GetMother();
1548
1549       AliAODMCParticle * grandmomP = 0x0;
1550
1551       if(grandMomLabel >=0 )
1552       {
1553         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1554         pdg    = grandmomP->GetPdgCode();
1555         status = grandmomP->GetStatus();
1556
1557         grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1558         greatMomLabel =  grandmomP->GetMother();
1559
1560       }
1561     }
1562   }
1563
1564   ok = kTRUE;
1565
1566   return grandmom;
1567 }
1568
1569 //_______________________________________________________________________________________________________________
1570 void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader,
1571                                                         Float_t & asym, Float_t & angle, Bool_t & ok)
1572 {
1573   //In case of an eta or pi0 decay into 2 photons, get the asymmetry  in the energy of the photons
1574
1575
1576   if(reader->ReadStack())
1577   {
1578     if(!reader->GetStack())
1579     {
1580       if (fDebug >=0)
1581         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1582
1583       ok = kFALSE;
1584     }
1585     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1586     {
1587       TParticle * momP = reader->GetStack()->Particle(label);
1588
1589       Int_t grandmomLabel = momP->GetFirstMother();
1590       Int_t grandmomPDG   = -1;
1591       TParticle * grandmomP = 0x0;
1592       while (grandmomLabel >=0 )
1593       {
1594         grandmomP   = reader->GetStack()->Particle(grandmomLabel);
1595         grandmomPDG = grandmomP->GetPdgCode();
1596
1597         if(grandmomPDG==pdg) break;
1598
1599         grandmomLabel =  grandmomP->GetFirstMother();
1600
1601       }
1602
1603       if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1604       {
1605         TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1606         TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1607         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1608         {
1609           asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
1610           TLorentzVector lvd1, lvd2;
1611           d1->Momentum(lvd1);
1612           d2->Momentum(lvd2);
1613           angle = lvd1.Angle(lvd2.Vect());
1614         }
1615       }
1616       else
1617       {
1618         ok=kFALSE;
1619         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1620       }
1621
1622       } // good label
1623   }
1624   else if(reader->ReadAODMCParticles())
1625   {
1626     TClonesArray* mcparticles = reader->GetAODMCParticles();
1627     if(!mcparticles)
1628     {
1629       if(fDebug >= 0)
1630         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1631
1632       ok=kFALSE;
1633       return;
1634     }
1635
1636     Int_t nprimaries = mcparticles->GetEntriesFast();
1637     if(label >= 0 && label < nprimaries)
1638     {
1639       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1640
1641       Int_t grandmomLabel = momP->GetMother();
1642       Int_t grandmomPDG   = -1;
1643       AliAODMCParticle * grandmomP = 0x0;
1644       while (grandmomLabel >=0 )
1645       {
1646         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1647         grandmomPDG = grandmomP->GetPdgCode();
1648
1649         if(grandmomPDG==pdg) break;
1650
1651         grandmomLabel =  grandmomP->GetMother();
1652
1653       }
1654
1655       if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1656       {
1657         AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1658         AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1659         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1660         {
1661           asym = (d1->E()-d2->E())/grandmomP->E();
1662           TLorentzVector lvd1, lvd2;
1663           lvd1.SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1664           lvd2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1665           angle = lvd1.Angle(lvd2.Vect());
1666         }
1667       }
1668       else
1669       {
1670         ok=kFALSE;
1671         printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg);
1672       }
1673
1674     } // good label
1675   }
1676
1677   ok = kTRUE;
1678
1679 }
1680
1681
1682 //_________________________________________________________________________________________________
1683 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1684 {
1685   // Return the the number of daughters of a given MC particle
1686
1687
1688   if(reader->ReadStack())
1689   {
1690     if(!reader->GetStack())
1691     {
1692       if (fDebug >=0)
1693         printf("AliMCAnalysisUtils::GetNDaughters() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1694
1695       ok=kFALSE;
1696       return -1;
1697     }
1698     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1699     {
1700       TParticle * momP = reader->GetStack()->Particle(label);
1701       ok=kTRUE;
1702       return momP->GetNDaughters();
1703     }
1704     else
1705     {
1706       ok = kFALSE;
1707       return -1;
1708     }
1709   }
1710   else if(reader->ReadAODMCParticles())
1711   {
1712     TClonesArray* mcparticles = reader->GetAODMCParticles();
1713     if(!mcparticles)
1714     {
1715       if(fDebug >= 0)
1716         printf("AliMCAnalysisUtils::GetNDaughters() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1717
1718       ok=kFALSE;
1719       return -1;
1720     }
1721
1722     Int_t nprimaries = mcparticles->GetEntriesFast();
1723     if(label >= 0 && label < nprimaries)
1724     {
1725       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1726       ok = kTRUE;
1727       return momP->GetNDaughters();
1728     }
1729     else
1730     {
1731       ok = kFALSE;
1732       return -1;
1733     }
1734   }
1735
1736   ok = kFALSE;
1737
1738   return -1;
1739 }
1740
1741 //_________________________________________________________________________________
1742 Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
1743                                        Int_t mctag, Int_t mesonLabel,
1744                                        AliCaloTrackReader * reader, Int_t *overpdg)
1745 {
1746   // Compare the primary depositing more energy with the rest,
1747   // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing
1748   // Give as input the meson label in case it was a pi0 or eta merged cluster
1749   // Init overpdg with nlabels
1750
1751   Int_t ancPDG = 0, ancStatus = -1;
1752   TLorentzVector momentum; TVector3 prodVertex;
1753   Int_t ancLabel = 0;
1754   Int_t noverlaps = 0;
1755   Bool_t ok = kFALSE;
1756
1757   for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1758   {
1759     ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,momentum,prodVertex);
1760
1761     //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1762     //       ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1763
1764     Bool_t overlap = kFALSE;
1765
1766     if     ( ancLabel < 0 )
1767     {
1768       overlap = kTRUE;
1769       //printf("\t \t \t No Label = %d\n",ancLabel);
1770     }
1771     else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) ||  CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
1772     {
1773       //printf("\t \t  meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1774       overlap = kTRUE;
1775     }
1776     else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1777     {
1778       //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1779       overlap = kTRUE ;
1780     }
1781
1782     if( !overlap ) continue ;
1783
1784     // We have at least one overlap
1785
1786     //printf("Overlap!!!!!!!!!!!!!!\n");
1787
1788     noverlaps++;
1789
1790     // What is the origin of the overlap?
1791     Bool_t  mOK = 0,      gOK = 0;
1792     Int_t   mpdg = -999999,  gpdg = -1;
1793     Int_t   mstatus = -1, gstatus = -1;
1794     Int_t   gLabel = -1, ggLabel = -1;
1795     TLorentzVector mother      = GetMother     (label[ilab],reader,mpdg,mstatus,mOK);
1796     TLorentzVector grandmother = GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
1797
1798     //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
1799
1800     if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
1801         ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
1802        gLabel >=0 )
1803     {
1804       Int_t labeltmp = gLabel;
1805       while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
1806       {
1807         mpdg=gpdg;
1808         grandmother = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
1809         labeltmp=gLabel;
1810       }
1811     }
1812     overpdg[noverlaps-1] = mpdg;
1813   }
1814
1815   return noverlaps ;
1816
1817 }
1818
1819 //________________________________________________________
1820 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1821 {
1822   //Print some relevant parameters set for the analysis
1823
1824   if(! opt)
1825     return;
1826
1827   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1828
1829   printf("Debug level    = %d\n",fDebug);
1830   printf("MC Generator   = %s\n",fMCGenerator.Data());
1831   printf(" \n");
1832
1833 }
1834
1835 //__________________________________________________
1836 void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
1837 {
1838   // print the assigned origins to this particle
1839
1840   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",
1841          tag,
1842          CheckTagBit(tag,kMCPhoton),
1843          CheckTagBit(tag,kMCConversion),
1844          CheckTagBit(tag,kMCPrompt),
1845          CheckTagBit(tag,kMCFragmentation),
1846          CheckTagBit(tag,kMCISR),
1847          CheckTagBit(tag,kMCPi0Decay),
1848          CheckTagBit(tag,kMCEtaDecay),
1849          CheckTagBit(tag,kMCOtherDecay),
1850          CheckTagBit(tag,kMCPi0),
1851          CheckTagBit(tag,kMCEta),
1852          CheckTagBit(tag,kMCElectron),
1853          CheckTagBit(tag,kMCMuon),
1854          CheckTagBit(tag,kMCPion),
1855          CheckTagBit(tag,kMCProton),
1856          CheckTagBit(tag,kMCAntiNeutron),
1857          CheckTagBit(tag,kMCKaon),
1858          CheckTagBit(tag,kMCAntiProton),
1859          CheckTagBit(tag,kMCAntiNeutron),
1860          CheckTagBit(tag,kMCOther),
1861          CheckTagBit(tag,kMCUnknown),
1862          CheckTagBit(tag,kMCBadLabel)
1863          );
1864 }
1865
1866