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