]> 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
create fill histograms with DCA cut when not applied before in the reader
[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 {
dbb79a05 56 fhPtPhotonLeading[i] = fhPtPi0Leading[i] = 0;
57 fhPtPhotonLeadingIsolated[i] = fhPtPi0LeadingIsolated[i] = 0;
58 for(Int_t j = 0; j < 2; j++)
59 {
60 fhZHardPhoton[j][i] = fhZHardPi0[j][i] = 0;
61 fhZHardPhotonIsolated[j][i] = fhZHardPi0Isolated[j][i] = 0;
62 fhZPartonPhoton[j][i] = fhZPartonPi0[j][i] = 0;
63 fhZPartonPhotonIsolated[j][i] = fhZPartonPi0Isolated[j][i] = 0;
64 fhZJetPhoton[j][i] = fhZJetPi0[j][i] = 0;
65 fhZJetPhotonIsolated[j][i] = fhZJetPi0Isolated[j][i] = 0;
66 fhXEPhoton[j][i] = fhXEPi0[j][i] = 0;
67 fhXEPhotonIsolated[j][i] = fhXEPi0Isolated[j][i] = 0;
68 fhXEUEPhoton[j][i] = fhXEUEPi0[j][i] = 0;
69 fhXEUEPhotonIsolated[j][i] = fhXEUEPi0Isolated[j][i] = 0;
70
71 fhPtPartonTypeNearPhoton[j][i] = fhPtPartonTypeNearPi0[j][i] = 0;
72 fhPtPartonTypeNearPhotonIsolated[j][i] = fhPtPartonTypeNearPi0Isolated[j][i] = 0;
73
74 fhPtPartonTypeAwayPhoton[j][i] = fhPtPartonTypeAwayPi0[j][i] = 0;
75 fhPtPartonTypeAwayPhotonIsolated[j][i] = fhPtPartonTypeAwayPi0Isolated[j][i] = 0;
76 }
7b2086c3 77 }
78
79}
80
b94e038e 81//___________________________________________________________________________
82Bool_t AliAnaGeneratorKine::CorrelateWithPartonOrJet(TLorentzVector trigger,
83 Int_t indexTrig,
84 Int_t pdgTrig,
85 Bool_t leading[4],
86 Bool_t isolated[4],
2292cf03 87 Int_t & iparton )
7b2086c3 88{
89 //Correlate trigger with partons or jets, get z
90
91 //Get the index of the mother
92 iparton = (fStack->Particle(indexTrig))->GetFirstMother();
93 TParticle * mother = fStack->Particle(iparton);
94 while (iparton > 7)
95 {
96 iparton = mother->GetFirstMother();
2292cf03 97 if(iparton < 0) { printf("AliAnaGeneratorKine::CorrelateWithPartonOrJet() - Negative index, skip event\n"); return kFALSE; }
7b2086c3 98 mother = fStack->Particle(iparton);
99 }
100
101 //printf("Mother is parton %d with pdg %d\n",imom,mother->GetPdgCode());
102
103 if(iparton < 6)
104 {
105 //printf("This particle is not from hard process - pdg %d, parton index %d\n",pdgTrig, iparton);
2292cf03 106 return kFALSE;
7b2086c3 107 }
108
109 Float_t ptTrig = trigger.Pt();
110 Float_t partonPt = fParton6->Pt();
111 Float_t jetPt = fJet6.Pt();
112 if(iparton==7)
113 {
114 partonPt = fParton6->Pt();
115 jetPt = fJet6.Pt();
116 }
117
23fa04c5 118 //Get id of parton in near and away side
119
120 Int_t away = -1;
121 Int_t near = -1;
122 Int_t nearPDG = -1;
123 Int_t awayPDG = -1;
124
125 //printf("parton 6 pdg = %d, parton 7 pdg = %d\n",fParton6->GetPdgCode(),fParton7->GetPdgCode());
126
127 if(iparton==6)
128 {
129 nearPDG = fParton6->GetPdgCode();
130 awayPDG = fParton7->GetPdgCode();
131 }
132 else
133 {
134 nearPDG = fParton7->GetPdgCode();
135 awayPDG = fParton6->GetPdgCode();
136 }
137
138 if (nearPDG == 22) near = 0;
139 else if(nearPDG == 21) near = 1;
140 else near = 2;
141
142 if (awayPDG == 22) away = 0;
143 else if(awayPDG == 21) away = 1;
144 else away = 2;
145
146 for( Int_t i = 0; i < 4; i++ )
147 {
148 if(pdgTrig==111)
149 {
dbb79a05 150
151 fhPtPartonTypeNearPi0[leading[i]][i]->Fill(ptTrig,near);
152 fhPtPartonTypeAwayPi0[leading[i]][i]->Fill(ptTrig,away);
153 if(isolated[i])
154 {
155 fhPtPartonTypeNearPi0Isolated[leading[i]][i]->Fill(ptTrig,near);
156 fhPtPartonTypeAwayPi0Isolated[leading[i]][i]->Fill(ptTrig,away);
23fa04c5 157 }
dbb79a05 158
23fa04c5 159 }// pi0
160 else if(pdgTrig==22)
161 {
dbb79a05 162 fhPtPartonTypeNearPhoton[leading[i]][i]->Fill(ptTrig,near);
163 fhPtPartonTypeAwayPhoton[leading[i]][i]->Fill(ptTrig,away);
164 if(isolated[i])
165 {
166 fhPtPartonTypeNearPhotonIsolated[leading[i]][i]->Fill(ptTrig,near);
167 fhPtPartonTypeAwayPhotonIsolated[leading[i]][i]->Fill(ptTrig,away);
23fa04c5 168 }
dbb79a05 169
23fa04c5 170 } // photon
dbb79a05 171 } // conditions loop
23fa04c5 172
173
174 // RATIOS
7b2086c3 175
176 fhPtPartonPtHard->Fill(fPtHard, partonPt/fPtHard);
177 fhPtJetPtHard ->Fill(fPtHard, jetPt/fPtHard);
178 fhPtJetPtParton ->Fill(fPtHard, jetPt/partonPt);
179
180 Float_t zHard = ptTrig / fPtHard;
181 Float_t zPart = ptTrig / partonPt;
182 Float_t zJet = ptTrig / jetPt;
183
184 //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard);
185
186 //printf("Z : hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
187
188 for( Int_t i = 0; i < 4; i++ )
189 {
190 if(pdgTrig==111)
191 {
dbb79a05 192 fhZHardPi0[leading[i]][i] ->Fill(ptTrig,zHard);
193 fhZPartonPi0[leading[i]][i]->Fill(ptTrig,zPart);
194 fhZJetPi0[leading[i]][i] ->Fill(ptTrig,zJet );
195
196 if(isolated[i])
197 {
198 fhZHardPi0Isolated[leading[i]][i] ->Fill(ptTrig,zHard);
199 fhZPartonPi0Isolated[leading[i]][i]->Fill(ptTrig,zPart);
200 fhZJetPi0Isolated[leading[i]][i] ->Fill(ptTrig,zJet);
7b2086c3 201 }
dbb79a05 202
7b2086c3 203 }// pi0
204 else if(pdgTrig==22)
205 {
dbb79a05 206
207 fhZHardPhoton[leading[i]][i] ->Fill(ptTrig,zHard);
208 fhZPartonPhoton[leading[i]][i]->Fill(ptTrig,zPart);
209 fhZJetPhoton[leading[i]][i] ->Fill(ptTrig,zJet );
210
211 if(isolated[i])
212 {
213 fhZHardPhotonIsolated[leading[i]][i] ->Fill(ptTrig,zHard);
214 fhZPartonPhotonIsolated[leading[i]][i]->Fill(ptTrig,zPart);
215 fhZJetPhotonIsolated[leading[i]][i] ->Fill(ptTrig,zJet);
216 }
217
7b2086c3 218 } // photon
dbb79a05 219 } // conditions loop
2292cf03 220
221 return kTRUE;
dbb79a05 222}
7b2086c3 223
224
225//____________________________________________________
226TList * AliAnaGeneratorKine::GetCreateOutputObjects()
227{
228 // Create histograms to be saved in output file
229
230 TList * outputContainer = new TList() ;
231 outputContainer->SetName("GenKineHistos") ;
232
233 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
234 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
235 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
236
237 fhPtHard = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax);
238 fhPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
239 outputContainer->Add(fhPtHard);
240
241 fhPtParton = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax);
242 fhPtParton->SetXTitle("p_{T}^{parton} (GeV/c)");
243 outputContainer->Add(fhPtParton);
244
245 fhPtJet = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax);
246 fhPtJet->SetXTitle("p_{T}^{jet} (GeV/c)");
247 outputContainer->Add(fhPtJet);
248
249 fhPtPartonPtHard = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
250 fhPtPartonPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
251 fhPtPartonPtHard->SetYTitle("p_{T}^{parton}/p_{T}^{hard}");
252 outputContainer->Add(fhPtPartonPtHard);
253
254 fhPtJetPtHard = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
255 fhPtJetPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
256 fhPtJetPtHard->SetYTitle("p_{T}^{jet}/p_{T}^{hard}");
257 outputContainer->Add(fhPtJetPtHard);
258
259 fhPtJetPtParton = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
260 fhPtJetPtParton->SetXTitle("p_{T}^{hard} (GeV/c)");
261 fhPtJetPtParton->SetYTitle("p_{T}^{jet}/p_{T}^{parton}");
262 outputContainer->Add(fhPtJetPtParton);
263
264
265 fhPtPhoton = new TH1F("hPtPhoton","Input Photon",nptbins,ptmin,ptmax);
266 fhPtPhoton->SetXTitle("p_{T} (GeV/c)");
267 outputContainer->Add(fhPtPhoton);
268
269 fhPtPi0 = new TH1F("hPtPi0","Input Pi0",nptbins,ptmin,ptmax);
270 fhPtPi0->SetXTitle("p_{T} (GeV/c)");
271 outputContainer->Add(fhPtPi0);
272
dbb79a05 273 TString name [] = {"","_EMC","_Photon","_EMC_Photon"};
274 TString title [] = {"",", neutral in EMCal",", neutral only photon like",", neutral in EMCal and only photon like"};
275 TString leading[] = {"NotLeading","Leading"};
276
7b2086c3 277 for(Int_t i = 0; i < 4; i++)
278 {
279
280 // Pt
281
282 fhPtPhotonLeading[i] = new TH1F(Form("hPtPhotonLeading%s",name[i].Data()),
dbb79a05 283 Form("Photon : Leading of all particles%s",title[i].Data()),
284 nptbins,ptmin,ptmax);
7b2086c3 285 fhPtPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
286 outputContainer->Add(fhPtPhotonLeading[i]);
287
288 fhPtPi0Leading[i] = new TH1F(Form("hPtPi0Leading%s",name[i].Data()),
dbb79a05 289 Form("Pi0 : Leading of all particles%s",title[i].Data()),
290 nptbins,ptmin,ptmax);
7b2086c3 291 fhPtPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
dbb79a05 292 outputContainer->Add(fhPtPi0Leading[i]);
7b2086c3 293
294 fhPtPhotonLeadingIsolated[i] = new TH1F(Form("hPtPhotonLeadingIsolated%s",name[i].Data()),
dbb79a05 295 Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
296 nptbins,ptmin,ptmax);
7b2086c3 297 fhPtPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
298 outputContainer->Add(fhPtPhotonLeadingIsolated[i]);
299
300 fhPtPi0LeadingIsolated[i] = new TH1F(Form("hPtPi0LeadingIsolated%s",name[i].Data()),
dbb79a05 301 Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
302 nptbins,ptmin,ptmax);
7b2086c3 303 fhPtPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
dbb79a05 304 outputContainer->Add(fhPtPi0LeadingIsolated[i]);
7b2086c3 305
dbb79a05 306 // Leading or not loop
307 for(Int_t j = 0; j < 2; j++)
308 {
309 // Near side parton
310
311 fhPtPartonTypeNearPhoton[j][i] = new TH2F(Form("hPtPartonTypeNearPhoton%s%s",leading[j].Data(),name[i].Data()),
312 Form("Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
313 nptbins,ptmin,ptmax,3,0,3);
314 fhPtPartonTypeNearPhoton[j][i]->SetXTitle("p_{T} (GeV/c)");
315 fhPtPartonTypeNearPhoton[j][i]->SetYTitle("Parton type");
316 fhPtPartonTypeNearPhoton[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
317 fhPtPartonTypeNearPhoton[j][i]->GetYaxis()->SetBinLabel(2,"g");
318 fhPtPartonTypeNearPhoton[j][i]->GetYaxis()->SetBinLabel(3,"q");
319 outputContainer->Add(fhPtPartonTypeNearPhoton[j][i]);
320
321 fhPtPartonTypeNearPi0[j][i] = new TH2F(Form("hPtPartonTypeNearPi0%s%s",leading[j].Data(),name[i].Data()),
322 Form("Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
323 nptbins,ptmin,ptmax,3,0,3);
324 fhPtPartonTypeNearPi0[j][i]->SetXTitle("p_{T} (GeV/c)");
325 fhPtPartonTypeNearPi0[j][i]->SetYTitle("Parton type");
326 fhPtPartonTypeNearPi0[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
327 fhPtPartonTypeNearPi0[j][i]->GetYaxis()->SetBinLabel(2,"g");
328 fhPtPartonTypeNearPi0[j][i]->GetYaxis()->SetBinLabel(3,"q");
329 outputContainer->Add(fhPtPartonTypeNearPi0[j][i]);
330
331 fhPtPartonTypeNearPhotonIsolated[j][i] = new TH2F(Form("hPtPartonTypeNearPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
332 Form("Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
333 nptbins,ptmin,ptmax,3,0,3);
334 fhPtPartonTypeNearPhotonIsolated[j][i]->SetXTitle("p_{T} (GeV/c)");
335 fhPtPartonTypeNearPhotonIsolated[j][i]->SetYTitle("Parton type");
336 fhPtPartonTypeNearPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
337 fhPtPartonTypeNearPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
338 fhPtPartonTypeNearPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
339 outputContainer->Add(fhPtPartonTypeNearPhotonIsolated[j][i]);
340
341 fhPtPartonTypeNearPi0Isolated[j][i] = new TH2F(Form("hPtPartonTypeNearPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
342 Form("Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
343 nptbins,ptmin,ptmax,3,0,3);
344 fhPtPartonTypeNearPi0Isolated[j][i]->SetXTitle("p_{T} (GeV/c)");
345 fhPtPartonTypeNearPi0Isolated[j][i]->SetYTitle("Parton type");
346 fhPtPartonTypeNearPi0Isolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
347 fhPtPartonTypeNearPi0Isolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
348 fhPtPartonTypeNearPi0Isolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
349 outputContainer->Add(fhPtPartonTypeNearPi0Isolated[j][i]);
350
351
352 // Away side parton
353
354 fhPtPartonTypeAwayPhoton[j][i] = new TH2F(Form("hPtPartonTypeAwayPhoton%s%s",leading[j].Data(),name[i].Data()),
355 Form("Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
356 nptbins,ptmin,ptmax,3,0,3);
357 fhPtPartonTypeAwayPhoton[j][i]->SetXTitle("p_{T} (GeV/c)");
358 fhPtPartonTypeAwayPhoton[j][i]->SetYTitle("Parton type");
359 fhPtPartonTypeAwayPhoton[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
360 fhPtPartonTypeAwayPhoton[j][i]->GetYaxis()->SetBinLabel(2,"g");
361 fhPtPartonTypeAwayPhoton[j][i]->GetYaxis()->SetBinLabel(3,"q");
362 outputContainer->Add(fhPtPartonTypeAwayPhoton[j][i]);
363
364 fhPtPartonTypeAwayPi0[j][i] = new TH2F(Form("hPtPartonTypeAwayPi0%s%s",leading[j].Data(),name[i].Data()),
365 Form("Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
366 nptbins,ptmin,ptmax,3,0,3);
367 fhPtPartonTypeAwayPi0[j][i]->SetXTitle("p_{T} (GeV/c)");
368 fhPtPartonTypeAwayPi0[j][i]->SetYTitle("Parton type");
369 fhPtPartonTypeAwayPi0[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
370 fhPtPartonTypeAwayPi0[j][i]->GetYaxis()->SetBinLabel(2,"g");
371 fhPtPartonTypeAwayPi0[j][i]->GetYaxis()->SetBinLabel(3,"q");
372 outputContainer->Add(fhPtPartonTypeAwayPi0[j][i]);
373
374 fhPtPartonTypeAwayPhotonIsolated[j][i] = new TH2F(Form("hPtPartonTypeAwayPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
375 Form("Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
376 nptbins,ptmin,ptmax,3,0,3);
377 fhPtPartonTypeAwayPhotonIsolated[j][i]->SetXTitle("p_{T} (GeV/c)");
378 fhPtPartonTypeAwayPhotonIsolated[j][i]->SetYTitle("Parton type");
379 fhPtPartonTypeAwayPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
380 fhPtPartonTypeAwayPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
381 fhPtPartonTypeAwayPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
382 outputContainer->Add(fhPtPartonTypeAwayPhotonIsolated[j][i]);
383
384 fhPtPartonTypeAwayPi0Isolated[j][i] = new TH2F(Form("hPtPartonTypeAwayPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
385 Form("Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
386 nptbins,ptmin,ptmax,3,0,3);
387 fhPtPartonTypeAwayPi0Isolated[j][i]->SetXTitle("p_{T} (GeV/c)");
388 fhPtPartonTypeAwayPi0Isolated[j][i]->SetYTitle("Parton type");
389 fhPtPartonTypeAwayPi0Isolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
390 fhPtPartonTypeAwayPi0Isolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
391 fhPtPartonTypeAwayPi0Isolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
392 outputContainer->Add(fhPtPartonTypeAwayPi0Isolated[j][i]);
393
394 // zHard
395
396 fhZHardPhoton[j][i] = new TH2F(Form("hZHardPhoton%s%s",leading[j].Data(),name[i].Data()),
397 Form("Z-Hard of Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
398 nptbins,ptmin,ptmax,200,0,2);
399 fhZHardPhoton[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
400 fhZHardPhoton[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
401 outputContainer->Add(fhZHardPhoton[j][i]);
402
403 fhZHardPi0[j][i] = new TH2F(Form("hZHardPi0%s%s",leading[j].Data(),name[i].Data()),
404 Form("Z-Hard of Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
405 nptbins,ptmin,ptmax,200,0,2);
406 fhZHardPi0[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
407 fhZHardPi0[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
408 outputContainer->Add(fhZHardPi0[j][i]);
409
410 fhZHardPhotonIsolated[j][i] = new TH2F(Form("hZHardPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
411 Form("Z-Hard of Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
412 nptbins,ptmin,ptmax,200,0,2);
413 fhZHardPhotonIsolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
414 fhZHardPhotonIsolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
415 outputContainer->Add(fhZHardPhotonIsolated[j][i]);
416
417 fhZHardPi0Isolated[j][i] = new TH2F(Form("hZHardPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
418 Form("Z-Hard of Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
419 nptbins,ptmin,ptmax,200,0,2);
420 fhZHardPi0Isolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
421 fhZHardPi0Isolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
422 outputContainer->Add(fhZHardPi0Isolated[j][i]);
423
424 // zHard
425
426 fhZPartonPhoton[j][i] = new TH2F(Form("hZPartonPhoton%s%s",leading[j].Data(),name[i].Data()),
427 Form("Z-Parton of Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
428 nptbins,ptmin,ptmax,200,0,2);
429 fhZPartonPhoton[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
430 fhZPartonPhoton[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
431 outputContainer->Add(fhZPartonPhoton[j][i]);
432
433 fhZPartonPi0[j][i] = new TH2F(Form("hZPartonPi0%s%s",leading[j].Data(),name[i].Data()),
434 Form("Z-Parton of Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
435 nptbins,ptmin,ptmax,200,0,2);
436 fhZPartonPi0[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
437 fhZPartonPi0[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
438 outputContainer->Add(fhZPartonPi0[j][i]);
439
440 fhZPartonPhotonIsolated[j][i] = new TH2F(Form("hZPartonPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
441 Form("Z-Parton of Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
442 nptbins,ptmin,ptmax,200,0,2);
443 fhZPartonPhotonIsolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
444 fhZPartonPhotonIsolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
445 outputContainer->Add(fhZPartonPhotonIsolated[j][i]);
446
447 fhZPartonPi0Isolated[j][i] = new TH2F(Form("hZPartonPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
448 Form("Z-Parton of Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
449 nptbins,ptmin,ptmax,200,0,2);
450 fhZPartonPi0Isolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
451 fhZPartonPi0Isolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
452 outputContainer->Add(fhZPartonPi0Isolated[j][i]);
453
454
455 // zJet
456
457 fhZJetPhoton[j][i] = new TH2F(Form("hZJetPhoton%s%s",leading[j].Data(),name[i].Data()),
458 Form("Z-Jet of Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
459 nptbins,ptmin,ptmax,200,0,2);
460 fhZJetPhoton[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
461 fhZJetPhoton[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
462 outputContainer->Add(fhZJetPhoton[j][i]);
463
464 fhZJetPi0[j][i] = new TH2F(Form("hZJetPi0%s%s",leading[j].Data(),name[i].Data()),
465 Form("Z-Jet of Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
466 nptbins,ptmin,ptmax,200,0,2);
467 fhZJetPi0[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
468 fhZJetPi0[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
469 outputContainer->Add(fhZJetPi0[j][i]);
470
471 fhZJetPhotonIsolated[j][i] = new TH2F(Form("hZJetPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
472 Form("Z-Jet of Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
473 nptbins,ptmin,ptmax,200,0,2);
474 fhZJetPhotonIsolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
475 fhZJetPhotonIsolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
476 outputContainer->Add(fhZJetPhotonIsolated[j][i]);
477
478 fhZJetPi0Isolated[j][i] = new TH2F(Form("hZJetPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
479 Form("Z-Jet of Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
480 nptbins,ptmin,ptmax,200,0,2);
481 fhZJetPi0Isolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
482 fhZJetPi0Isolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
483 outputContainer->Add(fhZJetPi0Isolated[j][i]);
484
485
486 // XE
487
488 fhXEPhoton[j][i] = new TH2F(Form("hXEPhoton%s%s",leading[j].Data(),name[i].Data()),
489 Form("Z-Jet of Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
490 nptbins,ptmin,ptmax,200,0,2);
491 fhXEPhoton[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
492 fhXEPhoton[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
493 outputContainer->Add(fhXEPhoton[j][i]);
494
495 fhXEPi0[j][i] = new TH2F(Form("hXEPi0%s%s",leading[j].Data(),name[i].Data()),
496 Form("Z-Jet of Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
497 nptbins,ptmin,ptmax,200,0,2);
498 fhXEPi0[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
499 fhXEPi0[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
500 outputContainer->Add(fhXEPi0[j][i]);
501
502 fhXEPhotonIsolated[j][i] = new TH2F(Form("hXEPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
503 Form("Z-Jet of Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
504 nptbins,ptmin,ptmax,200,0,2);
505 fhXEPhotonIsolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
506 fhXEPhotonIsolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
507 outputContainer->Add(fhXEPhotonIsolated[j][i]);
508
509 fhXEPi0Isolated[j][i] = new TH2F(Form("hXEPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
510 Form("Z-Jet of Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
511 nptbins,ptmin,ptmax,200,0,2);
512 fhXEPi0Isolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
513 fhXEPi0Isolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
514 outputContainer->Add(fhXEPi0Isolated[j][i]);
515
516
517 // XE from UE
518
519 fhXEUEPhoton[j][i] = new TH2F(Form("hXEUEPhoton%s%s",leading[j].Data(),name[i].Data()),
520 Form("Z-Jet of Photon : %s of all particles%s",leading[j].Data(),title[i].Data()),
521 nptbins,ptmin,ptmax,200,0,2);
522 fhXEUEPhoton[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
523 fhXEUEPhoton[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
524 outputContainer->Add(fhXEUEPhoton[j][i]);
525
526 fhXEUEPi0[j][i] = new TH2F(Form("hXEUEPi0%s%s",leading[j].Data(),name[i].Data()),
527 Form("Z-Jet of Pi0 : %s of all particles%s",leading[j].Data(),title[i].Data()),
528 nptbins,ptmin,ptmax,200,0,2);
529 fhXEUEPi0[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
530 fhXEUEPi0[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
531 outputContainer->Add(fhXEUEPi0[j][i]);
532
533 fhXEUEPhotonIsolated[j][i] = new TH2F(Form("hXEUEPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
534 Form("Z-Jet of Photon : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
535 nptbins,ptmin,ptmax,200,0,2);
536 fhXEUEPhotonIsolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
537 fhXEUEPhotonIsolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
538 outputContainer->Add(fhXEUEPhotonIsolated[j][i]);
539
540 fhXEUEPi0Isolated[j][i] = new TH2F(Form("hXEUEPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
541 Form("Z-Jet of Pi0 : %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
7b2086c3 542 nptbins,ptmin,ptmax,200,0,2);
dbb79a05 543 fhXEUEPi0Isolated[j][i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
544 fhXEUEPi0Isolated[j][i]->SetXTitle("p_{T}^{particle} (GeV/c)");
545 outputContainer->Add(fhXEUEPi0Isolated[j][i]);
546 }
7b2086c3 547 }
548
549 return outputContainer;
550
551}
552
553//____________________________________________
554void AliAnaGeneratorKine::GetPartonsAndJets()
555{
556 // Fill data members with partons,jets and generated pt hard
557
558 fStack = GetMCStack() ;
559
560 if(!fStack)
561 {
562 printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - No Stack available, STOP\n");
563 abort();
564 }
565
566 fParton2 = fStack->Particle(2) ;
567 fParton3 = fStack->Particle(3) ;
568 fParton6 = fStack->Particle(6) ;
569 fParton7 = fStack->Particle(7) ;
570
571 Float_t p6phi = fParton6->Phi();
572 if(p6phi < 0) p6phi +=TMath::TwoPi();
573 Float_t p7phi = fParton7->Phi();
574 if(p7phi < 0) p7phi +=TMath::TwoPi();
575
576 //printf("parton6: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton6->Pt(),fParton6->Eta(),p6phi, fParton6->GetPdgCode());
577 //printf("parton7: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton7->Pt(),fParton7->Eta(),p7phi, fParton7->GetPdgCode());
578
579 // Get the jets, only for pythia
580 if(!strcmp(GetReader()->GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
581 {
582 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetReader()->GetGenEventHeader();
583
584 fPtHard = pygeh->GetPtHard();
585
586 //printf("pt Hard %2.2f\n",fPtHard);
587
588 const Int_t nTriggerJets = pygeh->NTriggerJets();
589
590 Float_t tmpjet[]={0,0,0,0};
591
592 // select the closest jet to parton
593 Float_t jet7R = 100;
594 Float_t jet6R = 100;
595
596 for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
597 {
598 pygeh->TriggerJet(ijet, tmpjet);
599
600 TLorentzVector jet(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
601 Float_t jphi = jet.Phi();
602 if(jphi < 0) jphi +=TMath::TwoPi();
603
604 Double_t radius6 = GetIsolationCut()->Radius(fParton6->Eta(), p6phi, jet.Eta() , jphi) ;
605 Double_t radius7 = GetIsolationCut()->Radius(fParton7->Eta(), p7phi, jet.Eta() , jphi) ;
606
607 //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);
608
609 if (radius6 < jet6R)
610 {
611 jet6R = radius6;
612 fJet6 = jet;
613
614 }
615 if (radius7 < jet7R)
616 {
617 jet7R = radius7;
618 fJet7 = jet;
619 }
620
621 } // jet loop
622
623 //printf("jet6: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet6.Pt(),fJet6.Eta(),fJet6.Phi());
624 //printf("jet7: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet7.Pt(),fJet7.Eta(),fJet6.Phi());
625
626 } // pythia header
627
628 fhPtHard ->Fill(fPtHard);
629 fhPtJet ->Fill(fJet6.Pt());
630 fhPtJet ->Fill(fJet7.Pt());
631 fhPtParton ->Fill(fParton6->Pt());
632 fhPtParton ->Fill(fParton7->Pt());
633
634}
635
b94e038e 636//_____________________________________________________
637void AliAnaGeneratorKine::GetXE(TLorentzVector trigger,
638 Int_t indexTrig,
639 Int_t pdgTrig,
640 Bool_t leading[4],
641 Bool_t isolated[4],
642 Int_t iparton)
7b2086c3 643{
644
645 // Calculate the real XE and the UE XE
646
647 Float_t ptThresTrack = 0.2;
648
649 Float_t ptTrig = trigger.Pt();
650 Float_t etaTrig = trigger.Eta();
651 Float_t phiTrig = trigger.Phi();
652 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
653
654 //Loop on primaries, start from position 8, no partons
655 for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
656 {
7b2086c3 657 TParticle * particle = fStack->Particle(ipr) ;
658
659 if(ipr==indexTrig) continue;
660
23fa04c5 661
7b2086c3 662 Int_t pdg = particle->GetPdgCode();
663 Int_t status = particle->GetStatusCode();
664
665 // Compare trigger with final state particles
666 if( status != 1) continue ;
667
668 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
669
670 if(charge==0) continue; // construct xe only with charged
671
672 Float_t pt = particle->Pt();
673 Float_t eta = particle->Eta();
674 Float_t phi = particle->Phi();
675 if(phi < 0 ) phi += TMath::TwoPi();
676
677 if( pt < ptThresTrack) continue ;
678
679 if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
23fa04c5 680
7b2086c3 681 //Isolation
682 Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
683
684 Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
685
686 //Get the index of the mother
687 Int_t ipartonAway = particle->GetFirstMother();
23fa04c5 688 if(ipartonAway < 0) return;
7b2086c3 689 TParticle * mother = fStack->Particle(ipartonAway);
690 while (ipartonAway > 7)
691 {
692 ipartonAway = mother->GetFirstMother();
23fa04c5 693 if(ipartonAway < 0) break;
7b2086c3 694 mother = fStack->Particle(ipartonAway);
695 }
696
697 if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway)
698 {
699 //printf("xE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
700 if(radius > 1 ) continue; // avoid particles too far from trigger
701
702 for( Int_t i = 0; i < 4; i++ )
703 {
704 if(pdgTrig==111)
705 {
dbb79a05 706
707 fhXEPi0[leading[i]][i] ->Fill(ptTrig,xe);
708
709 if(isolated[i])
710 {
711 fhXEPi0Isolated[leading[i]][i] ->Fill(ptTrig,xe);
7b2086c3 712 }
dbb79a05 713
7b2086c3 714 }// pi0
715 else if(pdgTrig==22)
716 {
dbb79a05 717
718 fhXEPhoton[leading[i]][i] ->Fill(ptTrig,xe);
719
720 if(isolated[i])
721 {
722 fhXEPhotonIsolated[leading[i]][i] ->Fill(ptTrig,xe);
723 }
724
7b2086c3 725 } // photon
dbb79a05 726 } // conditions loop
7b2086c3 727 } // Away side
728
dbb79a05 729 if(ipartonAway!=6 && ipartonAway!=7)
7b2086c3 730 {
731
732 //printf("xE UE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
dbb79a05 733
7b2086c3 734 for( Int_t i = 0; i < 4; i++ )
735 {
736 if(pdgTrig==111)
737 {
dbb79a05 738
739 fhXEUEPi0[leading[i]][i] ->Fill(ptTrig,xe);
740
741 if(isolated[i])
742 {
743 fhXEUEPi0Isolated[leading[i]][i] ->Fill(ptTrig,xe);
7b2086c3 744 }
dbb79a05 745
7b2086c3 746 }// pi0
747 else if(pdgTrig==22)
748 {
dbb79a05 749
750 fhXEUEPhoton[leading[i]][i] ->Fill(ptTrig,xe);
751
752 if(isolated[i])
753 {
754 fhXEUEPhotonIsolated[leading[i]][i] ->Fill(ptTrig,xe);
755 }
756
7b2086c3 757 } // photon
758 } // conditions loop
759 } // Away side
760
761 } // primary loop
762
763
764}
765
766//________________________________________
767void AliAnaGeneratorKine::InitParameters()
768{
769
770 //Initialize the parameters of the analysis.
771 AddToHistogramsName("AnaGenKine_");
772
773}
774
b94e038e 775//_____________________________________________________________________
776void AliAnaGeneratorKine::IsLeadingAndIsolated(TLorentzVector trigger,
777 Int_t indexTrig,
778 Int_t pdgTrig,
7b2086c3 779 Bool_t leading[4],
780 Bool_t isolated[4])
781{
dbb79a05 782 // Check if the trigger is the leading particle and if it is isolated
7b2086c3 783 // In case of neutral particles check all neutral or neutral in EMCAL acceptance
784
785 Float_t ptMaxCharged = 0; // all charged
786 Float_t ptMaxNeutral = 0; // all neutral
787 Float_t ptMaxNeutEMCAL = 0; // for neutral, select them in EMCAL acceptance
788 Float_t ptMaxNeutPhot = 0; // for neutral, take only photons
789 Float_t ptMaxNeutEMCALPhot = 0; // for neutral, take only photons in EMCAL acceptance
790
791 leading[0] = 0;
792 leading[1] = 0;
793 leading[2] = 0;
794 leading[3] = 0;
795
796 isolated[0] = 0;
797 isolated[1] = 0;
798 isolated[2] = 0;
799 isolated[3] = 0;
800
801 Float_t ptTrig = trigger.Pt();
802 Float_t etaTrig = trigger.Eta();
803 Float_t phiTrig = trigger.Phi();
804 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
805
806 // Minimum track or cluster energy
807 Float_t ptThresTrack = 0.2;
808 Float_t ptThresCalo = 0.3;
809
810 //Isolation cuts
811 Float_t ptThresIC = 0.5;
812 Float_t rThresIC = 0.4;
813 Int_t nICTrack = 0;
814 Int_t nICNeutral = 0;
815 Int_t nICNeutEMCAL = 0;
816 Int_t nICNeutPhot = 0;
817 Int_t nICNeutEMCALPhot = 0;
818
819 //Loop on primaries, start from position 8, no partons
820 for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
821 {
822 if(ipr == indexTrig) continue;
823 TParticle * particle = fStack->Particle(ipr) ;
824
825 Int_t imother= particle->GetFirstMother();
826 //printf("Leading ipr %d - mother %d\n",ipr, imother);
827
828 if(imother==indexTrig) continue ;
829
830 Int_t pdg = particle->GetPdgCode();
831 Int_t status = particle->GetStatusCode();
832
833 // Compare trigger with final state particles
834 if( status != 1) continue ;
835
836 Float_t pt = particle->Pt();
837 Float_t eta = particle->Eta();
838 Float_t phi = particle->Phi();
839 if(phi < 0 ) phi += TMath::TwoPi();
840
841 if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
842
843 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
844
845 //Isolation
846 Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
847
848 if(charge==0)
849 {
850 if(pt < ptThresCalo) continue ;
851
852 if( ptMaxNeutral < pt ) ptMaxNeutral = pt;
853 if( pt > ptThresIC && radius < rThresIC ) nICNeutral++ ;
854
855 Bool_t phPDG = kFALSE;
856 if(pdg==22 || pdg==111) phPDG = kTRUE;
857
858 //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());
859 if(phPDG)
860 {
861 if( ptMaxNeutPhot < pt) ptMaxNeutPhot = pt;
862 if( pt > ptThresIC && radius < rThresIC ) nICNeutPhot++ ;
863 }
864
865 //EMCAL acceptance
dbb79a05 866 Bool_t inEMCAL = GetFiducialCut()->IsInFiducialCut(trigger,"EMCAL") ;
7b2086c3 867
868 if(inEMCAL)
869 {
870 if( ptMaxNeutEMCAL < pt ) ptMaxNeutEMCAL = pt;
871 if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCAL++ ;
872
873 if(phPDG)
874 {
875 if( ptMaxNeutEMCALPhot < pt ) ptMaxNeutEMCALPhot = pt;
876 if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCALPhot++ ;
877 }
878 }
879 }
880 else
881 {
882 if( pt < ptThresTrack) continue ;
883
884 if( ptMaxCharged < pt ) ptMaxCharged = pt;
885
886 if( pt > ptThresIC && radius < rThresIC )
887 {
888 //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",
889 // ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
890 nICTrack++ ;
891 }
892 }
893
894 } // particle loop
895
896 // Leding decision
897 if(ptTrig > ptMaxCharged)
898 {
899 //printf("pt charged %2.2f, pt neutral %2.2f, pt neutral emcal %2.2f, pt photon %2.2f, pt photon emcal %2.2f\n",
900 // ptMaxCharged, ptMaxNeutral, ptMaxNeutEMCAL,ptMaxNeutPhot, ptMaxNeutEMCALPhot);
901 if(ptTrig > ptMaxNeutral ) leading[0] = kTRUE ;
902 if(ptTrig > ptMaxNeutEMCAL ) leading[1] = kTRUE ;
903 if(ptTrig > ptMaxNeutPhot ) leading[2] = kTRUE ;
904 if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
905 }
906
907 //printf("N in cone over threshold : tracks %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n",
908 // nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
909
910 // Isolation decision
911 if(nICTrack == 0)
912 {
913 if(nICNeutral == 0 ) isolated[0] = kTRUE ;
914 if(nICNeutEMCAL == 0 ) isolated[1] = kTRUE ;
915 if(nICNeutPhot == 0 ) isolated[2] = kTRUE ;
916 if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
917 }
918
919 // Fill histograms if conditions apply for all 4 cases
920 for( Int_t i = 0; i < 4; i++ )
921 {
922 if(pdgTrig==111)
923 {
924 if(leading[i])
925 {
926 fhPtPi0Leading[i]->Fill(ptTrig);
927 if(isolated[i]) fhPtPi0LeadingIsolated[i]->Fill(ptTrig);
928 }
929 }// pi0
930 else if(pdgTrig==22)
931 {
932 if(leading[i])
933 {
934 fhPtPhotonLeading[i]->Fill(ptTrig);
935 if(isolated[i]) fhPtPhotonLeadingIsolated[i]->Fill(ptTrig);
936 }
937 } // photon
938 } // conditions loop
939
940}
941
942//_____________________________________________________
943void AliAnaGeneratorKine::MakeAnalysisFillHistograms()
944{
945 //Particle-Parton Correlation Analysis, fill histograms
946
947 TLorentzVector trigger;
948
949 GetPartonsAndJets();
950
951 for(Int_t ipr = 0; ipr < fStack->GetNprimary(); ipr ++ )
952 {
953 TParticle * particle = fStack->Particle(ipr) ;
954
955 Int_t pdgTrig = particle->GetPdgCode();
956 Int_t statusTrig = particle->GetStatusCode();
957 Int_t imother = particle->GetFirstMother();
958 Float_t ptTrig = particle->Pt();
959
960 // Select final state photons (prompt, fragmentation) or pi0s
961
962 //Check the origin of the photon, accept if prompt or initial/final state radiation
963 if(pdgTrig == 22 && statusTrig == 1)
964 {
965 if(imother > 8) continue;
966 }
967 // If not photon, trigger on pi0
968 else if(pdgTrig != 111) continue;
969
970 // Acceptance and kinematical cuts
2292cf03 971 if( ptTrig < GetMinPt() ) continue ;
7b2086c3 972
973 //EMCAL acceptance, a bit less
974 if(TMath::Abs(particle->Eta()) > 0.6) continue ;
975 if(particle->Phi() > TMath::DegToRad()*176) continue ;
976 if(particle->Phi() < TMath::DegToRad()*74 ) continue ;
dbb79a05 977
978 particle->Momentum(trigger);
7b2086c3 979
dbb79a05 980// Bool_t in = GetFiducialCut()->IsInFiducialCu(trigger,"EMCAL") ;
981// if(! in ) continue ;
982
7b2086c3 983// printf("Particle %d : pdg %d status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f \n",
984// ipr, pdgTrig, statusTrig, imother, ptTrig, particle->Eta(), particle->Phi()*TMath::RadToDeg());
985
986// if(pdgTrig==111)
987// {
988// printf("\t pi0 daughters %d, %d\n", particle->GetDaughter(0), particle->GetDaughter(1));
989// }
7b2086c3 990
991 if (pdgTrig==22 ) fhPtPhoton->Fill(ptTrig);
992 else if(pdgTrig==111) fhPtPi0 ->Fill(ptTrig);
993
994 // Check if it is leading
995 Bool_t leading[4] ;
996 Bool_t isolated[4] ;
997
7b2086c3 998 IsLeadingAndIsolated(trigger, ipr, pdgTrig, leading, isolated);
999
1000 Int_t iparton = -1;
2292cf03 1001 Int_t ok = CorrelateWithPartonOrJet(trigger, ipr, pdgTrig, leading, isolated, iparton);
1002 if(!ok) continue;
1003
7b2086c3 1004 GetXE(trigger,ipr,pdgTrig,leading,isolated,iparton) ;
1005
1006 }
1007
1008 if(GetDebug() > 1) printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - End fill histograms \n");
1009
1010}