bookkeep pt hard info when generating on-the-fly pythia events
[u/mrichter/AliRoot.git] / 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 "AliStack.h"
28 #include "AliVCluster.h"
29 #include "AliVEvent.h"
30 #include "AliGenPythiaEventHeader.h"
31
32 ClassImp(AliJetEmbeddingFromGenTask)
33
34 //________________________________________________________________________
35 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask() : 
36   AliJetModelBaseTask("AliJetEmbeddingFromGenTask"),
37   fGen(0),
38   fHistTrials(0),
39   fHistXsection(0),
40   fHistPtHard(0)
41 {
42   // Default constructor.
43   SetSuffix("EmbeddedFromGen");
44 }
45
46 //________________________________________________________________________
47 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask(const char *name, Bool_t drawqa) :
48   AliJetModelBaseTask(name,drawqa),
49   fGen(0),
50   fHistTrials(0),
51   fHistXsection(0),
52   fHistPtHard(0)
53 {
54   // Standard constructor.
55   SetSuffix("EmbeddedFromGen");
56 }
57
58 //________________________________________________________________________
59 AliJetEmbeddingFromGenTask::~AliJetEmbeddingFromGenTask()
60 {
61   // Destructor
62 }
63
64 //________________________________________________________________________
65 void AliJetEmbeddingFromGenTask::UserCreateOutputObjects()
66 {
67   // Create user output.
68
69   if (!fQAhistos)
70     return;
71
72   AliJetModelBaseTask::UserCreateOutputObjects();
73
74   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
75   fHistTrials->GetYaxis()->SetTitle("trials");
76   fOutput->Add(fHistTrials);
77
78   fHistXsection = new TProfile("fHistXsection", "fHistXsection", 1, 0, 1);
79   fHistXsection->GetYaxis()->SetTitle("xsection");
80   fOutput->Add(fHistXsection);
81
82   fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", 500, 0., 500.);
83   fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
84   fHistPtHard->GetYaxis()->SetTitle("counts");
85   fOutput->Add(fHistPtHard);
86
87   PostData(1, fOutput);
88 }
89
90 //________________________________________________________________________
91 Bool_t AliJetEmbeddingFromGenTask::ExecOnce() 
92 {
93   // Exec only once.
94
95   if (!gAlice) {
96     new AliRun("gAlice","The ALICE Off-line Simulation Framework");
97     delete gRandom;
98     gRandom = new TRandom3(0);
99   }
100
101   TFolder *folder = new TFolder(GetName(),GetName());
102   AliRunLoader *rl = new AliRunLoader(folder);
103   rl->MakeHeader();
104   rl->MakeStack();
105   AliStack *stack = rl->Stack();
106   fGen->SetStack(stack);
107   fGen->Init();
108
109   if (!(InputEvent()->FindListObject(fTracksName))) {
110     fOutTracks = new TClonesArray("AliPicoTrack", 1000);
111     fOutTracks->SetName(fTracksName);
112     InputEvent()->AddObject(fOutTracks);
113     fNTracks = 0;
114   }
115   
116   return kTRUE;
117 }
118
119 //________________________________________________________________________
120 void AliJetEmbeddingFromGenTask::Run() 
121 {
122   // Embed particles.
123
124   if (fCopyArray) 
125     CopyTracks();
126
127   AliStack *stack = fGen->GetStack();
128   stack->Reset();
129   fGen->Generate();
130   const Int_t nprim = stack->GetNprimary();
131   for (Int_t i=0;i<nprim;++i) {
132     if (!stack->IsPhysicalPrimary(i))
133       continue;
134     TParticle *part = stack->Particle(i);
135     TParticlePDG *pdg = part->GetPDG(1);
136     if (!pdg) 
137       continue;
138     Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
139     if (c==0) 
140       continue;
141     Double_t pt = part->Pt();
142     Double_t eta = part->Eta();
143     Double_t phi = part->Phi();
144     if (eta<fEtaMin)
145       continue;
146     if (eta>fEtaMax)
147       continue;
148     if (phi<fPhiMin)
149       continue;
150     if (phi>fPhiMax)
151       continue;
152     if (pt<fPtMin)
153       continue;
154     if (pt>fPtMax)
155       continue;
156     AddTrack(pt, eta, phi,0,0,0,0,0,0,part->GetMass());
157   }
158
159   FillPythiaHistograms();
160 }
161
162 //________________________________________________________________________
163 void AliJetEmbeddingFromGenTask::FillPythiaHistograms() {
164   //Get PYTHIA info: pt-hard, x-section, trials
165
166   if (!fQAhistos)
167     return;
168
169   AliRunLoader *rl = AliRunLoader::Instance();
170   AliGenPythiaEventHeader *genPH = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
171   if(genPH) {
172     Printf("found pythia event header. pThard: %f Trials: %d xsec: %f",genPH->GetPtHard(),genPH->Trials(),genPH->GetXsection());
173
174     Float_t xsec = genPH->GetXsection();
175     Int_t trials = genPH->Trials();
176     Float_t pthard = genPH->GetPtHard();
177
178     fHistXsection->Fill(0.5,xsec);
179     fHistTrials->Fill(0.5,trials);
180     fHistPtHard->Fill(pthard);
181   }
182 }