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