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