AliPhysics  9c7580a (9c7580a)
AliJetEmbeddingFromGenTask.cxx
Go to the documentation of this file.
1 //
2 // Jet embedding task.
3 //
4 // Author: S.Aiola, C.Loizides
5 
7 
8 #include <TClonesArray.h>
9 #include <TFolder.h>
10 #include <TLorentzVector.h>
11 #include <TParticle.h>
12 #include <TParticlePDG.h>
13 #include <TRandom3.h>
14 #include <TProfile.h>
15 #include "AliAnalysisManager.h"
16 #include "AliEMCALDigit.h"
17 #include "AliEMCALGeometry.h"
18 #include "AliEMCALRecPoint.h"
19 #include "AliGenerator.h"
20 #include "AliHeader.h"
21 #include "AliLog.h"
22 #include "AliPicoTrack.h"
23 #include "AliRun.h"
24 #include "AliRunLoader.h"
25 #include "AliStack.h"
26 #include "AliVCluster.h"
27 #include "AliVEvent.h"
28 #include "AliGenPythiaEventHeader.h"
29 #include "AliGenHerwigEventHeader.h"
30 #include "AliEmcalPythiaInfo.h"
31 #include "AliPythiaRndm.h"
32 #include "AliHerwigRndm.h"
34 
35 //________________________________________________________________________
38  fGen(0),
39  fMassless(kFALSE),
40  fChargedOnly(kFALSE),
41  fToyModelFragmentation(kFALSE),
42  fToyModelFraction(0.9),
43  fHistPt(0),
44  fHistEtaPhi(0),
45  fHistTrials(0),
46  fHistXsection(0),
47  fHistPtHard(0)
48 {
49  // Default constructor.
50  SetSuffix("EmbeddedFromGen");
51 }
52 
53 //________________________________________________________________________
55  AliJetModelBaseTask(name,drawqa),
56  fGen(0),
57  fMassless(kFALSE),
58  fChargedOnly(kFALSE),
59  fToyModelFragmentation(kFALSE),
60  fToyModelFraction(0.9),
61  fHistPt(0),
62  fHistEtaPhi(0),
63  fHistTrials(0),
64  fHistXsection(0),
65  fHistPtHard(0)
66 {
67  // Standard constructor.
68  SetSuffix("EmbeddedFromGen");
69 }
70 
71 //________________________________________________________________________
73 {
74  // Destructor
75 }
76 
77 //________________________________________________________________________
79 {
80  // Create user output.
81 
82  if (!fQAhistos)
83  return;
84 
86 
87  fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
88  fOutput->Add(fHistPt);
89 
90  fHistEtaPhi = new TH2F("fHistEtapHI","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
91  fOutput->Add(fHistEtaPhi);
92 
93  fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
94  fHistTrials->GetYaxis()->SetTitle("trials");
95  fOutput->Add(fHistTrials);
96 
97  fHistXsection = new TProfile("fHistXsection", "fHistXsection", 1, 0, 1);
98  fHistXsection->GetYaxis()->SetTitle("xsection");
99  fOutput->Add(fHistXsection);
100 
101  fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", 500, 0., 500.);
102  fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
103  fHistPtHard->GetYaxis()->SetTitle("counts");
104  fOutput->Add(fHistPtHard);
105 
106  PostData(1, fOutput);
107 }
108 
109 //________________________________________________________________________
111 {
112  // Exec only once.
113 
114  if (!gAlice) {
115  new AliRun("gAlice","The ALICE Off-line Simulation Framework");
116  delete gRandom;
117  gRandom = new TRandom3(0);
118  }
119 
120  TFolder *folder = new TFolder(GetName(),GetName());
121  AliRunLoader *rl = new AliRunLoader(folder);
122  gAlice->SetRunLoader(rl);
123  rl->MakeHeader();
124  rl->MakeStack();
125  AliStack *stack = rl->Stack();
126  fGen->SetStack(stack);
127  fGen->Init();
128 
129  if (!(InputEvent()->FindListObject(fTracksName))) {
130  fOutTracks = new TClonesArray("AliPicoTrack", 1000);
131  fOutTracks->SetName(fTracksName);
132  InputEvent()->AddObject(fOutTracks);
133  fNTracks = 0;
134  }
135 
136  if(!fPythiaInfoName.IsNull()) {
137  if (!(InputEvent()->FindListObject(fPythiaInfoName))) {
138  fPythiaInfo = new AliEmcalPythiaInfo("PythiaInfo");
139  fPythiaInfo->SetName(fPythiaInfoName);
140  InputEvent()->AddObject(fPythiaInfo);
141  }
142  }
143  return kTRUE;
144 }
145 
146 //________________________________________________________________________
148 {
149  // Embed particles.
150 
151  if (fCopyArray)
152  CopyTracks();
153  if(fGenType<3){
154  AliPythiaRndm::SetPythiaRandom(new TRandom3());
155  AliPythiaRndm::GetPythiaRandom()->SetSeed(clock()+gSystem->GetPid());}
156  if(fGenType==3){
157  AliHerwigRndm::SetHerwigRandom(new TRandom3());
158  AliHerwigRndm::GetHerwigRandom()->SetSeed(clock()+gSystem->GetPid());}
159 
160  AliStack *stack = fGen->GetStack();
161  stack->Reset();
162  fGen->Generate();
163  const Int_t nprim = stack->GetNprimary();
164  // reject if partons are missing from stack for some reason
165  if(nprim < 8) return;
166  if(fPythiaInfo) {
167  TParticle *part6 = stack->Particle(6);
168  TParticle *part7 = stack->Particle(7);
169 
170  fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
171  fPythiaInfo->SetParton6(part6->Pt(), part6->Eta(), part6->Phi(), part6->GetMass());
172 
173  fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
174  fPythiaInfo->SetParton7(part7->Pt(), part7->Eta(), part7->Phi(), part7->GetMass());
175  }
176 
177  for (Int_t i=0;i<nprim;++i) {
178  if (!stack->IsPhysicalPrimary(i))
179  continue;
180  TParticle *part = stack->Particle(i);
181  TParticlePDG *pdg = part->GetPDG(1);
182  if (!pdg)
183  continue;
184  Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
185  if (fChargedOnly && c==0) continue;
186  Double_t pt = part->Pt();
187  Double_t eta = part->Eta();
188  Double_t phi = part->Phi();
189  if (eta<fEtaMin)
190  continue;
191  if (eta>fEtaMax)
192  continue;
193  if (phi<fPhiMin)
194  continue;
195  if (phi>fPhiMax)
196  continue;
197  if (pt<fPtMin)
198  continue;
199  if (pt>fPtMax)
200  continue;
201  Double_t mass = part->GetMass();
202  if(fMassless) mass = 0.;
204  fHistPt->Fill(pt);
205  fHistEtaPhi->Fill(eta,phi);
206  AddTrack(pt, eta, phi,0,0,0,0,0,0,c,mass);
207  }
208 
210 }
211 
212 //________________________________________________________________________
214  //Get PYTHIA info: pt-hard, x-section, trials
215 
216  if (!fQAhistos)
217  return;
218 
219  AliRunLoader *rl = AliRunLoader::Instance();
220  AliGenPythiaEventHeader *genPP=0x0;
221  AliGenHerwigEventHeader *genPH=0x0;
222  if(fGenType<3) genPP = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
223  if(fGenType==3) genPH = dynamic_cast<AliGenHerwigEventHeader*>(rl->GetHeader()->GenEventHeader());
224  if(fGenType<3 && genPP) {
225  Float_t xsec = genPP->GetXsection();
226  Int_t trials = genPP->Trials();
227  Float_t pthard = genPP->GetPtHard();
228  Float_t ptWeight=genPP->EventWeight();
229 
231  fHistXsection->Fill(0.5,xsec);
232  fHistTrials->Fill(0.5,trials);
233  fHistPtHard->Fill(pthard);
234  }
235 
236  if(fGenType==3 && genPH) {
237  Float_t xsec = genPH->Weight();
238  Int_t trials = genPH->Trials();
239  Float_t ptWeight=genPH->Weight();
241  fHistXsection->Fill(0.5,xsec);
242  fHistTrials->Fill(0.5,trials);
243 
244  }
245 
246 
247 }
Int_t pdg
Int_t fGenType
generator type. 0=pythia, 1=qpythia,2=pyquen, 3=herwig6.5
double Double_t
Definition: External.C:58
void SetParton7(Float_t pt, Float_t eta, Float_t phi, Float_t mass=0)
Definition: External.C:236
Float_t fPtMin
pt minimum value
Bool_t ExecOnce()
generate a particle with random eta,phi, and correlated pt,mass values
void SetPartonFlag7(Int_t flag7)
Double_t mass
TList * fOutput
! output list for QA histograms
AliRunLoader * rl
TSystem * gSystem
TCanvas * c
Definition: TestFitELoss.C:172
Base class for embedding into an event.
AliStack * stack
TH2F * fHistEtaPhi
pT spectrum of generated particles
TRandom * gRandom
Float_t fEtaMax
eta maximum value
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
void SetParton6(Float_t pt, Float_t eta, Float_t phi, Float_t mass=0)
Class for embedding a generated monte carlo event into a data event.
TClonesArray * fOutTracks
! output track collection
TProfile * fHistXsection
trials from generator
Int_t fNTracks
how many tracks are being processed
Float_t fPhiMin
phi minimum value
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
TString fTracksName
name of track collection
Bool_t fQAhistos
draw QA histograms
Bool_t fCopyArray
whether or not the array will be copied to a new one before modelling
void SetPartonFlag6(Int_t flag6)
TH1F * fHistTrials
eta-phi of generated particles
AliEmcalPythiaInfo * fPythiaInfo
! Info on original partons:PDG,pt, eta, phi and pythia event weight
TH1 * fHistPtHard
x-section from generator
Float_t fEtaMin
eta minimum value
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
TString fPythiaInfoName
name of pythia info
Declaration of class AliEmcalPythiaInfo.
void SetSuffix(const char *s)
bool Bool_t
Definition: External.C:53
void SetPythiaEventWeight(Float_t ptWeight)
AliPicoTrack * AddTrack(Double_t pt=-999, Double_t eta=-999, Double_t phi=-999, Byte_t type=0, Double_t etaemc=0, Double_t phiemc=0, Double_t ptemc=0, Bool_t ise=kFALSE, Int_t label=0, Short_t charge=1, Double_t mass=0.1396)
add a cluster (copy)
Float_t fPhiMax
phi maximum value
Float_t fPtMax
pt maximum value