PWG4/PartCorrBase/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 /* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
  16
  17 //_________________________________________________________________________
  18 // Class for analysis utils for MC data
  19 // stored in stack or event header.
  20 // Contains:
  21 //  - method to check the origin of a given track/cluster
  22 //  - method to obtain the generated jets
  23 //
  24 //*-- Author: Gustavo Conesa (LNF-INFN)
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27
  28 // --- ROOT system ---
  29 #include <TMath.h>
  30 #include <TList.h>
  31 #include "TParticle.h"
  32 #include "TDatabasePDG.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(), fCurrentEvent(-1), fDebug(-1),
  46 fJetsList(new TList), fMCGenerator("PYTHIA")
  47 {
  48   //Ctor
  49 }
  50 /*
  51  //____________________________________________________________________________
  52  AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
  53  TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
  54  fJetsList(new TList), fMCGenerator(mcutils.fMCGenerator)
  55  {
  56  // cpy ctor
  57
  58  }
  59  */
  60
  61 //_________________________________________________________________________
  62 //AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
  63 //{
  64 //  // assignment operator
  65 //
  66 //  if(&mcutils == this) return *this;
  67 //  fCurrentEvent = mcutils.fCurrentEvent ;
  68 //  fDebug        = mcutils.fDebug;
  69 //  fJetsList     = mcutils.fJetsList;
  70 //  fMCGenerator  = mcutils.fMCGenerator;
  71 //
  72 //  return *this;
  73 //}
  74
  75 //____________________________________________________________________________
  76 AliMCAnalysisUtils::~AliMCAnalysisUtils()
  77 {
  78   // Remove all pointers.
  79
  80   if (fJetsList) {
  81     fJetsList->Clear();
  82     delete fJetsList ;
  83   }
  84 }
  85
  86
  87 //_________________________________________________________________________
  88 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, const Int_t nlabels, AliCaloTrackReader* reader, const Int_t input = 0) {
  89         //Play with the montecarlo particles if available
  90         Int_t tag = 0;
  91
  92   if(nlabels<=0) {
  93     printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
  94     return kMCBadLabel;
  95   }
  96
  97         //Select where the information is, ESD-galice stack or AOD mcparticles branch
  98         if(reader->ReadStack()){
  99                 tag = CheckOriginInStack(label, nlabels, reader->GetStack());
 100         }
 101         else if(reader->ReadAODMCParticles()){
 102                 tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(input));
 103         }
 104
 105         return tag ;
 106 }
 107
 108 //_________________________________________________________________________
 109 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) {
 110         //Play with the montecarlo particles if available
 111         Int_t tag = 0;
 112
 113   if(label<0) {
 114     printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
 115     return kMCBadLabel;
 116   }
 117
 118         Int_t labels[]={label};
 119
 120         //Select where the information is, ESD-galice stack or AOD mcparticles branch
 121         if(reader->ReadStack()){
 122                 tag = CheckOriginInStack(labels, 1,reader->GetStack());
 123         }
 124         else if(reader->ReadAODMCParticles()){
 125                 tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(input));
 126         }
 127
 128         return tag ;
 129 }
 130
 131 //_________________________________________________________________________
 132 Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, const Int_t nlabels, AliStack* stack) {
 133   // Play with the MC stack if available. Tag particles depending on their origin.
 134   // Do same things as in CheckOriginInAOD but different input.
 135
 136   //generally speaking, label is the MC label of a reconstructed
 137   //entity (track, cluster, etc) for which we want to know something
 138   //about its heritage, but one can also use it directly with stack
 139   //particles not connected to reconstructed entities
 140
 141   if(!stack) {
 142     if (fDebug >=0)
 143       printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
 144     return -1;
 145   }
 146
 147   Int_t tag = 0;
 148   Int_t label=labels[0];//Most significant particle contributing to the cluster
 149
 150   if(label >= 0 && label < stack->GetNtrack()){
 151     //MC particle of interest is the "mom" of the entity
 152     TParticle * mom = stack->Particle(label);
 153     Int_t iMom = label;
 154     Int_t mPdg = TMath::Abs(mom->GetPdgCode());
 155     Int_t mStatus =  mom->GetStatusCode() ;
 156     Int_t iParent =  mom->GetFirstMother() ;
 157     if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 158
 159     //GrandParent of the entity
 160     TParticle * parent = NULL;
 161     Int_t pPdg = -1;
 162     Int_t pStatus =-1;
 163     if(iParent >= 0){
 164       parent = stack->Particle(iParent);
 165       if(parent){
 166         pPdg = TMath::Abs(parent->GetPdgCode());
 167         pStatus = parent->GetStatusCode();
 168       }
 169     }
 170     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
 171
 172     if(fDebug > 2 ) {
 173                   printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
 174                   printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 175                   printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 176     }
 177
 178     //Check if "mother" of entity is converted, if not, get the first non converted mother
 179     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
 180       SetTagBit(tag,kMCConversion);
 181       //Check if the mother is photon or electron with status not stable
 182       while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
 183         //Mother
 184         iMom = mom->GetFirstMother();
 185         mom = stack->Particle(iMom);
 186         mPdg = TMath::Abs(mom->GetPdgCode());
 187         mStatus =  mom->GetStatusCode() ;
 188         iParent =  mom->GetFirstMother() ;
 189         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 190
 191         //GrandParent
 192         if(iParent >= 0){
 193           parent = stack->Particle(iParent);
 194           if(parent){
 195             pPdg = TMath::Abs(parent->GetPdgCode());
 196             pStatus = parent->GetStatusCode();
 197           }
 198         }
 199         else {// in case of gun/box simulations
 200           pPdg    = 0;
 201           pStatus = 0;
 202           break;
 203         }
 204       }//while
 205       if(fDebug > 2 ) {
 206         printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
 207         printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 208         printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
 209       }
 210
 211     }//mother and parent are electron or photon and have status 0
 212     else if((mPdg == 22 || mPdg == 11) && mStatus == 0){
 213       //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
 214       if(pPdg == 2112 ||  pPdg == 211  ||  pPdg == 321 ||
 215          pPdg == 2212 ||  pPdg == 130  ||  pPdg == 13 ) {
 216         SetTagBit(tag,kMCConversion);
 217         iMom = mom->GetFirstMother();
 218         mom  = stack->Particle(iMom);
 219         mPdg = TMath::Abs(mom->GetPdgCode());
 220
 221         if(fDebug > 2 ) {
 222           printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
 223           printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
 224         }
 225       }//hadron converted
 226
 227       //Comment for the next lines, we do not check the parent of the hadron for the moment.
 228       //iParent =  mom->GetFirstMother() ;
 229       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
 230
 231       //GrandParent
 232       //if(iParent >= 0){
 233       //        parent = stack->Particle(iParent);
 234       //        pPdg = TMath::Abs(parent->GetPdgCode());
 235       //}
 236     }
 237     // conversion into electrons/photons checked
 238
 239     //first check for typical charged particles
 240     if(mPdg == 13) SetTagBit(tag,kMCMuon);
 241     else if(mPdg == 211) SetTagBit(tag,kMCPion);
 242     else if(mPdg == 321) SetTagBit(tag,kMCKaon);
 243     else if(mPdg == 2212) SetTagBit(tag,kMCProton);
 244     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 245     else if(mPdg == 111)  {
 246       SetTagBit(tag,kMCPi0Decay);
 247       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
 248       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 249     }
 250     else if(mPdg == 221) {
 251       SetTagBit(tag,kMCEtaDecay);
 252       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
 253       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
 254     }
 255     //Photons
 256     else if(mPdg == 22){
 257       SetTagBit(tag,kMCPhoton);
 258       if(mStatus == 1){ //undecayed particle
 259         if(fMCGenerator == "PYTHIA"){
 260           if(iParent < 8 && iParent > 5) {//outgoing partons
 261             if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 262             else SetTagBit(tag,kMCFragmentation);
 263           }//Outgoing partons
 264           else  if(iParent <= 5) {
 265             SetTagBit(tag, kMCISR); //Initial state radiation
 266           }
 267           else if(pStatus == 11){//Decay
 268             if(pPdg == 111) {
 269               SetTagBit(tag,kMCPi0Decay);
 270               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 271               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 272             }
 273             else if (pPdg == 221) {
 274               SetTagBit(tag, kMCEtaDecay);
 275               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 276               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 277             }
 278             else SetTagBit(tag,kMCOtherDecay);
 279           }//Decay
 280           else {
 281             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",
 282                                    mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
 283             if(pPdg == 111) {
 284               SetTagBit(tag,kMCPi0Decay);
 285               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
 286               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 287             }
 288             else if (pPdg == 221) {
 289               SetTagBit(tag, kMCEtaDecay);
 290               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
 291               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
 292             }
 293             else SetTagBit(tag,kMCOtherDecay);
 294           }
 295         }//PYTHIA
 296
 297         else if(fMCGenerator == "HERWIG"){
 298           if(pStatus < 197){//Not decay
 299             while(1){
 300               if(parent){
 301                 if(parent->GetFirstMother()<=5) break;
 302                 iParent = parent->GetFirstMother();
 303                 parent=stack->Particle(iParent);
 304                 pStatus= parent->GetStatusCode();
 305                 pPdg = TMath::Abs(parent->GetPdgCode());
 306               } else break;
 307             }//Look for the parton
 308
 309             if(iParent < 8 && iParent > 5) {
 310               if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 311               else SetTagBit(tag,kMCFragmentation);
 312             }
 313             else SetTagBit(tag,kMCISR);//Initial state radiation
 314           }//Not decay
 315           else{//Decay
 316             if(pPdg == 111) {
 317               SetTagBit(tag,kMCPi0Decay);
 318               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
 319               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 320             }
 321             else if (pPdg == 221) {
 322               SetTagBit(tag,kMCEtaDecay);
 323               if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
 324               CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
 325             }
 326             else SetTagBit(tag,kMCOtherDecay);
 327           }//Decay
 328         }//HERWIG
 329
 330         else SetTagBit(tag,kMCUnknown);
 331
 332       }//Status 1 : created by event generator
 333
 334       else if(mStatus == 0){ // geant
 335         if(pPdg == 111) {
 336           SetTagBit(tag,kMCPi0Decay);
 337           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
 338           CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
 339         }
 340         else if (pPdg == 221) {
 341           SetTagBit(tag,kMCEtaDecay);
 342           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
 343           CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
 344         }
 345         else  SetTagBit(tag,kMCOtherDecay);
 346       }//status 0 : geant generated
 347
 348     }//Mother Photon
 349
 350     //Electron check.  Where did that electron come from?
 351     else if(mPdg == 11){ //electron
 352       if(pPdg == 11 && parent){
 353         Int_t iGrandma = parent->GetFirstMother();
 354         if(iGrandma >= 0) {
 355           TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
 356           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 357
 358           if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 359           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 360         }
 361       }
 362       SetTagBit(tag,kMCElectron);
 363       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
 364       if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
 365       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
 366       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
 367       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
 368         if(parent){
 369           Int_t iGrandma = parent->GetFirstMother();
 370           if(iGrandma >= 0) {
 371             TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
 372             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 373             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
 374             else SetTagBit(tag,kMCEFromC); //c-->e
 375           } else SetTagBit(tag,kMCEFromC); //c-->e
 376         }//parent
 377       } else {
 378         //if it is not from any of the above, where is it from?
 379         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 380         else SetTagBit(tag,kMCOtherDecay);
 381         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);
 382       }
 383     }//electron check
 384     //Cluster was made by something else
 385     else {
 386       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
 387       SetTagBit(tag,kMCUnknown);
 388     }
 389   }//Good label value
 390   else{// Bad label
 391
 392     if(label < 0 && (fDebug >= 0))
 393       printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***:  label %d \n", label);
 394     if(label >=  stack->GetNtrack() &&  (fDebug >= 0))
 395       printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 396     SetTagBit(tag,kMCUnknown);
 397   }//Bad label
 398
 399   return tag;
 400
 401 }
 402
 403
 404 //_________________________________________________________________________
 405 Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, const Int_t nlabels, TClonesArray *mcparticles) {
 406   // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
 407   // Do same things as in CheckOriginInStack but different input.
 408   if(!mcparticles) {
 409     if(fDebug >= 0)
 410       printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
 411     return -1;
 412   }
 413
 414   Int_t tag = 0;
 415   Int_t label=labels[0];//Most significant particle contributing to the cluster
 416
 417   Int_t nprimaries = mcparticles->GetEntriesFast();
 418   if(label >= 0 && label < nprimaries){
 419     //Mother
 420     AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
 421     Int_t iMom = label;
 422     Int_t mPdg = TMath::Abs(mom->GetPdgCode());
 423     Int_t iParent =  mom->GetMother() ;
 424     if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 425
 426     //GrandParent
 427     AliAODMCParticle * parent = NULL ;
 428     Int_t pPdg = -1;
 429     if(iParent >= 0){
 430       parent = (AliAODMCParticle *) mcparticles->At(iParent);
 431       pPdg = TMath::Abs(parent->GetPdgCode());
 432     }
 433     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
 434
 435     if(fDebug > 2 ) {
 436                   printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
 437                   printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 438                   if(parent)
 439         printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 440     }
 441
 442     //Check if mother is converted, if not, get the first non converted mother
 443     if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
 444       SetTagBit(tag,kMCConversion);
 445       //Check if the mother is photon or electron with status not stable
 446       while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
 447         //Mother
 448         iMom = mom->GetMother();
 449         mom = (AliAODMCParticle *) mcparticles->At(iMom);
 450         mPdg = TMath::Abs(mom->GetPdgCode());
 451         iParent =  mom->GetMother() ;
 452         if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 453
 454         //GrandParent
 455         if(iParent >= 0 && parent){
 456           parent = (AliAODMCParticle *) mcparticles->At(iParent);
 457           pPdg = TMath::Abs(parent->GetPdgCode());
 458         }
 459         // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 460         // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 461
 462       }//while
 463
 464       if(fDebug > 2 ) {
 465         printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
 466         printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 467         if(parent)
 468           printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 469       }
 470
 471     }//mother and parent are electron or photon and have status 0 and parent is photon or electron
 472     else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){
 473       //Still a conversion but only one electron/photon generated. Just from hadrons
 474       if(pPdg == 2112 ||  pPdg == 211 ||  pPdg == 321 ||
 475          pPdg == 2212 ||  pPdg == 130 ||  pPdg == 13 ) {
 476         SetTagBit(tag,kMCConversion);
 477         iMom = mom->GetMother();
 478         mom = (AliAODMCParticle *) mcparticles->At(iMom);
 479         mPdg = TMath::Abs(mom->GetPdgCode());
 480
 481         if(fDebug > 2 ) {
 482           printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
 483           printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
 484         }
 485       }//hadron converted
 486
 487       //Comment for next lines, we do not check the parent of the hadron for the moment.
 488       //iParent =  mom->GetMother() ;
 489       //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
 490
 491       //GrandParent
 492       //if(iParent >= 0){
 493       //        parent = (AliAODMCParticle *) mcparticles->At(iParent);
 494       //        pPdg = TMath::Abs(parent->GetPdgCode());
 495       //}
 496     }
 497
 498     // conversion into electrons/photons checked
 499
 500     //first check for typical charged particles
 501     if(mPdg == 13) SetTagBit(tag,kMCMuon);
 502     else if(mPdg == 211) SetTagBit(tag,kMCPion);
 503     else if(mPdg == 321) SetTagBit(tag,kMCKaon);
 504     else if(mPdg == 2212) SetTagBit(tag,kMCProton);
 505     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
 506     else if(mPdg == 111)  {
 507       SetTagBit(tag,kMCPi0Decay);
 508       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
 509       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 510     }
 511     else if(mPdg == 221)  {
 512       SetTagBit(tag,kMCEtaDecay);
 513       if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
 514       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 515     }
 516     //Photons
 517     else if(mPdg == 22){
 518       SetTagBit(tag,kMCPhoton);
 519       if(mom->IsPhysicalPrimary()){ //undecayed particle
 520         if(iParent < 8 && iParent > 5) {//outgoing partons
 521           if(pPdg == 22) SetTagBit(tag,kMCPrompt);
 522           else SetTagBit(tag,kMCFragmentation);
 523         }//Outgoing partons
 524         else if(iParent <= 5) {
 525           SetTagBit(tag, kMCISR); //Initial state radiation
 526         }
 527         else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
 528           if(pPdg == 111){
 529             SetTagBit(tag,kMCPi0Decay);
 530             if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
 531             CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 532           }
 533           else if (pPdg == 221) {
 534             SetTagBit(tag, kMCEtaDecay);
 535             if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
 536             CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 537           }
 538           else SetTagBit(tag,kMCOtherDecay);
 539         }//Decay
 540         else {
 541           if(parent)printf("AliMCAnalysisUtils::ChecOrigingInAOD() - 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",
 542                  mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
 543           SetTagBit(tag,kMCOtherDecay);//Check
 544         }
 545       }// Pythia generated
 546       else if(!mom->IsPrimary()){       //Decays
 547         if(pPdg == 111){
 548           SetTagBit(tag,kMCPi0Decay);
 549           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
 550           CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
 551         }
 552         else if (pPdg == 221) {
 553           SetTagBit(tag,kMCEtaDecay);
 554           if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n");
 555           CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
 556         }
 557         else  SetTagBit(tag,kMCOtherDecay);
 558       }//not primary : geant generated, decays
 559       else  {
 560         //printf("UNKNOWN 1, mom  pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
 561         //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
 562         SetTagBit(tag,kMCUnknown);
 563       }
 564     }//Mother Photon
 565
 566     //Electron check.  Where did that electron come from?
 567     else if(mPdg == 11){ //electron
 568       if(pPdg == 11 && parent){
 569         Int_t iGrandma = parent->GetMother();
 570         if(iGrandma >= 0) {
 571           AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
 572           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 573
 574           if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
 575           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
 576         }
 577       }
 578       SetTagBit(tag,kMCElectron);
 579       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
 580       if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
 581       else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
 582       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
 583       else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
 584         if(parent){
 585           Int_t iGrandma = parent->GetMother();
 586           if(iGrandma >= 0) {
 587             AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
 588             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
 589             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
 590             else SetTagBit(tag,kMCEFromC); //c-hadron decay
 591           } else SetTagBit(tag,kMCEFromC); //c-hadron decay
 592         }//parent
 593       } else { //prompt or other decay
 594         TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
 595         TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
 596         if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg);
 597         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
 598         else SetTagBit(tag,kMCOtherDecay);
 599       }
 600     }//electron check
 601     //cluster was made by something else
 602     else {
 603       if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg);
 604       SetTagBit(tag,kMCUnknown);
 605     }
 606   }//Good label value
 607   else{//Bad label
 608
 609     if(label < 0 && (fDebug >= 0) )
 610       printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***:  label %d \n", label);
 611     if(label >=  mcparticles->GetEntriesFast() &&  (fDebug >= 0) )
 612       printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***:  label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
 613     SetTagBit(tag,kMCUnknown);
 614
 615   }//Bad label
 616
 617   return tag;
 618
 619 }
 620
 621 //_________________________________________________________________________
 622 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
 623                                                     AliStack *stack, Int_t &tag) {
 624         //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
 625
 626         if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
 627                 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 628                           labels[0],stack->GetNtrack(), nlabels);
 629                 return;
 630         }
 631
 632         TParticle * meson = stack->Particle(mesonIndex);
 633         Int_t mesonPdg    = meson->GetPdgCode();
 634         if(mesonPdg!=111 && mesonPdg!=221){
 635                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 636                 return;
 637         }
 638
 639         if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
 640
 641         //Check if meson decayed into 2 daughters or if both were kept.
 642         if(meson->GetNDaughters() != 2){
 643                 if(fDebug > 2)
 644                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 645                 return;
 646         }
 647
 648         //Get the daughters
 649         Int_t iPhoton0 = meson->GetDaughter(0);
 650         Int_t iPhoton1 = meson->GetDaughter(1);
 651         TParticle *photon0 = stack->Particle(iPhoton0);
 652         TParticle *photon1 = stack->Particle(iPhoton1);
 653
 654         //Check if both daughters are photons
 655         if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
 656                 if(fDebug > 2)
 657                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 658                 return;
 659         }
 660
 661         if(fDebug > 2)
 662                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 663
 664         //Check if both photons contribute to the cluster
 665         Bool_t okPhoton0 = kFALSE;
 666         Bool_t okPhoton1 = kFALSE;
 667
 668         if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
 669
 670         for(Int_t i = 0; i < nlabels; i++){
 671                 if(fDebug > 3) printf("\t  at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
 672
 673                 //If we already found both, break the loop
 674                 if(okPhoton0 && okPhoton1) break;
 675
 676                 Int_t index =   labels[i];
 677                 if      (iPhoton0 == index) {
 678                         okPhoton0 = kTRUE;
 679                         continue;
 680                 }
 681                 else if (iPhoton1 == index) {
 682                         okPhoton1 = kTRUE;
 683                         continue;
 684                 }
 685
 686                 //Trace back the mother in case it was a conversion
 687                 TParticle * daught = stack->Particle(index);
 688                 Int_t tmpindex = daught->GetFirstMother();
 689                 if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
 690                 while(tmpindex>=0){
 691                         //MC particle of interest is the mother
 692                         if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
 693                         daught   = stack->Particle(tmpindex);
 694                         if      (iPhoton0 == tmpindex) {
 695                                 okPhoton0 = kTRUE;
 696                                 break;
 697                         }
 698                         else if (iPhoton1 == tmpindex) {
 699                                 okPhoton1 = kTRUE;
 700                                 break;
 701                         }
 702                         tmpindex = daught->GetFirstMother();
 703                 }//While to check if pi0/eta daughter was one of these contributors to the cluster
 704
 705                 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
 706                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
 707
 708         }//loop on list of labels
 709
 710         //If both photons contribute tag as the corresponding meson.
 711         if(okPhoton0 && okPhoton1){
 712                 if(fDebug > 2)
 713                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
 714
 715                 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 716                 else  SetTagBit(tag,kMCEta);
 717         }
 718
 719 }
 720
 721 //_________________________________________________________________________
 722 void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
 723                                                     TClonesArray *mcparticles, Int_t &tag) {
 724         //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
 725
 726         if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1) {
 727                 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
 728                           labels[0],mcparticles->GetEntriesFast(), nlabels);
 729                 return;
 730         }
 731
 732
 733   AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
 734         Int_t mesonPdg = meson->GetPdgCode();
 735         if(mesonPdg != 111 && mesonPdg != 221) {
 736                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
 737                 return;
 738         }
 739
 740         if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
 741
 742
 743         //Get the daughters
 744         if(meson->GetNDaughters() != 2){
 745                 if(fDebug > 2)
 746                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
 747                 return;
 748         }
 749         Int_t iPhoton0 = meson->GetDaughter(0);
 750         Int_t iPhoton1 = meson->GetDaughter(1);
 751         //if((iPhoton0 == -1) || (iPhoton1 == -1)){
 752         //      if(fDebug > 2)
 753         //              printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overalapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
 754         //      return;
 755         //}
 756         AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
 757         AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
 758
 759         //Check if both daughters are photons
 760         if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22){
 761                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
 762                 return;
 763         }
 764
 765         if(fDebug > 2)
 766                 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
 767
 768         //Check if both photons contribute to the cluster
 769         Bool_t okPhoton0 = kFALSE;
 770         Bool_t okPhoton1 = kFALSE;
 771
 772         if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
 773
 774         for(Int_t i = 0; i < nlabels; i++){
 775                 if(fDebug > 3) printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
 776
 777                 //If we already found both, break the loop
 778                 if(okPhoton0 && okPhoton1) break;
 779
 780                 Int_t index =   labels[i];
 781                 if      (iPhoton0 == index) {
 782                         okPhoton0 = kTRUE;
 783                         continue;
 784                 }
 785                 else if (iPhoton1 == index) {
 786                         okPhoton1 = kTRUE;
 787                         continue;
 788                 }
 789
 790                 //Trace back the mother in case it was a conversion
 791                 AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
 792                 Int_t tmpindex = daught->GetMother();
 793                 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
 794
 795                 while(tmpindex>=0){
 796
 797                         //MC particle of interest is the mother
 798                         if(fDebug > 3) printf("\t parent index %d\n",tmpindex);
 799                         daught   = (AliAODMCParticle*) mcparticles->At(tmpindex);
 800                         //printf("tmpindex %d\n",tmpindex);
 801                         if      (iPhoton0 == tmpindex) {
 802                                 okPhoton0 = kTRUE;
 803                                 break;
 804                         }
 805                         else if (iPhoton1 == tmpindex) {
 806                                 okPhoton1 = kTRUE;
 807                                 break;
 808                         }
 809                         tmpindex = daught->GetMother();
 810                 }//While to check if pi0/eta daughter was one of these contributors to the cluster
 811
 812                 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
 813                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
 814
 815         }//loop on list of labels
 816
 817         //If both photons contribute tag as the corresponding meson.
 818         if(okPhoton0 && okPhoton1){
 819                 if(fDebug > 2)
 820                         printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
 821
 822                 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
 823                 else  SetTagBit(tag,kMCEta);
 824         }
 825
 826 }
 827
 828 //_________________________________________________________________________
 829 TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * const reader){
 830   //Return list of jets (TParticles) and index of most likely parton that originated it.
 831   AliStack * stack = reader->GetStack();
 832   Int_t iEvent = reader->GetEventNumber();
 833   AliGenEventHeader * geh = reader->GetGenEventHeader();
 834   if(fCurrentEvent!=iEvent){
 835     fCurrentEvent = iEvent;
 836     fJetsList = new TList;
 837     Int_t nTriggerJets = 0;
 838     Float_t tmpjet[]={0,0,0,0};
 839
 840     //printf("Event %d %d\n",fCurrentEvent,iEvent);
 841     //Get outgoing partons
 842     if(stack->GetNtrack() < 8) return fJetsList;
 843     TParticle * parton1 =  stack->Particle(6);
 844     TParticle * parton2 =  stack->Particle(7);
 845     if(fDebug > 2){
 846       printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 847              parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
 848       printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 849              parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
 850                 }
 851     //          //Trace the jet from the mother parton
 852     //          Float_t pt  = 0;
 853     //          Float_t pt1 = 0;
 854     //          Float_t pt2 = 0;
 855     //          Float_t e   = 0;
 856     //          Float_t e1  = 0;
 857     //          Float_t e2  = 0;
 858     //          TParticle * tmptmp = new TParticle;
 859     //          for(Int_t i = 0; i< stack->GetNprimary(); i++){
 860     //                  tmptmp = stack->Particle(i);
 861
 862     //                  if(tmptmp->GetStatusCode() == 1){
 863     //                          pt = tmptmp->Pt();
 864     //                          e =  tmptmp->Energy();
 865     //                          Int_t imom = tmptmp->GetFirstMother();
 866     //                          Int_t imom1 = 0;
 867     //                          //printf("1st imom %d\n",imom);
 868     //                          while(imom > 5){
 869     //                                  imom1=imom;
 870     //                                  tmptmp = stack->Particle(imom);
 871     //                                  imom = tmptmp->GetFirstMother();
 872     //                                  //printf("imom %d       \n",imom);
 873     //                          }
 874     //                          //printf("Last imom %d %d\n",imom1, imom);
 875     //                          if(imom1 == 6) {
 876     //                                  pt1+=pt;
 877     //                                  e1+=e;
 878     //                          }
 879     //                          else if (imom1 == 7){
 880     //                                  pt2+=pt;
 881     //                                  e2+=e;                                  }
 882     //                  }// status 1
 883
 884     //          }// for
 885
 886     //          printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
 887
 888                 //Get the jet, different way for different generator
 889                 //PYTHIA
 890     if(fMCGenerator == "PYTHIA"){
 891       TParticle * jet =  0x0;
 892       AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
 893       nTriggerJets =  pygeh->NTriggerJets();
 894       if(fDebug > 1)
 895         printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
 896
 897       Int_t iparton = -1;
 898       for(Int_t i = 0; i< nTriggerJets; i++){
 899         iparton=-1;
 900         pygeh->TriggerJet(i, tmpjet);
 901         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 902         //Assign an outgoing parton as mother
 903         Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
 904         Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
 905         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
 906         else  jet->SetFirstMother(6);
 907         //jet->Print();
 908         if(fDebug > 1)
 909           printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 910                  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
 911         fJetsList->Add(jet);
 912       }
 913     }//Pythia triggered jets
 914     //HERWIG
 915     else if (fMCGenerator=="HERWIG"){
 916       Int_t pdg = -1;
 917       //Check parton 1
 918       TParticle * tmp = parton1;
 919       if(parton1->GetPdgCode()!=22){
 920         while(pdg != 94){
 921           if(tmp->GetFirstDaughter()==-1) return fJetsList;
 922           tmp = stack->Particle(tmp->GetFirstDaughter());
 923           pdg = tmp->GetPdgCode();
 924         }//while
 925
 926         //Add found jet to list
 927         TParticle *jet1 = new TParticle(*tmp);
 928         jet1->SetFirstMother(6);
 929         fJetsList->Add(jet1);
 930         //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
 931         //tmp = stack->Particle(tmp->GetFirstDaughter());
 932         //tmp->Print();
 933         //jet1->Print();
 934         if(fDebug > 1)
 935           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 936                  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
 937       }//not photon
 938
 939       //Check parton 2
 940       pdg = -1;
 941       tmp = parton2;
 942       Int_t i = -1;
 943       if(parton2->GetPdgCode()!=22){
 944         while(pdg != 94){
 945           if(tmp->GetFirstDaughter()==-1) return fJetsList;
 946           i = tmp->GetFirstDaughter();
 947           tmp = stack->Particle(tmp->GetFirstDaughter());
 948           pdg = tmp->GetPdgCode();
 949         }//while
 950         //Add found jet to list
 951         TParticle *jet2 = new TParticle(*tmp);
 952         jet2->SetFirstMother(7);
 953         fJetsList->Add(jet2);
 954         //jet2->Print();
 955         if(fDebug > 1)
 956           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 957                  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
 958         //Int_t first =  tmp->GetFirstDaughter();
 959         //Int_t last  =  tmp->GetLastDaughter();
 960         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
 961                                 //      for(Int_t d = first ; d < last+1; d++){
 962         //                                              tmp = stack->Particle(d);
 963         //                                              if(i == tmp->GetFirstMother())
 964         //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 965         //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
 966         //                         }
 967         //tmp->Print();
 968       }//not photon
 969     }//Herwig generated jets
 970   }
 971
 972   return fJetsList;
 973 }
 974
 975
 976 //________________________________________________________________
 977 void AliMCAnalysisUtils::Print(const Option_t * opt) const
 978 {
 979   //Print some relevant parameters set for the analysis
 980
 981   if(! opt)
 982     return;
 983
 984   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 985
 986   printf("Debug level    = %d\n",fDebug);
 987   printf("MC Generator   = %s\n",fMCGenerator.Data());
 988   printf(" \n");
 989
 990 }
 991
 992