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