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