]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
task to do the isolation and correlation analysis only at generator level with the...
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaGeneratorKine.cxx
CommitLineData
7b2086c3 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
16//_________________________________________________________________________
17// Do photon/pi0 analysis for isolation and correlation
18// at the generator level. Only for kine stack (ESDs)
19//
20//
21// -- Author: Gustavo Conesa (LPSC-CNRS-Grenoble)
22//////////////////////////////////////////////////////////////////////////////
23
24// --- ROOT system ---
25#include "TH2F.h"
26#include "TParticle.h"
27#include "TDatabasePDG.h"
28
29//---- ANALYSIS system ----
30#include "AliAnaGeneratorKine.h"
31#include "AliStack.h"
32#include "AliGenPythiaEventHeader.h"
33
34ClassImp(AliAnaGeneratorKine)
35
36
37//__________________________________________
38AliAnaGeneratorKine::AliAnaGeneratorKine() :
39AliAnaCaloTrackCorrBaseClass(),
40fStack(0),
41fParton2(0), fParton3(0),
42fParton6(0), fParton7(0),
43fJet6(), fJet7(),
44fPtHard(0),
45fhPtHard(0), fhPtParton(0), fhPtJet(0),
46fhPtPartonPtHard(0), fhPtJetPtHard(0), fhPtJetPtParton(0),
47fhPtPhoton(0), fhPtPi0(0)
48{
49 //Default Ctor
50
51 //Initialize parameters
52 InitParameters();
53
54 for(Int_t i = 0; i < 4; i++)
55 {
56 fhPtPhotonLeading[i] = fhPtPi0Leading[i] = 0;
57 fhPtPhotonLeadingIsolated[i] = fhPtPi0LeadingIsolated[i] = 0;
58 fhZHardPhotonLeading[i] = fhZHardPi0Leading[i] = 0;
59 fhZHardPhotonLeadingIsolated[i] = fhZHardPi0LeadingIsolated[i] = 0;
60 fhZPartonPhotonLeading[i] = fhZPartonPi0Leading[i] = 0;
61 fhZPartonPhotonLeadingIsolated[i] = fhZPartonPi0LeadingIsolated[i] = 0;
62 fhZJetPhotonLeading[i] = fhZJetPi0Leading[i] = 0;
63 fhZJetPhotonLeadingIsolated[i] = fhZJetPi0LeadingIsolated[i] = 0;
64 fhXEPhotonLeading[i] = fhXEPi0Leading[i] = 0;
65 fhXEPhotonLeadingIsolated[i] = fhXEPi0LeadingIsolated[i] = 0;
66 fhXEUEPhotonLeading[i] = fhXEUEPi0Leading[i] = 0;
67 fhXEUEPhotonLeadingIsolated[i] = fhXEUEPi0LeadingIsolated[i] = 0;
68 }
69
70}
71
72//_______________________________________________________________________________
73void AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger,
74 const Int_t indexTrig,
75 const Int_t pdgTrig,
76 const Bool_t leading[4],
77 const Bool_t isolated[4],
78 Int_t & iparton )
79{
80 //Correlate trigger with partons or jets, get z
81
82 //Get the index of the mother
83 iparton = (fStack->Particle(indexTrig))->GetFirstMother();
84 TParticle * mother = fStack->Particle(iparton);
85 while (iparton > 7)
86 {
87 iparton = mother->GetFirstMother();
88 mother = fStack->Particle(iparton);
89 }
90
91 //printf("Mother is parton %d with pdg %d\n",imom,mother->GetPdgCode());
92
93 if(iparton < 6)
94 {
95 //printf("This particle is not from hard process - pdg %d, parton index %d\n",pdgTrig, iparton);
96 return;
97 }
98
99 Float_t ptTrig = trigger.Pt();
100 Float_t partonPt = fParton6->Pt();
101 Float_t jetPt = fJet6.Pt();
102 if(iparton==7)
103 {
104 partonPt = fParton6->Pt();
105 jetPt = fJet6.Pt();
106 }
107
108
109 fhPtPartonPtHard->Fill(fPtHard, partonPt/fPtHard);
110 fhPtJetPtHard ->Fill(fPtHard, jetPt/fPtHard);
111 fhPtJetPtParton ->Fill(fPtHard, jetPt/partonPt);
112
113 Float_t zHard = ptTrig / fPtHard;
114 Float_t zPart = ptTrig / partonPt;
115 Float_t zJet = ptTrig / jetPt;
116
117 //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard);
118
119 //printf("Z : hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
120
121 for( Int_t i = 0; i < 4; i++ )
122 {
123 if(pdgTrig==111)
124 {
125 if(leading[i])
126 {
127 fhZHardPi0Leading[i] ->Fill(ptTrig,zHard);
128 fhZPartonPi0Leading[i]->Fill(ptTrig,zPart);
129 fhZJetPi0Leading[i] ->Fill(ptTrig,zJet );
130
131 if(isolated[i])
132 {
133 fhZHardPi0LeadingIsolated[i] ->Fill(ptTrig,zHard);
134 fhZPartonPi0LeadingIsolated[i]->Fill(ptTrig,zPart);
135 fhZJetPi0LeadingIsolated[i] ->Fill(ptTrig,zJet);
136 }
137 }
138 }// pi0
139 else if(pdgTrig==22)
140 {
141 if(leading[i])
142 {
143 fhZHardPhotonLeading[i] ->Fill(ptTrig,zHard);
144 fhZPartonPhotonLeading[i]->Fill(ptTrig,zPart);
145 fhZJetPhotonLeading[i] ->Fill(ptTrig,zJet );
146
147 if(isolated[i])
148 {
149 fhZHardPhotonLeadingIsolated[i] ->Fill(ptTrig,zHard);
150 fhZPartonPhotonLeadingIsolated[i]->Fill(ptTrig,zPart);
151 fhZJetPhotonLeadingIsolated[i] ->Fill(ptTrig,zJet);
152 }
153 }
154 } // photon
155 } // conditions loop
156}
157
158
159//____________________________________________________
160TList * AliAnaGeneratorKine::GetCreateOutputObjects()
161{
162 // Create histograms to be saved in output file
163
164 TList * outputContainer = new TList() ;
165 outputContainer->SetName("GenKineHistos") ;
166
167 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
168 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
169 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
170
171 fhPtHard = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax);
172 fhPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
173 outputContainer->Add(fhPtHard);
174
175 fhPtParton = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax);
176 fhPtParton->SetXTitle("p_{T}^{parton} (GeV/c)");
177 outputContainer->Add(fhPtParton);
178
179 fhPtJet = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax);
180 fhPtJet->SetXTitle("p_{T}^{jet} (GeV/c)");
181 outputContainer->Add(fhPtJet);
182
183 fhPtPartonPtHard = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
184 fhPtPartonPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
185 fhPtPartonPtHard->SetYTitle("p_{T}^{parton}/p_{T}^{hard}");
186 outputContainer->Add(fhPtPartonPtHard);
187
188 fhPtJetPtHard = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
189 fhPtJetPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
190 fhPtJetPtHard->SetYTitle("p_{T}^{jet}/p_{T}^{hard}");
191 outputContainer->Add(fhPtJetPtHard);
192
193 fhPtJetPtParton = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
194 fhPtJetPtParton->SetXTitle("p_{T}^{hard} (GeV/c)");
195 fhPtJetPtParton->SetYTitle("p_{T}^{jet}/p_{T}^{parton}");
196 outputContainer->Add(fhPtJetPtParton);
197
198
199 fhPtPhoton = new TH1F("hPtPhoton","Input Photon",nptbins,ptmin,ptmax);
200 fhPtPhoton->SetXTitle("p_{T} (GeV/c)");
201 outputContainer->Add(fhPtPhoton);
202
203 fhPtPi0 = new TH1F("hPtPi0","Input Pi0",nptbins,ptmin,ptmax);
204 fhPtPi0->SetXTitle("p_{T} (GeV/c)");
205 outputContainer->Add(fhPtPi0);
206
207 TString name[] = {"","_EMC","_Photon","_EMC_Photon"};
208 TString title[] = {"",", neutral in EMCal",", neutral only photon like",", neutral in EMCal and only photon like"};
209
210 for(Int_t i = 0; i < 4; i++)
211 {
212
213 // Pt
214
215 fhPtPhotonLeading[i] = new TH1F(Form("hPtPhotonLeading%s",name[i].Data()),
216 Form("Photon : Leading of all particles%s",title[i].Data()),
217 nptbins,ptmin,ptmax);
218 fhPtPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
219 outputContainer->Add(fhPtPhotonLeading[i]);
220
221 fhPtPi0Leading[i] = new TH1F(Form("hPtPi0Leading%s",name[i].Data()),
222 Form("Pi0 : Leading of all particles%s",title[i].Data()),
223 nptbins,ptmin,ptmax);
224 fhPtPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
225 outputContainer->Add(fhPtPi0Leading[i]);
226
227 fhPtPhotonLeadingIsolated[i] = new TH1F(Form("hPtPhotonLeadingIsolated%s",name[i].Data()),
228 Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
229 nptbins,ptmin,ptmax);
230 fhPtPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
231 outputContainer->Add(fhPtPhotonLeadingIsolated[i]);
232
233 fhPtPi0LeadingIsolated[i] = new TH1F(Form("hPtPi0LeadingIsolated%s",name[i].Data()),
234 Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
235 nptbins,ptmin,ptmax);
236 fhPtPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
237 outputContainer->Add(fhPtPi0LeadingIsolated[i]);
238
239
240
241 // zHard
242
243 fhZHardPhotonLeading[i] = new TH2F(Form("hZHardPhotonLeading%s",name[i].Data()),
244 Form("Z-Hard of Photon : Leading of all particles%s",title[i].Data()),
245 nptbins,ptmin,ptmax,200,0,2);
246 fhZHardPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
247 fhZHardPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
248 outputContainer->Add(fhZHardPhotonLeading[i]);
249
250 fhZHardPi0Leading[i] = new TH2F(Form("hZHardPi0Leading%s",name[i].Data()),
251 Form("Z-Hard of Pi0 : Leading of all particles%s",title[i].Data()),
252 nptbins,ptmin,ptmax,200,0,2);
253 fhZHardPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
254 fhZHardPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
255 outputContainer->Add(fhZHardPi0Leading[i]);
256
257 fhZHardPhotonLeadingIsolated[i] = new TH2F(Form("hZHardPhotonLeadingIsolated%s",name[i].Data()),
258 Form("Z-Hard of Photon : Leading of all particles%s, isolated",title[i].Data()),
259 nptbins,ptmin,ptmax,200,0,2);
260 fhZHardPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
261 fhZHardPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
262 outputContainer->Add(fhZHardPhotonLeadingIsolated[i]);
263
264 fhZHardPi0LeadingIsolated[i] = new TH2F(Form("hZHardPi0LeadingIsolated%s",name[i].Data()),
265 Form("Z-Hard of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
266 nptbins,ptmin,ptmax,200,0,2);
267 fhZHardPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
268 fhZHardPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
269 outputContainer->Add(fhZHardPi0LeadingIsolated[i]);
270
271 // zHard
272
273 fhZPartonPhotonLeading[i] = new TH2F(Form("hZPartonPhotonLeading%s",name[i].Data()),
274 Form("Z-Parton of Photon : Leading of all particles%s",title[i].Data()),
275 nptbins,ptmin,ptmax,200,0,2);
276 fhZPartonPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
277 fhZPartonPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
278 outputContainer->Add(fhZPartonPhotonLeading[i]);
279
280 fhZPartonPi0Leading[i] = new TH2F(Form("hZPartonPi0Leading%s",name[i].Data()),
281 Form("Z-Parton of Pi0 : Leading of all particles%s",title[i].Data()),
282 nptbins,ptmin,ptmax,200,0,2);
283 fhZPartonPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
284 fhZPartonPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
285 outputContainer->Add(fhZPartonPi0Leading[i]);
286
287 fhZPartonPhotonLeadingIsolated[i] = new TH2F(Form("hZPartonPhotonLeadingIsolated%s",name[i].Data()),
288 Form("Z-Parton of Photon : Leading of all particles%s, isolated",title[i].Data()),
289 nptbins,ptmin,ptmax,200,0,2);
290 fhZPartonPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
291 fhZPartonPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
292 outputContainer->Add(fhZPartonPhotonLeadingIsolated[i]);
293
294 fhZPartonPi0LeadingIsolated[i] = new TH2F(Form("hZPartonPi0LeadingIsolated%s",name[i].Data()),
295 Form("Z-Parton of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
296 nptbins,ptmin,ptmax,200,0,2);
297 fhZPartonPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
298 fhZPartonPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
299 outputContainer->Add(fhZPartonPi0LeadingIsolated[i]);
300
301
302 // zJet
303
304 fhZJetPhotonLeading[i] = new TH2F(Form("hZJetPhotonLeading%s",name[i].Data()),
305 Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
306 nptbins,ptmin,ptmax,200,0,2);
307 fhZJetPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
308 fhZJetPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
309 outputContainer->Add(fhZJetPhotonLeading[i]);
310
311 fhZJetPi0Leading[i] = new TH2F(Form("hZJetPi0Leading%s",name[i].Data()),
312 Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
313 nptbins,ptmin,ptmax,200,0,2);
314 fhZJetPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
315 fhZJetPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
316 outputContainer->Add(fhZJetPi0Leading[i]);
317
318 fhZJetPhotonLeadingIsolated[i] = new TH2F(Form("hZJetPhotonLeadingIsolated%s",name[i].Data()),
319 Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
320 nptbins,ptmin,ptmax,200,0,2);
321 fhZJetPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
322 fhZJetPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
323 outputContainer->Add(fhZJetPhotonLeadingIsolated[i]);
324
325 fhZJetPi0LeadingIsolated[i] = new TH2F(Form("hZJetPi0LeadingIsolated%s",name[i].Data()),
326 Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
327 nptbins,ptmin,ptmax,200,0,2);
328 fhZJetPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
329 fhZJetPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
330 outputContainer->Add(fhZJetPi0LeadingIsolated[i]);
331
332
333 // XE
334
335 fhXEPhotonLeading[i] = new TH2F(Form("hXEPhotonLeading%s",name[i].Data()),
336 Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
337 nptbins,ptmin,ptmax,200,0,2);
338 fhXEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
339 fhXEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
340 outputContainer->Add(fhXEPhotonLeading[i]);
341
342 fhXEPi0Leading[i] = new TH2F(Form("hXEPi0Leading%s",name[i].Data()),
343 Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
344 nptbins,ptmin,ptmax,200,0,2);
345 fhXEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
346 fhXEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
347 outputContainer->Add(fhXEPi0Leading[i]);
348
349 fhXEPhotonLeadingIsolated[i] = new TH2F(Form("hXEPhotonLeadingIsolated%s",name[i].Data()),
350 Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
351 nptbins,ptmin,ptmax,200,0,2);
352 fhXEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
353 fhXEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
354 outputContainer->Add(fhXEPhotonLeadingIsolated[i]);
355
356 fhXEPi0LeadingIsolated[i] = new TH2F(Form("hXEPi0LeadingIsolated%s",name[i].Data()),
357 Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
358 nptbins,ptmin,ptmax,200,0,2);
359 fhXEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
360 fhXEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
361 outputContainer->Add(fhXEPi0LeadingIsolated[i]);
362
363
364 // XE from UE
365
366 fhXEUEPhotonLeading[i] = new TH2F(Form("hXEUEPhotonLeading%s",name[i].Data()),
367 Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
368 nptbins,ptmin,ptmax,200,0,2);
369 fhXEUEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
370 fhXEUEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
371 outputContainer->Add(fhXEUEPhotonLeading[i]);
372
373 fhXEUEPi0Leading[i] = new TH2F(Form("hXEUEPi0Leading%s",name[i].Data()),
374 Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
375 nptbins,ptmin,ptmax,200,0,2);
376 fhXEUEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
377 fhXEUEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
378 outputContainer->Add(fhXEUEPi0Leading[i]);
379
380 fhXEUEPhotonLeadingIsolated[i] = new TH2F(Form("hXEUEPhotonLeadingIsolated%s",name[i].Data()),
381 Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
382 nptbins,ptmin,ptmax,200,0,2);
383 fhXEUEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
384 fhXEUEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
385 outputContainer->Add(fhXEUEPhotonLeadingIsolated[i]);
386
387 fhXEUEPi0LeadingIsolated[i] = new TH2F(Form("hXEUEPi0LeadingIsolated%s",name[i].Data()),
388 Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
389 nptbins,ptmin,ptmax,200,0,2);
390 fhXEUEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
391 fhXEUEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
392 outputContainer->Add(fhXEUEPi0LeadingIsolated[i]);
393
394 }
395
396 return outputContainer;
397
398}
399
400//____________________________________________
401void AliAnaGeneratorKine::GetPartonsAndJets()
402{
403 // Fill data members with partons,jets and generated pt hard
404
405 fStack = GetMCStack() ;
406
407 if(!fStack)
408 {
409 printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - No Stack available, STOP\n");
410 abort();
411 }
412
413 fParton2 = fStack->Particle(2) ;
414 fParton3 = fStack->Particle(3) ;
415 fParton6 = fStack->Particle(6) ;
416 fParton7 = fStack->Particle(7) ;
417
418 Float_t p6phi = fParton6->Phi();
419 if(p6phi < 0) p6phi +=TMath::TwoPi();
420 Float_t p7phi = fParton7->Phi();
421 if(p7phi < 0) p7phi +=TMath::TwoPi();
422
423 //printf("parton6: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton6->Pt(),fParton6->Eta(),p6phi, fParton6->GetPdgCode());
424 //printf("parton7: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton7->Pt(),fParton7->Eta(),p7phi, fParton7->GetPdgCode());
425
426 // Get the jets, only for pythia
427 if(!strcmp(GetReader()->GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
428 {
429 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetReader()->GetGenEventHeader();
430
431 fPtHard = pygeh->GetPtHard();
432
433 //printf("pt Hard %2.2f\n",fPtHard);
434
435 const Int_t nTriggerJets = pygeh->NTriggerJets();
436
437 Float_t tmpjet[]={0,0,0,0};
438
439 // select the closest jet to parton
440 Float_t jet7R = 100;
441 Float_t jet6R = 100;
442
443 for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
444 {
445 pygeh->TriggerJet(ijet, tmpjet);
446
447 TLorentzVector jet(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
448 Float_t jphi = jet.Phi();
449 if(jphi < 0) jphi +=TMath::TwoPi();
450
451 Double_t radius6 = GetIsolationCut()->Radius(fParton6->Eta(), p6phi, jet.Eta() , jphi) ;
452 Double_t radius7 = GetIsolationCut()->Radius(fParton7->Eta(), p7phi, jet.Eta() , jphi) ;
453
454 //printf("jet %d: pt %2.2f, eta %2.2f, phi %2.2f, r6 %2.2f, r7 %2.2f\n",ijet,jet.Pt(),jet.Eta(),jphi,radius6, radius7);
455
456 if (radius6 < jet6R)
457 {
458 jet6R = radius6;
459 fJet6 = jet;
460
461 }
462 if (radius7 < jet7R)
463 {
464 jet7R = radius7;
465 fJet7 = jet;
466 }
467
468 } // jet loop
469
470 //printf("jet6: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet6.Pt(),fJet6.Eta(),fJet6.Phi());
471 //printf("jet7: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet7.Pt(),fJet7.Eta(),fJet6.Phi());
472
473 } // pythia header
474
475 fhPtHard ->Fill(fPtHard);
476 fhPtJet ->Fill(fJet6.Pt());
477 fhPtJet ->Fill(fJet7.Pt());
478 fhPtParton ->Fill(fParton6->Pt());
479 fhPtParton ->Fill(fParton7->Pt());
480
481}
482
483//___________________________________________________________
484void AliAnaGeneratorKine::GetXE(const TLorentzVector trigger,
485 const Int_t indexTrig,
486 const Int_t pdgTrig,
487 const Bool_t leading[4],
488 const Bool_t isolated[4],
489 const Int_t iparton)
490{
491
492 // Calculate the real XE and the UE XE
493
494 Float_t ptThresTrack = 0.2;
495
496 Float_t ptTrig = trigger.Pt();
497 Float_t etaTrig = trigger.Eta();
498 Float_t phiTrig = trigger.Phi();
499 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
500
501 //Loop on primaries, start from position 8, no partons
502 for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
503 {
504
505 TParticle * particle = fStack->Particle(ipr) ;
506
507 if(ipr==indexTrig) continue;
508
509 Int_t pdg = particle->GetPdgCode();
510 Int_t status = particle->GetStatusCode();
511
512 // Compare trigger with final state particles
513 if( status != 1) continue ;
514
515 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
516
517 if(charge==0) continue; // construct xe only with charged
518
519 Float_t pt = particle->Pt();
520 Float_t eta = particle->Eta();
521 Float_t phi = particle->Phi();
522 if(phi < 0 ) phi += TMath::TwoPi();
523
524 if( pt < ptThresTrack) continue ;
525
526 if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
527
528 //Isolation
529 Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
530
531 Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
532
533 //Get the index of the mother
534 Int_t ipartonAway = particle->GetFirstMother();
535 TParticle * mother = fStack->Particle(ipartonAway);
536 while (ipartonAway > 7)
537 {
538 ipartonAway = mother->GetFirstMother();
539 mother = fStack->Particle(ipartonAway);
540 }
541
542 if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway)
543 {
544 //printf("xE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
545 if(radius > 1 ) continue; // avoid particles too far from trigger
546
547 for( Int_t i = 0; i < 4; i++ )
548 {
549 if(pdgTrig==111)
550 {
551 if(leading[i])
552 {
553 fhXEPi0Leading[i] ->Fill(ptTrig,xe);
554
555 if(isolated[i])
556 {
557 fhXEPi0LeadingIsolated[i] ->Fill(ptTrig,xe);
558 }
559 }
560 }// pi0
561 else if(pdgTrig==22)
562 {
563 if(leading[i])
564 {
565 fhXEPhotonLeading[i] ->Fill(ptTrig,xe);
566
567 if(isolated[i])
568 {
569 fhXEPhotonLeadingIsolated[i] ->Fill(ptTrig,xe);
570 }
571 }
572 } // photon
573 } // conditions loop
574 } // Away side
575
576 if(ipartonAway!=6 && ipartonAway!=7)
577 {
578
579 //printf("xE UE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
580
581 for( Int_t i = 0; i < 4; i++ )
582 {
583 if(pdgTrig==111)
584 {
585 if(leading[i])
586 {
587 fhXEUEPi0Leading[i] ->Fill(ptTrig,xe);
588
589 if(isolated[i])
590 {
591 fhXEUEPi0LeadingIsolated[i] ->Fill(ptTrig,xe);
592 }
593 }
594 }// pi0
595 else if(pdgTrig==22)
596 {
597 if(leading[i])
598 {
599 fhXEUEPhotonLeading[i] ->Fill(ptTrig,xe);
600
601 if(isolated[i])
602 {
603 fhXEUEPhotonLeadingIsolated[i] ->Fill(ptTrig,xe);
604 }
605 }
606 } // photon
607 } // conditions loop
608 } // Away side
609
610 } // primary loop
611
612
613}
614
615//________________________________________
616void AliAnaGeneratorKine::InitParameters()
617{
618
619 //Initialize the parameters of the analysis.
620 AddToHistogramsName("AnaGenKine_");
621
622}
623
624//___________________________________________________________________________
625void AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
626 const Int_t indexTrig,
627 const Int_t pdgTrig,
628 Bool_t leading[4],
629 Bool_t isolated[4])
630{
631 // Check if the trigger is the leading particle
632 // In case of neutral particles check all neutral or neutral in EMCAL acceptance
633
634 Float_t ptMaxCharged = 0; // all charged
635 Float_t ptMaxNeutral = 0; // all neutral
636 Float_t ptMaxNeutEMCAL = 0; // for neutral, select them in EMCAL acceptance
637 Float_t ptMaxNeutPhot = 0; // for neutral, take only photons
638 Float_t ptMaxNeutEMCALPhot = 0; // for neutral, take only photons in EMCAL acceptance
639
640 leading[0] = 0;
641 leading[1] = 0;
642 leading[2] = 0;
643 leading[3] = 0;
644
645 isolated[0] = 0;
646 isolated[1] = 0;
647 isolated[2] = 0;
648 isolated[3] = 0;
649
650 Float_t ptTrig = trigger.Pt();
651 Float_t etaTrig = trigger.Eta();
652 Float_t phiTrig = trigger.Phi();
653 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
654
655 // Minimum track or cluster energy
656 Float_t ptThresTrack = 0.2;
657 Float_t ptThresCalo = 0.3;
658
659 //Isolation cuts
660 Float_t ptThresIC = 0.5;
661 Float_t rThresIC = 0.4;
662 Int_t nICTrack = 0;
663 Int_t nICNeutral = 0;
664 Int_t nICNeutEMCAL = 0;
665 Int_t nICNeutPhot = 0;
666 Int_t nICNeutEMCALPhot = 0;
667
668 //Loop on primaries, start from position 8, no partons
669 for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
670 {
671 if(ipr == indexTrig) continue;
672 TParticle * particle = fStack->Particle(ipr) ;
673
674 Int_t imother= particle->GetFirstMother();
675 //printf("Leading ipr %d - mother %d\n",ipr, imother);
676
677 if(imother==indexTrig) continue ;
678
679 Int_t pdg = particle->GetPdgCode();
680 Int_t status = particle->GetStatusCode();
681
682 // Compare trigger with final state particles
683 if( status != 1) continue ;
684
685 Float_t pt = particle->Pt();
686 Float_t eta = particle->Eta();
687 Float_t phi = particle->Phi();
688 if(phi < 0 ) phi += TMath::TwoPi();
689
690 if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
691
692 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
693
694 //Isolation
695 Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
696
697 if(charge==0)
698 {
699 if(pt < ptThresCalo) continue ;
700
701 if( ptMaxNeutral < pt ) ptMaxNeutral = pt;
702 if( pt > ptThresIC && radius < rThresIC ) nICNeutral++ ;
703
704 Bool_t phPDG = kFALSE;
705 if(pdg==22 || pdg==111) phPDG = kTRUE;
706
707 //if(pt > ptTrig) printf(" --- pdg %d, phPDG %d pT %2.2f, pTtrig %2.2f, eta %2.2f, phi %2.2f\n",pdg,phPDG,pt,ptTrig,particle->Eta(), particle->Phi()*TMath::RadToDeg());
708 if(phPDG)
709 {
710 if( ptMaxNeutPhot < pt) ptMaxNeutPhot = pt;
711 if( pt > ptThresIC && radius < rThresIC ) nICNeutPhot++ ;
712 }
713
714 //EMCAL acceptance
715 Bool_t inEMCAL = kTRUE;
716 if(TMath::Abs(particle->Eta()) > 0.7
717 || particle->Phi() > TMath::DegToRad()*180
718 || particle->Phi() < TMath::DegToRad()*80 ) inEMCAL = kFALSE ;
719
720 if(inEMCAL)
721 {
722 if( ptMaxNeutEMCAL < pt ) ptMaxNeutEMCAL = pt;
723 if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCAL++ ;
724
725 if(phPDG)
726 {
727 if( ptMaxNeutEMCALPhot < pt ) ptMaxNeutEMCALPhot = pt;
728 if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCALPhot++ ;
729 }
730 }
731 }
732 else
733 {
734 if( pt < ptThresTrack) continue ;
735
736 if( ptMaxCharged < pt ) ptMaxCharged = pt;
737
738 if( pt > ptThresIC && radius < rThresIC )
739 {
740 //printf("UE track? pTtrig %2.2f, pt %2.2f, etaTrig %2.2f, eta %2.2f, phiTrig %2.2f, phi %2.2f, radius %2.2f\n",
741 // ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
742 nICTrack++ ;
743 }
744 }
745
746 } // particle loop
747
748 // Leding decision
749 if(ptTrig > ptMaxCharged)
750 {
751 //printf("pt charged %2.2f, pt neutral %2.2f, pt neutral emcal %2.2f, pt photon %2.2f, pt photon emcal %2.2f\n",
752 // ptMaxCharged, ptMaxNeutral, ptMaxNeutEMCAL,ptMaxNeutPhot, ptMaxNeutEMCALPhot);
753 if(ptTrig > ptMaxNeutral ) leading[0] = kTRUE ;
754 if(ptTrig > ptMaxNeutEMCAL ) leading[1] = kTRUE ;
755 if(ptTrig > ptMaxNeutPhot ) leading[2] = kTRUE ;
756 if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
757 }
758
759 //printf("N in cone over threshold : tracks %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n",
760 // nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
761
762 // Isolation decision
763 if(nICTrack == 0)
764 {
765 if(nICNeutral == 0 ) isolated[0] = kTRUE ;
766 if(nICNeutEMCAL == 0 ) isolated[1] = kTRUE ;
767 if(nICNeutPhot == 0 ) isolated[2] = kTRUE ;
768 if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
769 }
770
771 // Fill histograms if conditions apply for all 4 cases
772 for( Int_t i = 0; i < 4; i++ )
773 {
774 if(pdgTrig==111)
775 {
776 if(leading[i])
777 {
778 fhPtPi0Leading[i]->Fill(ptTrig);
779 if(isolated[i]) fhPtPi0LeadingIsolated[i]->Fill(ptTrig);
780 }
781 }// pi0
782 else if(pdgTrig==22)
783 {
784 if(leading[i])
785 {
786 fhPtPhotonLeading[i]->Fill(ptTrig);
787 if(isolated[i]) fhPtPhotonLeadingIsolated[i]->Fill(ptTrig);
788 }
789 } // photon
790 } // conditions loop
791
792}
793
794//_____________________________________________________
795void AliAnaGeneratorKine::MakeAnalysisFillHistograms()
796{
797 //Particle-Parton Correlation Analysis, fill histograms
798
799 TLorentzVector trigger;
800
801 GetPartonsAndJets();
802
803 for(Int_t ipr = 0; ipr < fStack->GetNprimary(); ipr ++ )
804 {
805 TParticle * particle = fStack->Particle(ipr) ;
806
807 Int_t pdgTrig = particle->GetPdgCode();
808 Int_t statusTrig = particle->GetStatusCode();
809 Int_t imother = particle->GetFirstMother();
810 Float_t ptTrig = particle->Pt();
811
812 // Select final state photons (prompt, fragmentation) or pi0s
813
814 //Check the origin of the photon, accept if prompt or initial/final state radiation
815 if(pdgTrig == 22 && statusTrig == 1)
816 {
817 if(imother > 8) continue;
818 }
819 // If not photon, trigger on pi0
820 else if(pdgTrig != 111) continue;
821
822 // Acceptance and kinematical cuts
823 if( ptTrig < 8 ) continue ;
824
825 //EMCAL acceptance, a bit less
826 if(TMath::Abs(particle->Eta()) > 0.6) continue ;
827 if(particle->Phi() > TMath::DegToRad()*176) continue ;
828 if(particle->Phi() < TMath::DegToRad()*74 ) continue ;
829
830// printf("Particle %d : pdg %d status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f \n",
831// ipr, pdgTrig, statusTrig, imother, ptTrig, particle->Eta(), particle->Phi()*TMath::RadToDeg());
832
833// if(pdgTrig==111)
834// {
835// printf("\t pi0 daughters %d, %d\n", particle->GetDaughter(0), particle->GetDaughter(1));
836// }
837
838
839 if (pdgTrig==22 ) fhPtPhoton->Fill(ptTrig);
840 else if(pdgTrig==111) fhPtPi0 ->Fill(ptTrig);
841
842 // Check if it is leading
843 Bool_t leading[4] ;
844 Bool_t isolated[4] ;
845
846 particle->Momentum(trigger);
847
848 IsLeadingAndIsolated(trigger, ipr, pdgTrig, leading, isolated);
849
850 Int_t iparton = -1;
851 CorrelateWithPartonOrJet(trigger, ipr, pdgTrig, leading, isolated, iparton);
852
853 GetXE(trigger,ipr,pdgTrig,leading,isolated,iparton) ;
854
855 }
856
857 if(GetDebug() > 1) printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - End fill histograms \n");
858
859}