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