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         counter1++;
  96         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
  97         if(mom){
  98           label  = mom->GetMother() ;
  99           label1[counter1]=label;
 100         }
 101         //printf("\t counter %d, label %d\n", counter1,label);
 102       }
 103       //printf("Org label2=%d,\n",label2[0]);
 104       label=label2[0];
 105       if(label < 0) return -1;
 106       while(label > -1 && counter2 < 99){
 107         counter2++;
 108         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 109         if(mom){
 110           label  = mom->GetMother() ;
 111           label2[counter2]=label;
 112         }
 113         //printf("\t counter %d, label %d\n", counter2,label);
 114       }
 115     }//AOD MC
 116     else { //Kine stack from ESDs
 117       AliStack * stack = reader->GetStack();
 118       Int_t label=label1[0];
 119       while(label > -1 && counter1 < 99){
 120         counter1++;
 121         TParticle * mom = stack->Particle(label);
 122         if(mom){
 123           label  = mom->GetFirstMother() ;
 124           label1[counter1]=label;
 125         }
 126         //printf("\t counter %d, label %d\n", counter1,label);
 127       }
 128       //printf("Org label2=%d,\n",label2[0]);
 129       label=label2[0];
 130       while(label > -1 && counter2 < 99){
 131         counter2++;
 132         TParticle * mom = stack->Particle(label);
 133         if(mom){
 134           label  = mom->GetFirstMother() ;
 135           label2[counter2]=label;
 136         }
 137         //printf("\t counter %d, label %d\n", counter2,label);
 138       }
 139     }// Kine stack from ESDs
 140   }//First labels not the same
 141
 142   if((counter1==99 || counter2==99) && fDebug >=0) printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
 143   //printf("CheckAncestor:\n");
 144   Int_t commonparents = 0;
 145   Int_t ancLabel = -1;
 146   //printf("counters %d %d \n",counter1, counter2);
 147   for (Int_t c1 = 0; c1 < counter1; c1++) {
 148     for (Int_t c2 = 0; c2 < counter2; c2++) {
 149       if(label1[c1]==label2[c2] && label1[c1]>-1) {
 150         ancLabel = label1[c1];
 151         commonparents++;
 152         if(reader->ReadAODMCParticles()){
 153           AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]);
 154           if (mom) {
 155             ancPDG    = mom->GetPdgCode();
 156             ancStatus = mom->GetStatus();
 157             momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
 158             prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
 159           }
 160         }
 161         else {
 162           TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
 163           if (mom) {
 164             ancPDG    = mom->GetPdgCode();
 165             ancStatus = mom->GetStatusCode();
 166             mom->Momentum(momentum);
 167             prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz());
 168           }
 169         }
 170         //First ancestor found, end the loops
 171         counter1=0;
 172         counter2=0;
 173       }//Ancestor found
 174     }//second cluster loop
 175   }//first cluster loop
 176
 177   return ancLabel;
 178 }
 179
 180 //_____________________________________________________________________
 181 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label,
 182                                       const Int_t nlabels,
 183                                       const AliCaloTrackReader* reader)
 184 {
 185   //Play with the montecarlo particles if available
 186   Int_t tag = 0;
 187
 188   if(nlabels<=0) {
 189     printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
 190     return kMCBadLabel;
 191   }
 192
 193   //Select where the information is, ESD-galice stack or AOD mcparticles branch
 194   if(reader->ReadStack()){
 195     tag = CheckOriginInStack(label, nlabels, reader->GetStack());
 196   }
 197   else if(reader->ReadAODMCParticles()){
 198     tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles());
 199   }
 200
 201   return tag ;
 202 }
 203
 204 //_____________________________________________________________________
 205 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label,
 206                                       const AliCaloTrackReader* reader)
 207 {
 208   //Play with the montecarlo particles if available
 209   Int_t tag = 0;
 210
 211   if(label<0) {
 212     printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
 213     return kMCBadLabel;
 214   }
 215
 216   Int_t labels[]={label};
 217
 218   //Select where the information is, ESD-galice stack or AOD mcparticles branch
 219   if(reader->ReadStack()){
 220     tag = CheckOriginInStack(labels, 1,reader->GetStack());
 221   }
 222   else if(reader->ReadAODMCParticles()){
 223     tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles());
 224   }
 225
 226   return tag ;
 227 }
 228
 229 //_________________________________________________________________
 230 Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
 231                                              const Int_t nlabels,
 232                                              AliStack* stack)
 233 {
 234   // Play with the MC stack if available. Tag particles depending on their origin.
 235   // Do same things as in CheckOriginInAOD but different input.
 236
 237   //generally speaking, label is the MC label of a reconstructed
 238   //entity (track, cluster, etc) for which we want to know something
 239   //about its heritage, but one can also use it directly with stack
 240   //particles not connected to reconstructed entities
 241
 242   if(!stack) {
 243     if (fDebug >=0)
 244       printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
 245     return -1;
 246   }
 247
 248   Int_t tag = 0;
 249   Int_t label=labels[0];//Most significant particle contributing to the cluster
 250
 251   if(label >= 0 && label < stack->GetNtrack()){
 252     //MC particle of interest is the "mom" of the entity
 253     TParticle * mom = stack->Particle(label);
 254     Int_t iMom     = label;
 255     Int_t mPdgSign = mom->GetPdgCode();
 256     Int_t mPdg     = TMath::Abs(mPdgSign);
 257     Int_t mStatus  = mom->GetStatusCode() ;
 258     Int_t iParent  = mom->GetFirstMother() ;
 259     if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 260
 261     //GrandParent of the entity
 262     TParticle * parent = NULL;
 263     Int_t pPdg = -1;
 264     Int_t pStatus =-1;
 265     if(iParent >= 0){
 266       parent = stack->Particle(iParent);
 267       if(parent){
 268         pPdg = TMath::Abs(parent->GetPdgCode());
 269         pStatus = parent->GetStatusCode();
 270       }
 271     }
 272     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
 273
 274     if(fDebug > 2 ) {
 275       printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
 276       printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 277       printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 278     }
 279
 280     //Check if "mother" of entity is converted, if not, get the first non converted mother
 281     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
 282       SetTagBit(tag,kMCConversion);
 283       //Check if the mother is photon or electron with status not stable
 284       while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
 285         //Mother
 286         iMom     = mom->GetFirstMother();
 287         mom      = stack->Particle(iMom);
 288         mPdgSign = mom->GetPdgCode();
 289         mPdg     = TMath::Abs(mPdgSign);
 290         mStatus  = mom->GetStatusCode() ;
 291         iParent  = mom->GetFirstMother() ;
 292         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 293
 294         //GrandParent
 295         if(iParent >= 0){
 296           parent = stack->Particle(iParent);
 297           if(parent){
 298             pPdg = TMath::Abs(parent->GetPdgCode());
 299             pStatus = parent->GetStatusCode();
 300           }
 301         }
 302         else {// in case of gun/box simulations
 303           pPdg    = 0;
 304           pStatus = 0;
 305           break;
 306         }
 307       }//while
 308       if(fDebug > 2 ) {
 309         printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
 310         printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 311         printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 312       }
 313
 314     }//mother and parent are electron or photon and have status 0
 315     else if((mPdg == 22 || mPdg == 11) && mStatus == 0){
 316       //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
 317       if(pPdg == 2112 ||  pPdg == 211  ||  pPdg == 321 ||
 318          pPdg == 2212 ||  pPdg == 130  ||  pPdg == 13 ) {
 319         SetTagBit(tag,kMCConversion);
 320         iMom     = mom->GetFirstMother();
 321         mom      = stack->Particle(iMom);
 322         mPdgSign = mom->GetPdgCode();
 323         mPdg     = TMath::Abs(mPdgSign);
 324
 325         if(fDebug > 2 ) {
 326           printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
 327           printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 328         }
 329       }//hadron converted
 330
 331       //Comment for the next lines, we do not check the parent of the hadron for the moment.
 332       //iParent =  mom->GetFirstMother() ;
 333       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 334
 335       //GrandParent
 336       //if(iParent >= 0){
 337       //        parent = stack->Particle(iParent);
 338       //        pPdg = TMath::Abs(parent->GetPdgCode());
 339       //}
 340     }
 341     // conversion into electrons/photons checked
 342
 343     //first check for typical charged particles
 344     if     (mPdg     ==    13) SetTagBit(tag,kMCMuon);
 345     else if(mPdg     ==   211) SetTagBit(tag,kMCPion);
 346     else if(mPdg     ==   321) SetTagBit(tag,kMCKaon);
 347     else if(mPdgSign ==  2212) SetTagBit(tag,kMCProton);
 348     else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
 349     else if(mPdgSign ==  2112) SetTagBit(tag,kMCNeutron);
 350     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
 351
 352     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 353     else if(mPdg == 111)  {
 354       SetTagBit(tag,kMCPi0Decay);
 355       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
 356       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 357     }
 358     else if(mPdg == 221) {
 359       SetTagBit(tag,kMCEtaDecay);
 360       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
 361       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
 362     }
 363     //Photons
 364     else if(mPdg == 22){
 365       SetTagBit(tag,kMCPhoton);
 366       if(mStatus == 1){ //undecayed particle
 367         if(fMCGenerator == "PYTHIA"){
 368           if(iParent < 8 && iParent > 5) {//outgoing partons
 369             if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 370             else SetTagBit(tag,kMCFragmentation);
 371           }//Outgoing partons
 372           else  if(iParent <= 5) {
 373             SetTagBit(tag, kMCISR); //Initial state radiation
 374           }
 375           else if(pStatus == 11){//Decay
 376             if(pPdg == 111) {
 377               SetTagBit(tag,kMCPi0Decay);
 378               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 379               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 380             }
 381             else if (pPdg == 221) {
 382               SetTagBit(tag, kMCEtaDecay);
 383               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 384               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 385             }
 386             else SetTagBit(tag,kMCOtherDecay);
 387           }//Decay
 388           else {
 389             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",
 390                                             mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
 391             if(pPdg == 111) {
 392               SetTagBit(tag,kMCPi0Decay);
 393               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 394               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 395             }
 396             else if (pPdg == 221) {
 397               SetTagBit(tag, kMCEtaDecay);
 398               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 399               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 400             }
 401             else SetTagBit(tag,kMCOtherDecay);
 402           }
 403         }//PYTHIA
 404
 405         else if(fMCGenerator == "HERWIG"){
 406           if(pStatus < 197){//Not decay
 407             while(1){
 408               if(parent){
 409                 if(parent->GetFirstMother()<=5) break;
 410                 iParent = parent->GetFirstMother();
 411                 parent=stack->Particle(iParent);
 412                 pStatus= parent->GetStatusCode();
 413                 pPdg = TMath::Abs(parent->GetPdgCode());
 414               } else break;
 415             }//Look for the parton
 416
 417             if(iParent < 8 && iParent > 5) {
 418               if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 419               else SetTagBit(tag,kMCFragmentation);
 420             }
 421             else SetTagBit(tag,kMCISR);//Initial state radiation
 422           }//Not decay
 423           else{//Decay
 424             if(pPdg == 111) {
 425               SetTagBit(tag,kMCPi0Decay);
 426               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
 427               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 428             }
 429             else if (pPdg == 221) {
 430               SetTagBit(tag,kMCEtaDecay);
 431               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
 432               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
 433             }
 434             else SetTagBit(tag,kMCOtherDecay);
 435           }//Decay
 436         }//HERWIG
 437
 438         else SetTagBit(tag,kMCUnknown);
 439
 440       }//Status 1 : created by event generator
 441
 442       else if(mStatus == 0){ // geant
 443         if(pPdg == 111) {
 444           SetTagBit(tag,kMCPi0Decay);
 445           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
 446           CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 447         }
 448         else if (pPdg == 221) {
 449           SetTagBit(tag,kMCEtaDecay);
 450           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
 451           CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
 452         }
 453         else  SetTagBit(tag,kMCOtherDecay);
 454       }//status 0 : geant generated
 455
 456     }//Mother Photon
 457
 458     //Electron check.  Where did that electron come from?
 459     else if(mPdg == 11){ //electron
 460       if(pPdg == 11 && parent){
 461         Int_t iGrandma = parent->GetFirstMother();
 462         if(iGrandma >= 0) {
 463           TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
 464           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 465
 466           if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 467           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 468         }
 469       }
 470       SetTagBit(tag,kMCElectron);
 471       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
 472       if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
 473       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
 474       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
 475       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
 476         if(parent){
 477           Int_t iGrandma = parent->GetFirstMother();
 478           if(iGrandma >= 0) {
 479             TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
 480             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 481             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
 482             else SetTagBit(tag,kMCEFromC); //c-->e
 483           } else SetTagBit(tag,kMCEFromC); //c-->e
 484         }//parent
 485       } else {
 486         //if it is not from any of the above, where is it from?
 487         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 488         else SetTagBit(tag,kMCOtherDecay);
 489         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);
 490       }
 491     }//electron check
 492     //Cluster was made by something else
 493     else {
 494       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
 495       SetTagBit(tag,kMCUnknown);
 496     }
 497   }//Good label value
 498   else{// Bad label
 499
 500     if(label < 0 && (fDebug >= 0))
 501       printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***:  label %d \n", label);
 502     if(label >=  stack->GetNtrack() &&  (fDebug >= 0))
 503       printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 504     SetTagBit(tag,kMCUnknown);
 505   }//Bad label
 506
 507   return tag;
 508
 509 }
 510
 511
 512 //_________________________________________________________________________
 513 Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
 514                                            const Int_t nlabels,
 515                                            const TClonesArray *mcparticles)
 516 {
 517   // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
 518   // Do same things as in CheckOriginInStack but different input.
 519   if(!mcparticles) {
 520     if(fDebug >= 0)
 521       printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
 522     return -1;
 523   }
 524
 525   Int_t tag = 0;
 526   Int_t label=labels[0];//Most significant particle contributing to the cluster
 527
 528   Int_t nprimaries = mcparticles->GetEntriesFast();
 529   if(label >= 0 && label < nprimaries){
 530     //Mother
 531     AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 532     Int_t iMom     = label;
 533     Int_t mPdgSign = mom->GetPdgCode();
 534     Int_t mPdg     = TMath::Abs(mPdgSign);
 535     Int_t iParent  = mom->GetMother() ;
 536     if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 537
 538     //GrandParent
 539     AliAODMCParticle * parent = NULL ;
 540     Int_t pPdg = -1;
 541     if(iParent >= 0){
 542       parent = (AliAODMCParticle *) mcparticles->At(iParent);
 543       pPdg = TMath::Abs(parent->GetPdgCode());
 544     }
 545     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
 546
 547     if(fDebug > 2 ) {
 548       printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
 549       printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 550       if(parent)
 551         printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 552     }
 553
 554     //Check if mother is converted, if not, get the first non converted mother
 555     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
 556       SetTagBit(tag,kMCConversion);
 557       //Check if the mother is photon or electron with status not stable
 558       while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
 559         //Mother
 560         iMom     = mom->GetMother();
 561         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
 562         mPdgSign = mom->GetPdgCode();
 563         mPdg     = TMath::Abs(mPdgSign);
 564         iParent  = mom->GetMother() ;
 565         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 566
 567         //GrandParent
 568         if(iParent >= 0 && parent){
 569           parent = (AliAODMCParticle *) mcparticles->At(iParent);
 570           pPdg = TMath::Abs(parent->GetPdgCode());
 571         }
 572         // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 573         // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 574
 575       }//while
 576
 577       if(fDebug > 2 ) {
 578         printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
 579         printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 580         if(parent)
 581           printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 582       }
 583
 584     }//mother and parent are electron or photon and have status 0 and parent is photon or electron
 585     else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){
 586       //Still a conversion but only one electron/photon generated. Just from hadrons
 587       if(pPdg == 2112 ||  pPdg == 211 ||  pPdg == 321 ||
 588          pPdg == 2212 ||  pPdg == 130 ||  pPdg == 13 ) {
 589         SetTagBit(tag,kMCConversion);
 590         iMom     = mom->GetMother();
 591         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
 592         mPdgSign = mom->GetPdgCode();
 593         mPdg     = TMath::Abs(mPdgSign);
 594
 595         if(fDebug > 2 ) {
 596           printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
 597           printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 598         }
 599       }//hadron converted
 600
 601       //Comment for next lines, we do not check the parent of the hadron for the moment.
 602       //iParent =  mom->GetMother() ;
 603       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 604
 605       //GrandParent
 606       //if(iParent >= 0){
 607       //        parent = (AliAODMCParticle *) mcparticles->At(iParent);
 608       //        pPdg = TMath::Abs(parent->GetPdgCode());
 609       //}
 610     }
 611
 612     //printf("Final mother mPDG %d\n",mPdg);
 613
 614     // conversion into electrons/photons checked
 615
 616     //first check for typical charged particles
 617     if     (mPdg     ==    13) SetTagBit(tag,kMCMuon);
 618     else if(mPdg     ==   211) SetTagBit(tag,kMCPion);
 619     else if(mPdg     ==   321) SetTagBit(tag,kMCKaon);
 620     else if(mPdgSign ==  2212) SetTagBit(tag,kMCProton);
 621     else if(mPdgSign ==  2112) SetTagBit(tag,kMCNeutron);
 622     else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
 623     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
 624
 625     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 626     else if(mPdg == 111)  {
 627       SetTagBit(tag,kMCPi0Decay);
 628       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
 629       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 630     }
 631     else if(mPdg == 221)  {
 632       SetTagBit(tag,kMCEtaDecay);
 633       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
 634       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 635     }
 636     //Photons
 637     else if(mPdg == 22){
 638       SetTagBit(tag,kMCPhoton);
 639       if(mom->IsPhysicalPrimary() && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) //undecayed particle
 640       {
 641         if(iParent < 8 && iParent > 5 ) {//outgoing partons
 642           if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 643           else SetTagBit(tag,kMCFragmentation);
 644         }//Outgoing partons
 645         else if(iParent <= 5 && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) {
 646           SetTagBit(tag, kMCISR); //Initial state radiation
 647         }
 648         else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
 649           if(pPdg == 111){
 650             SetTagBit(tag,kMCPi0Decay);
 651             if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
 652             CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 653           }
 654           else if (pPdg == 221) {
 655             SetTagBit(tag, kMCEtaDecay);
 656             if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
 657             CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 658           }
 659           else SetTagBit(tag,kMCOtherDecay);
 660         }//Decay
 661         else {
 662           if(parent)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",
 663                            mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
 664           SetTagBit(tag,kMCOtherDecay);//Check
 665         }
 666       }//Physical primary
 667       else if(!mom->IsPrimary()){       //Decays
 668         if(pPdg == 111){
 669           SetTagBit(tag,kMCPi0Decay);
 670           if(fDebug > 2 )
 671             printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
 672           CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 673         }
 674         else if (pPdg == 221) {
 675           SetTagBit(tag,kMCEtaDecay);
 676           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n");
 677           CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 678         }
 679         else  SetTagBit(tag,kMCOtherDecay);
 680       }//not primary : geant generated, decays
 681       else  {
 682         //printf("UNKNOWN 1, mom  pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
 683         //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 684         SetTagBit(tag,kMCUnknown);
 685       }
 686     }//Mother Photon
 687
 688     //Electron check.  Where did that electron come from?
 689     else if(mPdg == 11){ //electron
 690       if(pPdg == 11 && parent){
 691         Int_t iGrandma = parent->GetMother();
 692         if(iGrandma >= 0) {
 693           AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
 694           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 695
 696           if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 697           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 698         }
 699       }
 700       SetTagBit(tag,kMCElectron);
 701       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
 702       if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
 703       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
 704       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
 705       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
 706         if(parent){
 707           Int_t iGrandma = parent->GetMother();
 708           if(iGrandma >= 0) {
 709             AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
 710             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 711             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
 712             else SetTagBit(tag,kMCEFromC); //c-hadron decay
 713           } else SetTagBit(tag,kMCEFromC); //c-hadron decay
 714         }//parent
 715       } else { //prompt or other decay
 716         TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
 717         TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
 718         if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg);
 719         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 720         else SetTagBit(tag,kMCOtherDecay);
 721       }
 722     }//electron check
 723     //cluster was made by something else
 724     else {
 725       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg);
 726       SetTagBit(tag,kMCUnknown);
 727     }
 728   }//Good label value
 729   else{//Bad label
 730
 731     if(label < 0 && (fDebug >= 0) )
 732       printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***:  label %d \n", label);
 733     if(label >=  mcparticles->GetEntriesFast() &&  (fDebug >= 0) )
 734       printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***:  label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
 735     SetTagBit(tag,kMCUnknown);
 736
 737   }//Bad label
 738
 739   return tag;
 740
 741 }
 742
 743 //_________________________________________________________________________
 744 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
 745                                                     const Int_t nlabels,
 746                                                     const Int_t mesonIndex,
 747                                                     AliStack *stack,
 748                                                     Int_t &tag)
 749 {
 750   //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
 751
 752   if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
 753     if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 754                           labels[0],stack->GetNtrack(), nlabels);
 755     return;
 756   }
 757
 758   TParticle * meson = stack->Particle(mesonIndex);
 759   Int_t mesonPdg    = meson->GetPdgCode();
 760   if(mesonPdg!=111 && mesonPdg!=221){
 761     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 762     return;
 763   }
 764
 765   if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
 766
 767   //Check if meson decayed into 2 daughters or if both were kept.
 768   if(meson->GetNDaughters() != 2){
 769     if(fDebug > 2)
 770       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 771     return;
 772   }
 773
 774   //Get the daughters
 775   Int_t iPhoton0 = meson->GetDaughter(0);
 776   Int_t iPhoton1 = meson->GetDaughter(1);
 777   TParticle *photon0 = stack->Particle(iPhoton0);
 778   TParticle *photon1 = stack->Particle(iPhoton1);
 779
 780   //Check if both daughters are photons
 781   if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
 782     if(fDebug > 2)
 783       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 784     return;
 785   }
 786
 787   if(fDebug > 2)
 788     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 789
 790   //Check if both photons contribute to the cluster
 791   Bool_t okPhoton0 = kFALSE;
 792   Bool_t okPhoton1 = kFALSE;
 793
 794   if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
 795
 796   for(Int_t i = 0; i < nlabels; i++){
 797     if(fDebug > 3) printf("\t  at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
 798
 799     //If we already found both, break the loop
 800     if(okPhoton0 && okPhoton1) break;
 801
 802     Int_t index =       labels[i];
 803     if      (iPhoton0 == index) {
 804       okPhoton0 = kTRUE;
 805       continue;
 806     }
 807     else if (iPhoton1 == index) {
 808       okPhoton1 = kTRUE;
 809       continue;
 810     }
 811
 812     //Trace back the mother in case it was a conversion
 813
 814     if(index >= stack->GetNtrack()){
 815       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",index,stack->GetNtrack());
 816       continue;
 817     }
 818
 819     TParticle * daught = stack->Particle(index);
 820     Int_t tmpindex = daught->GetFirstMother();
 821     if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
 822     while(tmpindex>=0){
 823       //MC particle of interest is the mother
 824       if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
 825       daught   = stack->Particle(tmpindex);
 826       if      (iPhoton0 == tmpindex) {
 827         okPhoton0 = kTRUE;
 828         break;
 829       }
 830       else if (iPhoton1 == tmpindex) {
 831         okPhoton1 = kTRUE;
 832         break;
 833       }
 834       tmpindex = daught->GetFirstMother();
 835     }//While to check if pi0/eta daughter was one of these contributors to the cluster
 836
 837     if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
 838       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
 839
 840   }//loop on list of labels
 841
 842   //If both photons contribute tag as the corresponding meson.
 843   if(okPhoton0 && okPhoton1){
 844     if(fDebug > 2)
 845       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
 846
 847     if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 848     else  SetTagBit(tag,kMCEta);
 849   }
 850
 851 }
 852
 853 //__________________________________________________________________________________
 854 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
 855                                                     const Int_t nlabels,
 856                                                     const Int_t mesonIndex,
 857                                                     const TClonesArray *mcparticles,
 858                                                     Int_t & tag )
 859 {
 860   //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
 861
 862   if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1)
 863   {
 864     if(fDebug > 2)
 865       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 866              labels[0],mcparticles->GetEntriesFast(), nlabels);
 867     return;
 868   }
 869
 870   AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
 871   Int_t mesonPdg = meson->GetPdgCode();
 872   if(mesonPdg != 111 && mesonPdg != 221)
 873   {
 874     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 875     return;
 876   }
 877
 878   if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
 879
 880
 881   //Get the daughters
 882   if(meson->GetNDaughters() != 2)
 883   {
 884     if(fDebug > 2)
 885       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 886     return;
 887   }
 888
 889   Int_t iPhoton0 = meson->GetDaughter(0);
 890   Int_t iPhoton1 = meson->GetDaughter(1);
 891   //if((iPhoton0 == -1) || (iPhoton1 == -1)){
 892   //    if(fDebug > 2)
 893   //            printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
 894   //    return;
 895   //}
 896   AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
 897   AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
 898
 899   //Check if both daughters are photons
 900   if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22)
 901   {
 902     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overlapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 903     return;
 904   }
 905
 906   if(fDebug > 2)
 907     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 908
 909   //Check if both photons contribute to the cluster
 910   Bool_t okPhoton0 = kFALSE;
 911   Bool_t okPhoton1 = kFALSE;
 912
 913   if(fDebug > 3)
 914     printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
 915
 916   for(Int_t i = 0; i < nlabels; i++)
 917   {
 918     if(fDebug > 3)
 919       printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
 920
 921     if(labels[i]<0) continue;
 922
 923     //If we already found both, break the loop
 924     if(okPhoton0 && okPhoton1) break;
 925
 926     Int_t index =       labels[i];
 927     if      (iPhoton0 == index) {
 928       okPhoton0 = kTRUE;
 929       continue;
 930     }
 931     else if (iPhoton1 == index) {
 932       okPhoton1 = kTRUE;
 933       continue;
 934     }
 935
 936     //Trace back the mother in case it was a conversion
 937
 938     if(index >= mcparticles->GetEntriesFast()){
 939       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast());
 940       continue;
 941     }
 942
 943     AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
 944     Int_t tmpindex = daught->GetMother();
 945     if(fDebug > 3)
 946       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
 947
 948     while(tmpindex>=0){
 949
 950       //MC particle of interest is the mother
 951       if(fDebug > 3)
 952         printf("\t parent index %d\n",tmpindex);
 953       daught   = (AliAODMCParticle*) mcparticles->At(tmpindex);
 954       //printf("tmpindex %d\n",tmpindex);
 955       if      (iPhoton0 == tmpindex) {
 956         okPhoton0 = kTRUE;
 957         break;
 958       }
 959       else if (iPhoton1 == tmpindex) {
 960         okPhoton1 = kTRUE;
 961         break;
 962       }
 963       tmpindex = daught->GetMother();
 964     }//While to check if pi0/eta daughter was one of these contributors to the cluster
 965
 966     if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
 967       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
 968
 969   }//loop on list of labels
 970
 971   //If both photons contribute tag as the corresponding meson.
 972   if(okPhoton0 && okPhoton1){
 973     if(fDebug > 2)
 974       printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
 975
 976     if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 977     else                SetTagBit(tag,kMCEta);
 978   }
 979
 980 }
 981
 982 //_________________________________________________________________________
 983 TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
 984 {
 985   //Return list of jets (TParticles) and index of most likely parton that originated it.
 986   AliStack * stack = reader->GetStack();
 987   Int_t iEvent = reader->GetEventNumber();
 988   AliGenEventHeader * geh = reader->GetGenEventHeader();
 989   if(fCurrentEvent!=iEvent){
 990     fCurrentEvent = iEvent;
 991     fJetsList = new TList;
 992     Int_t nTriggerJets = 0;
 993     Float_t tmpjet[]={0,0,0,0};
 994
 995     //printf("Event %d %d\n",fCurrentEvent,iEvent);
 996     //Get outgoing partons
 997     if(stack->GetNtrack() < 8) return fJetsList;
 998     TParticle * parton1 =  stack->Particle(6);
 999     TParticle * parton2 =  stack->Particle(7);
1000     if(fDebug > 2){
1001       printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1002              parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
1003       printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1004              parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
1005                 }
1006     //          //Trace the jet from the mother parton
1007     //          Float_t pt  = 0;
1008     //          Float_t pt1 = 0;
1009     //          Float_t pt2 = 0;
1010     //          Float_t e   = 0;
1011     //          Float_t e1  = 0;
1012     //          Float_t e2  = 0;
1013     //          TParticle * tmptmp = new TParticle;
1014     //          for(Int_t i = 0; i< stack->GetNprimary(); i++){
1015     //                  tmptmp = stack->Particle(i);
1016
1017     //                  if(tmptmp->GetStatusCode() == 1){
1018     //                          pt = tmptmp->Pt();
1019     //                          e =  tmptmp->Energy();
1020     //                          Int_t imom = tmptmp->GetFirstMother();
1021     //                          Int_t imom1 = 0;
1022     //                          //printf("1st imom %d\n",imom);
1023     //                          while(imom > 5){
1024     //                                  imom1=imom;
1025     //                                  tmptmp = stack->Particle(imom);
1026     //                                  imom = tmptmp->GetFirstMother();
1027     //                                  //printf("imom %d       \n",imom);
1028     //                          }
1029     //                          //printf("Last imom %d %d\n",imom1, imom);
1030     //                          if(imom1 == 6) {
1031     //                                  pt1+=pt;
1032     //                                  e1+=e;
1033     //                          }
1034     //                          else if (imom1 == 7){
1035     //                                  pt2+=pt;
1036     //                                  e2+=e;                                  }
1037     //                  }// status 1
1038
1039     //          }// for
1040
1041     //          printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
1042
1043                 //Get the jet, different way for different generator
1044                 //PYTHIA
1045     if(fMCGenerator == "PYTHIA"){
1046       TParticle * jet =  0x0;
1047       AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
1048       nTriggerJets =  pygeh->NTriggerJets();
1049       if(fDebug > 1)
1050         printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
1051
1052       Int_t iparton = -1;
1053       for(Int_t i = 0; i< nTriggerJets; i++){
1054         iparton=-1;
1055         pygeh->TriggerJet(i, tmpjet);
1056         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1057         //Assign an outgoing parton as mother
1058         Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1059         Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1060         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1061         else  jet->SetFirstMother(6);
1062         //jet->Print();
1063         if(fDebug > 1)
1064           printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1065                  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1066         fJetsList->Add(jet);
1067       }
1068     }//Pythia triggered jets
1069     //HERWIG
1070     else if (fMCGenerator=="HERWIG"){
1071       Int_t pdg = -1;
1072       //Check parton 1
1073       TParticle * tmp = parton1;
1074       if(parton1->GetPdgCode()!=22){
1075         while(pdg != 94){
1076           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1077           tmp = stack->Particle(tmp->GetFirstDaughter());
1078           pdg = tmp->GetPdgCode();
1079         }//while
1080
1081         //Add found jet to list
1082         TParticle *jet1 = new TParticle(*tmp);
1083         jet1->SetFirstMother(6);
1084         fJetsList->Add(jet1);
1085         //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1086         //tmp = stack->Particle(tmp->GetFirstDaughter());
1087         //tmp->Print();
1088         //jet1->Print();
1089         if(fDebug > 1)
1090           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1091                  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
1092       }//not photon
1093
1094       //Check parton 2
1095       pdg = -1;
1096       tmp = parton2;
1097       Int_t i = -1;
1098       if(parton2->GetPdgCode()!=22){
1099         while(pdg != 94){
1100           if(tmp->GetFirstDaughter()==-1) return fJetsList;
1101           i = tmp->GetFirstDaughter();
1102           tmp = stack->Particle(tmp->GetFirstDaughter());
1103           pdg = tmp->GetPdgCode();
1104         }//while
1105         //Add found jet to list
1106         TParticle *jet2 = new TParticle(*tmp);
1107         jet2->SetFirstMother(7);
1108         fJetsList->Add(jet2);
1109         //jet2->Print();
1110         if(fDebug > 1)
1111           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1112                  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1113         //Int_t first =  tmp->GetFirstDaughter();
1114         //Int_t last  =  tmp->GetLastDaughter();
1115         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1116                                 //      for(Int_t d = first ; d < last+1; d++){
1117         //                                              tmp = stack->Particle(d);
1118         //                                              if(i == tmp->GetFirstMother())
1119         //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1120         //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1121         //                         }
1122         //tmp->Print();
1123       }//not photon
1124     }//Herwig generated jets
1125   }
1126
1127   return fJetsList;
1128 }
1129
1130
1131 //_______________________________________________________________________________________________
1132 TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader,
1133                                              Bool_t & ok)
1134 {
1135   //Return the kinematics of the particle that generated the signal
1136
1137   Int_t pdg = -1; Int_t status = -1;
1138   return GetMother(label,reader,pdg,status, ok);
1139 }
1140
1141 //_______________________________________________________________________________________________
1142 TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader,
1143                                              Int_t & pdg, Int_t & status, Bool_t & ok)
1144 {
1145   //Return the kinematics of the particle that generated the signal, its pdg and its status
1146
1147   TLorentzVector mom(0,0,0,0);
1148
1149   if(reader->ReadStack())
1150   {
1151     if(!reader->GetStack())
1152     {
1153       if (fDebug >=0)
1154         printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1155
1156       ok=kFALSE;
1157       return mom;
1158     }
1159     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1160     {
1161       TParticle * momP = reader->GetStack()->Particle(label);
1162       momP->Momentum(mom);
1163       pdg    = momP->GetPdgCode();
1164       status = momP->GetStatusCode();
1165     }
1166     else
1167     {
1168       ok = kFALSE;
1169       return mom;
1170     }
1171   }
1172   else if(reader->ReadAODMCParticles())
1173   {
1174     TClonesArray* mcparticles = reader->GetAODMCParticles();
1175     if(!mcparticles)
1176     {
1177       if(fDebug >= 0)
1178         printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1179
1180       ok=kFALSE;
1181       return mom;
1182     }
1183
1184     Int_t nprimaries = mcparticles->GetEntriesFast();
1185     if(label >= 0 && label < nprimaries)
1186     {
1187       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1188       mom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1189       pdg    = momP->GetPdgCode();
1190       status = momP->GetStatus();
1191     }
1192     else
1193     {
1194       ok = kFALSE;
1195       return mom;
1196     }
1197   }
1198
1199   ok = kTRUE;
1200
1201   return mom;
1202 }
1203
1204
1205 //_____________________________________________________________________________________
1206 TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader, Bool_t & ok)
1207 {
1208   //Return the kinematics of the particle that generated the signal
1209
1210   TLorentzVector grandmom(0,0,0,0);
1211
1212
1213   if(reader->ReadStack())
1214   {
1215     if(!reader->GetStack())
1216     {
1217       if (fDebug >=0)
1218         printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1219
1220       ok = kFALSE;
1221       return grandmom;
1222     }
1223     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1224     {
1225       TParticle * momP = reader->GetStack()->Particle(label);
1226
1227       Int_t grandmomLabel = momP->GetFirstMother();
1228       Int_t grandmomPDG   = -1;
1229       TParticle * grandmomP = 0x0;
1230       while (grandmomLabel >=0 )
1231       {
1232         grandmomP   = reader->GetStack()->Particle(grandmomLabel);
1233         grandmomPDG = grandmomP->GetPdgCode();
1234         if(grandmomPDG==pdg)
1235         {
1236           grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1237           break;
1238         }
1239
1240         grandmomLabel =  grandmomP->GetFirstMother();
1241
1242       }
1243
1244       if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1245     }
1246   }
1247   else if(reader->ReadAODMCParticles())
1248   {
1249     TClonesArray* mcparticles = reader->GetAODMCParticles();
1250     if(!mcparticles)
1251     {
1252       if(fDebug >= 0)
1253         printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1254
1255       ok=kFALSE;
1256       return grandmom;
1257     }
1258
1259     Int_t nprimaries = mcparticles->GetEntriesFast();
1260     if(label >= 0 && label < nprimaries)
1261     {
1262       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1263
1264       Int_t grandmomLabel = momP->GetMother();
1265       Int_t grandmomPDG   = -1;
1266       AliAODMCParticle * grandmomP = 0x0;
1267       while (grandmomLabel >=0 )
1268       {
1269         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1270         grandmomPDG = grandmomP->GetPdgCode();
1271         if(grandmomPDG==pdg)
1272         {
1273           //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1274
1275           grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1276           break;
1277         }
1278
1279         grandmomLabel =  grandmomP->GetMother();
1280
1281       }
1282
1283       if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
1284
1285     }
1286   }
1287
1288   ok = kTRUE;
1289
1290   return grandmom;
1291 }
1292
1293 //_____________________________________________________________________________________
1294 Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader, Bool_t & ok)
1295 {
1296   //In case of an eta or pi0 decay into 2 photons, get the asymmetry  in the energy of the photons
1297
1298   Float_t asym = -2;
1299
1300   if(reader->ReadStack())
1301   {
1302     if(!reader->GetStack())
1303     {
1304       if (fDebug >=0)
1305         printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1306
1307       ok = kFALSE;
1308       return asym;
1309     }
1310     if(label >= 0 && label < reader->GetStack()->GetNtrack())
1311     {
1312       TParticle * momP = reader->GetStack()->Particle(label);
1313
1314       Int_t grandmomLabel = momP->GetFirstMother();
1315       Int_t grandmomPDG   = -1;
1316       TParticle * grandmomP = 0x0;
1317       while (grandmomLabel >=0 )
1318       {
1319         grandmomP   = reader->GetStack()->Particle(grandmomLabel);
1320         grandmomPDG = grandmomP->GetPdgCode();
1321
1322         if(grandmomPDG==pdg) break;
1323
1324         grandmomLabel =  grandmomP->GetFirstMother();
1325
1326       }
1327
1328       if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1329       {
1330         TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1331         TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1332         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1333         {
1334           asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
1335         }
1336       }
1337       else
1338       {
1339         ok=kFALSE;
1340         printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1341       }
1342
1343       } // good label
1344   }
1345   else if(reader->ReadAODMCParticles())
1346   {
1347     TClonesArray* mcparticles = reader->GetAODMCParticles();
1348     if(!mcparticles)
1349     {
1350       if(fDebug >= 0)
1351         printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1352
1353       ok=kFALSE;
1354       return asym;
1355     }
1356
1357     Int_t nprimaries = mcparticles->GetEntriesFast();
1358     if(label >= 0 && label < nprimaries)
1359     {
1360       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1361
1362       Int_t grandmomLabel = momP->GetMother();
1363       Int_t grandmomPDG   = -1;
1364       AliAODMCParticle * grandmomP = 0x0;
1365       while (grandmomLabel >=0 )
1366       {
1367         grandmomP   = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1368         grandmomPDG = grandmomP->GetPdgCode();
1369
1370         if(grandmomPDG==pdg) break;
1371
1372         grandmomLabel =  grandmomP->GetMother();
1373
1374       }
1375
1376       if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1377       {
1378         AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1379         AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1380         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1381         {
1382           asym = (d1->E()-d2->E())/grandmomP->E();
1383         }
1384       }
1385       else
1386       {
1387         ok=kFALSE;
1388         printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, not found! \n",pdg);
1389       }
1390
1391     } // good label
1392   }
1393
1394   ok = kTRUE;
1395
1396   return asym;
1397 }
1398
1399
1400 //________________________________________________________
1401 void AliMCAnalysisUtils::Print(const Option_t * opt) const
1402 {
1403   //Print some relevant parameters set for the analysis
1404
1405   if(! opt)
1406     return;
1407
1408   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1409
1410   printf("Debug level    = %d\n",fDebug);
1411   printf("MC Generator   = %s\n",fMCGenerator.Data());
1412   printf(" \n");
1413
1414 }
1415
1416 //________________________________________________________
1417 void AliMCAnalysisUtils::PrintMCTag(const Int_t tag) const
1418 {
1419   // print the assigned origins to this particle
1420
1421   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",
1422          tag,
1423          CheckTagBit(tag,kMCPhoton),
1424          CheckTagBit(tag,kMCConversion),
1425          CheckTagBit(tag,kMCPrompt),
1426          CheckTagBit(tag,kMCFragmentation),
1427          CheckTagBit(tag,kMCISR),
1428          CheckTagBit(tag,kMCPi0Decay),
1429          CheckTagBit(tag,kMCEtaDecay),
1430          CheckTagBit(tag,kMCOtherDecay),
1431          CheckTagBit(tag,kMCPi0),
1432          CheckTagBit(tag,kMCEta),
1433          CheckTagBit(tag,kMCElectron),
1434          CheckTagBit(tag,kMCMuon),
1435          CheckTagBit(tag,kMCPion),
1436          CheckTagBit(tag,kMCProton),
1437          CheckTagBit(tag,kMCAntiNeutron),
1438          CheckTagBit(tag,kMCKaon),
1439          CheckTagBit(tag,kMCAntiProton),
1440          CheckTagBit(tag,kMCAntiNeutron),
1441          CheckTagBit(tag,kMCOther),
1442          CheckTagBit(tag,kMCUnknown),
1443          CheckTagBit(tag,kMCBadLabel)
1444          );
1445
1446 }
1447
1448