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