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
  32 //---- ANALYSIS system ----
  33 #include "AliMCAnalysisUtils.h"
  34 #include "AliStack.h"
  35 #include "TParticle.h"
  36 #include "AliGenPythiaEventHeader.h"
  37
  38   ClassImp(AliMCAnalysisUtils)
  39
  40  //________________________________________________
  41   AliMCAnalysisUtils::AliMCAnalysisUtils() :
  42     TObject(), fCurrentEvent(-1), fDebug(-1),
  43     fJetsList(new TList), fMCGenerator("PYTHIA")
  44 {
  45   //Ctor
  46 }
  47
  48 //____________________________________________________________________________
  49 AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
  50   TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
  51   fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
  52 {
  53   // cpy ctor
  54
  55 }
  56
  57 //_________________________________________________________________________
  58 AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
  59 {
  60   // assignment operator
  61
  62   if(&mcutils == this) return *this;
  63   fCurrentEvent = mcutils.fCurrentEvent ;
  64   fDebug        = mcutils.fDebug;
  65   fJetsList     = mcutils.fJetsList;
  66   fMCGenerator  = mcutils.fMCGenerator;
  67
  68   return *this;
  69 }
  70
  71 //____________________________________________________________________________
  72 AliMCAnalysisUtils::~AliMCAnalysisUtils()
  73 {
  74   // Remove all pointers.
  75
  76   if (fJetsList) {
  77     fJetsList->Clear();
  78     delete fJetsList ;
  79   }
  80 }
  81
  82 //_________________________________________________________________________
  83 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliStack * stack) const {
  84   //Play with the MC stack if available
  85   //Check origin of the candidates, good for PYTHIA
  86
  87   if(!stack) {
  88     printf("AliMCAnalysisUtils::CheckOrigin() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
  89     abort();
  90   }
  91   //  printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary());
  92   //   for(Int_t i = 0; i< stack->GetNprimary(); i++){
  93   //      TParticle *particle =   stack->Particle(i);
  94   //                    //particle->Print();
  95   //   }
  96   if(label >= 0 && label <  stack->GetNtrack()){
  97     //Mother
  98     TParticle * mom = stack->Particle(label);
  99     Int_t mPdg = TMath::Abs(mom->GetPdgCode());
 100     Int_t mStatus =  mom->GetStatusCode() ;
 101     Int_t iParent =  mom->GetFirstMother() ;
 102     if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOrigin: Mother is parton %d\n",iParent);
 103
 104     //GrandParent
 105     TParticle * parent = new TParticle ;
 106     Int_t pPdg = -1;
 107     Int_t pStatus =-1;
 108     if(iParent > 0){
 109       parent = stack->Particle(iParent);
 110       pPdg = TMath::Abs(parent->GetPdgCode());
 111       pStatus = parent->GetStatusCode();
 112     }
 113     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent);
 114
 115     //return tag
 116     if(mPdg == 22){
 117       if(mStatus == 1){
 118         if(fMCGenerator == "PYTHIA"){
 119           if(iParent < 8 && iParent > 5) {//outgoing partons
 120             if(pPdg == 22) return kMCPrompt;
 121             else  return kMCFragmentation;
 122           }//Outgoing partons
 123           else if(pStatus == 11){//Decay
 124             if(pPdg == 111) return kMCPi0Decay ;
 125             else if (pPdg == 221)  return kMCEtaDecay ;
 126             else  return kMCOtherDecay ;
 127           }//Decay
 128           else return kMCISR; //Initial state radiation
 129         }//PYTHIA
 130
 131         else if(fMCGenerator == "HERWIG"){
 132           if(pStatus < 197){//Not decay
 133             while(1){
 134               if(parent->GetFirstMother()<=5) break;
 135               iParent = parent->GetFirstMother();
 136               parent=stack->Particle(iParent);
 137               pStatus= parent->GetStatusCode();
 138               pPdg = parent->GetPdgCode();
 139             }//Look for the parton
 140
 141             if(iParent < 8 && iParent > 5) {
 142               if(pPdg == 22) return kMCPrompt;
 143               else  return kMCFragmentation;
 144             }
 145             return kMCISR;//Initial state radiation
 146           }//Not decay
 147           else{//Decay
 148             if(pPdg == 111) return kMCPi0Decay ;
 149             else if (pPdg == 221)  return kMCEtaDecay ;
 150             else  return kMCOtherDecay ;
 151           }//Decay
 152         }//HERWIG
 153         else return  kMCUnknown;
 154       }//Status 1 : Pythia generated
 155       else if(mStatus == 0){
 156         if(pPdg ==22 || pPdg ==11|| pPdg == 2112 ||  pPdg == 211 ||
 157            pPdg == 321 ||  pPdg == 2212  ||  pPdg == 130  ||  pPdg == 13 )
 158           return kMCConversion ;
 159         if(pPdg == 111) return kMCPi0Decay ;
 160         else if (pPdg == 221)  return kMCEtaDecay ;
 161         else  return kMCOtherDecay ;
 162       }//status 0 : geant generated
 163     }//Mother Photon
 164     else if(mPdg == 111)  return kMCPi0 ;
 165     else if(mPdg == 221)  return kMCEta ;
 166     else if(mPdg ==11){
 167 //       if(pPdg !=22 &&mStatus == 0) {
 168 //      printf("Origin electron, pT %f, status %d, parent %s, pstatus %d, vx %f, vy %f, vz %f\n",
 169 //             mom->Pt(),mStatus, parent->GetName(),pStatus,mom->Vx(),mom->Vy(), mom->Vz());
 170
 171 //       }
 172
 173       if(mStatus == 0) {
 174         if(pPdg ==22) return kMCConversion ;
 175         else if (pStatus == 0) return kMCConversion;
 176         else return kMCElectron ;
 177       }
 178       else return kMCElectron ;
 179     }
 180     else return kMCUnknown;
 181   }//Good label value
 182   else{
 183     if(label < 0 ) printf("AliMCAnalysisUtils::CheckOrigin: *** bad label or no stack ***:  label %d \n", label);
 184     if(label >=  stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOrigin: *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 185     return kMCUnknown;
 186   }//Bad label
 187
 188   return kMCUnknown;
 189
 190 }
 191
 192 //_________________________________________________________________________
 193 TList * AliMCAnalysisUtils::GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) {
 194  //Return list of jets (TParticles) and index of most likely parton that originated it.
 195
 196   if(fCurrentEvent!=iEvent){
 197     fCurrentEvent = iEvent;
 198     fJetsList = new TList;
 199     Int_t nTriggerJets = 0;
 200     Float_t tmpjet[]={0,0,0,0};
 201
 202     //printf("Event %d %d\n",fCurrentEvent,iEvent);
 203     //Get outgoing partons
 204     if(stack->GetNtrack() < 8) return fJetsList;
 205     TParticle * parton1 =  stack->Particle(6);
 206     TParticle * parton2 =  stack->Particle(7);
 207     if(fDebug > 2){
 208       printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 209              parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
 210       printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 211              parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
 212                 }
 213 //              //Trace the jet from the mother parton
 214 //              Float_t pt  = 0;
 215 //              Float_t pt1 = 0;
 216 //              Float_t pt2 = 0;
 217 //              Float_t e   = 0;
 218 //              Float_t e1  = 0;
 219 //              Float_t e2  = 0;
 220 //              TParticle * tmptmp = new TParticle;
 221 //              for(Int_t i = 0; i< stack->GetNprimary(); i++){
 222 //                      tmptmp = stack->Particle(i);
 223
 224 //                      if(tmptmp->GetStatusCode() == 1){
 225 //                              pt = tmptmp->Pt();
 226 //                              e =  tmptmp->Energy();
 227 //                              Int_t imom = tmptmp->GetFirstMother();
 228 //                              Int_t imom1 = 0;
 229 //                              //printf("1st imom %d\n",imom);
 230 //                              while(imom > 5){
 231 //                                      imom1=imom;
 232 //                                      tmptmp = stack->Particle(imom);
 233 //                                      imom = tmptmp->GetFirstMother();
 234 //                                      //printf("imom %d       \n",imom);
 235 //                              }
 236 //                              //printf("Last imom %d %d\n",imom1, imom);
 237 //                              if(imom1 == 6) {
 238 //                                      pt1+=pt;
 239 //                                      e1+=e;
 240 //                              }
 241 //                              else if (imom1 == 7){
 242 //                                      pt2+=pt;
 243 //                                      e2+=e;                                  }
 244 //                      }// status 1
 245
 246 //              }// for
 247
 248 //              printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
 249
 250                 //Get the jet, different way for different generator
 251                 //PYTHIA
 252     if(fMCGenerator == "PYTHIA"){
 253       TParticle * jet =  new TParticle;
 254       AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
 255       nTriggerJets =  pygeh->NTriggerJets();
 256       if(fDebug > 1)
 257          printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
 258
 259       Int_t iparton = -1;
 260       for(Int_t i = 0; i< nTriggerJets; i++){
 261         iparton=-1;
 262         pygeh->TriggerJet(i, tmpjet);
 263         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 264         //Assign an outgoing parton as mother
 265         Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
 266         Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
 267         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
 268         else  jet->SetFirstMother(6);
 269         //jet->Print();
 270         if(fDebug > 1)
 271           printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 272                  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
 273         fJetsList->Add(jet);
 274       }
 275     }//Pythia triggered jets
 276     //HERWIG
 277     else if (fMCGenerator=="HERWIG"){
 278       Int_t pdg = -1;
 279       //Check parton 1
 280       TParticle * tmp = parton1;
 281       if(parton1->GetPdgCode()!=22){
 282         while(pdg != 94){
 283           if(tmp->GetFirstDaughter()==-1) return fJetsList;
 284           tmp = stack->Particle(tmp->GetFirstDaughter());
 285           pdg = tmp->GetPdgCode();
 286         }//while
 287
 288         //Add found jet to list
 289         TParticle *jet1 = new TParticle(*tmp);
 290         jet1->SetFirstMother(6);
 291         fJetsList->Add(jet1);
 292         //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
 293         //tmp = stack->Particle(tmp->GetFirstDaughter());
 294         //tmp->Print();
 295         //jet1->Print();
 296         if(fDebug > 1)
 297           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 298                  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
 299       }//not photon
 300
 301       //Check parton 2
 302       pdg = -1;
 303       tmp = parton2;
 304       Int_t i = -1;
 305       if(parton2->GetPdgCode()!=22){
 306         while(pdg != 94){
 307           if(tmp->GetFirstDaughter()==-1) return fJetsList;
 308           i = tmp->GetFirstDaughter();
 309           tmp = stack->Particle(tmp->GetFirstDaughter());
 310           pdg = tmp->GetPdgCode();
 311         }//while
 312         //Add found jet to list
 313         TParticle *jet2 = new TParticle(*tmp);
 314         jet2->SetFirstMother(7);
 315         fJetsList->Add(jet2);
 316         //jet2->Print();
 317         if(fDebug > 1)
 318           printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 319                  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
 320         //Int_t first =  tmp->GetFirstDaughter();
 321         //Int_t last  =  tmp->GetLastDaughter();
 322         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
 323                                 //      for(Int_t d = first ; d < last+1; d++){
 324 //                                              tmp = stack->Particle(d);
 325 //                                              if(i == tmp->GetFirstMother())
 326 //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
 327 //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
 328 //                         }
 329                                                    //tmp->Print();
 330       }//not photon
 331     }//Herwig generated jets
 332   }
 333
 334   return fJetsList;
 335 }
 336
 337 //________________________________________________________________
 338 void AliMCAnalysisUtils::Print(const Option_t * opt) const
 339 {
 340   //Print some relevant parameters set for the analysis
 341
 342  if(! opt)
 343    return;
 344
 345  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 346
 347  printf("Debug level    = %d\n",fDebug);
 348  printf("MC Generator   = %s\n",fMCGenerator.Data());
 349
 350  printf(" \n");
 351
 352 }
 353
 354