PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx

   1 // $Id$
   2 //
   3 // Jet embedding task.
   4 //
   5 // Author: S.Aiola, C.Loizides
   6
   7 #include "AliJetEmbeddingFromGenTask.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TFolder.h>
  11 #include <TLorentzVector.h>
  12 #include <TParticle.h>
  13 #include <TParticlePDG.h>
  14 #include <TRandom3.h>
  15 #include <TProfile.h>
  16 #include "AliAnalysisManager.h"
  17 #include "AliEMCALDigit.h"
  18 #include "AliEMCALGeometry.h"
  19 #include "AliEMCALRecPoint.h"
  20 #include "AliGenerator.h"
  21 #include "AliHeader.h"
  22 #include "AliLog.h"
  23 #include "AliPicoTrack.h"
  24 #include "AliRun.h"
  25 #include "AliRunLoader.h"
  26 #include "AliStack.h"
  27 #include "AliVCluster.h"
  28 #include "AliVEvent.h"
  29 #include "AliGenPythiaEventHeader.h"
  30 #include "AliPythiaInfo.h"
  31 #include "AliPythiaRndm.h"
  32 ClassImp(AliJetEmbeddingFromGenTask)
  33
  34 //________________________________________________________________________
  35 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask() :
  36   AliJetModelBaseTask("AliJetEmbeddingFromGenTask"),
  37   fGen(0),
  38   fMassless(kFALSE),
  39   fChargedOnly(kFALSE),
  40   fHistPt(0),
  41   fHistEtaPhi(0),
  42   fHistTrials(0),
  43   fHistXsection(0),
  44   fHistPtHard(0)
  45 {
  46   // Default constructor.
  47   SetSuffix("EmbeddedFromGen");
  48 }
  49
  50 //________________________________________________________________________
  51 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask(const char *name, Bool_t drawqa) :
  52   AliJetModelBaseTask(name,drawqa),
  53   fGen(0),
  54   fMassless(kFALSE),
  55   fChargedOnly(kFALSE),
  56   fHistPt(0),
  57   fHistEtaPhi(0),
  58   fHistTrials(0),
  59   fHistXsection(0),
  60   fHistPtHard(0)
  61 {
  62   // Standard constructor.
  63   SetSuffix("EmbeddedFromGen");
  64 }
  65
  66 //________________________________________________________________________
  67 AliJetEmbeddingFromGenTask::~AliJetEmbeddingFromGenTask()
  68 {
  69   // Destructor
  70 }
  71
  72 //________________________________________________________________________
  73 void AliJetEmbeddingFromGenTask::UserCreateOutputObjects()
  74 {
  75   // Create user output.
  76
  77   if (!fQAhistos)
  78     return;
  79
  80   AliJetModelBaseTask::UserCreateOutputObjects();
  81
  82   fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
  83   fOutput->Add(fHistPt);
  84
  85   fHistEtaPhi = new TH2F("fHistEtapHI","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
  86   fOutput->Add(fHistEtaPhi);
  87
  88   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
  89   fHistTrials->GetYaxis()->SetTitle("trials");
  90   fOutput->Add(fHistTrials);
  91
  92   fHistXsection = new TProfile("fHistXsection", "fHistXsection", 1, 0, 1);
  93   fHistXsection->GetYaxis()->SetTitle("xsection");
  94   fOutput->Add(fHistXsection);
  95
  96   fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", 500, 0., 500.);
  97   fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
  98   fHistPtHard->GetYaxis()->SetTitle("counts");
  99   fOutput->Add(fHistPtHard);
 100
 101   PostData(1, fOutput);
 102 }
 103
 104 //________________________________________________________________________
 105 Bool_t AliJetEmbeddingFromGenTask::ExecOnce()
 106 {
 107   // Exec only once.
 108
 109   if (!gAlice) {
 110     new AliRun("gAlice","The ALICE Off-line Simulation Framework");
 111     delete gRandom;
 112     gRandom = new TRandom3(0);
 113   }
 114
 115   TFolder *folder = new TFolder(GetName(),GetName());
 116   AliRunLoader *rl = new AliRunLoader(folder);
 117   gAlice->SetRunLoader(rl);
 118   rl->MakeHeader();
 119   rl->MakeStack();
 120   AliStack *stack = rl->Stack();
 121   fGen->SetStack(stack);
 122   fGen->Init();
 123
 124   if (!(InputEvent()->FindListObject(fTracksName))) {
 125     fOutTracks = new TClonesArray("AliPicoTrack", 1000);
 126     fOutTracks->SetName(fTracksName);
 127     InputEvent()->AddObject(fOutTracks);
 128     fNTracks = 0;
 129   }
 130
 131   if(!fPythiaInfoName.IsNull()) {
 132     if (!(InputEvent()->FindListObject(fPythiaInfoName))) {
 133       fPythiaInfo = new AliPythiaInfo("PythiaInfo");
 134       fPythiaInfo->SetName(fPythiaInfoName);
 135       InputEvent()->AddObject(fPythiaInfo);
 136     }
 137   }
 138   return kTRUE;
 139 }
 140
 141 //________________________________________________________________________
 142 void AliJetEmbeddingFromGenTask::Run()
 143 {
 144   // Embed particles.
 145
 146   if (fCopyArray)
 147     CopyTracks();
 148   AliPythiaRndm::SetPythiaRandom(new TRandom3());
 149   AliPythiaRndm::GetPythiaRandom()->SetSeed(clock()+gSystem->GetPid());
 150   AliStack *stack = fGen->GetStack();
 151   stack->Reset();
 152   fGen->Generate();
 153   const Int_t nprim = stack->GetNprimary();
 154   // reject if partons are missing from stack for some reason
 155   if(nprim < 8) return;
 156   if(fPythiaInfo) {
 157     TParticle *part6 = stack->Particle(6);
 158     TParticle *part7 = stack->Particle(7);
 159
 160     fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
 161     fPythiaInfo->SetPartonPt6(part6->Pt());
 162     fPythiaInfo->SetPartonEta6(part6->Eta());
 163     fPythiaInfo->SetPartonPhi6(part6->Phi());
 164
 165     fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
 166     fPythiaInfo->SetPartonPt7(part7->Pt());
 167     fPythiaInfo->SetPartonEta7(part7->Eta());
 168     fPythiaInfo->SetPartonPhi7(part7->Phi());
 169   }
 170
 171   for (Int_t i=0;i<nprim;++i) {
 172     if (!stack->IsPhysicalPrimary(i))
 173       continue;
 174     TParticle *part = stack->Particle(i);
 175     TParticlePDG *pdg = part->GetPDG(1);
 176     if (!pdg)
 177       continue;
 178     Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
 179     if (fChargedOnly && c==0)  continue;
 180     Double_t pt = part->Pt();
 181     Double_t eta = part->Eta();
 182     Double_t phi = part->Phi();
 183     if (eta<fEtaMin)
 184       continue;
 185     if (eta>fEtaMax)
 186       continue;
 187     if (phi<fPhiMin)
 188       continue;
 189     if (phi>fPhiMax)
 190       continue;
 191     if (pt<fPtMin)
 192       continue;
 193     if (pt>fPtMax)
 194       continue;
 195     Double_t mass = part->GetMass();
 196     if(fMassless) mass = 0.;
 197     fHistPt->Fill(pt);
 198     fHistEtaPhi->Fill(eta,phi);
 199     AddTrack(pt, eta, phi,0,0,0,0,0,0,c,mass);
 200   }
 201
 202   FillPythiaHistograms();
 203 }
 204
 205 //________________________________________________________________________
 206 void AliJetEmbeddingFromGenTask::FillPythiaHistograms() {
 207   //Get PYTHIA info: pt-hard, x-section, trials
 208
 209   if (!fQAhistos)
 210     return;
 211
 212   AliRunLoader *rl = AliRunLoader::Instance();
 213   AliGenPythiaEventHeader *genPH = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
 214   if(genPH) {
 215     Float_t xsec = genPH->GetXsection();
 216     Int_t trials = genPH->Trials();
 217     Float_t pthard = genPH->GetPtHard();
 218     Float_t ptWeight=genPH->EventWeight();
 219     if(fPythiaInfo) fPythiaInfo->SetPythiaEventWeight(ptWeight);
 220     fHistXsection->Fill(0.5,xsec);
 221     fHistTrials->Fill(0.5,trials);
 222     fHistPtHard->Fill(pthard);
 223   }
 224 }