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