AliPhysics  e59a9ba (e59a9ba)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaGeneratorKine.cxx
Go to the documentation of this file.
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 // --- ROOT system ---
17 #include "TH2F.h"
18 #include "TParticle.h"
19 #include "TDatabasePDG.h"
20 
21 //---- ANALYSIS system ----
22 #include "AliAnaGeneratorKine.h"
23 #include "AliStack.h"
24 #include "AliGenPythiaEventHeader.h"
25 #include "AliAODMCParticle.h"
26 
30 
31 //__________________________________________
33 //__________________________________________
36 fTriggerDetector(), fTriggerDetectorString(),
37 fFidCutTrigger(0),
38 fMinChargedPt(0), fMinNeutralPt(0),
39 fStack(0), fAODMCparticles(0),
40 //fParton2(0), fParton3(0),
41 fParton6(), fParton7(),
42 fParton6PDG(0), fParton7PDG(0),
43 fJet6(), fJet7(),
44 fTrigger(), fLVTmp(),
45 fNPrimaries(0), fPtHard(0),
46 fhPtHard(0), fhPtParton(0), fhPtJet(0),
47 fhPtPartonPtHard(0), fhPtJetPtHard(0), fhPtJetPtParton(0),
48 fhPtOtherDecayMesonId(0),
49 fhPtPi0Status(0), fhPtEtaStatus(0),
50 fhPtPi0DecayStatus(0), fhPtEtaDecayStatus(0), fhPtOtherDecayStatus(0),
51 fhPtPi0Not2Gamma(0), fhPtEtaNot2Gamma(0)
52 {
54 
55  for(Int_t p = 0; p < fgkNmcPrimTypes; p++)
56  {
57  fhPt [p] = 0;
58  fhPtOrigin[p] = 0;
59  fhPtOriginNotFinal[p] = 0;
60  fhPhi [p] = 0;
61  fhEta [p] = 0;
62  fhEtaPhi [p] = 0;
63 
64  fhPhiStatus[p] = 0;
65  fhEtaStatus[p] = 0;
66 
67  for(Int_t i = 0; i < fgkNIso; i++)
68  {
69  fhPtLeading[p][i] = 0;
70  fhPtLeadingSumPt[p][i] = 0;
71  fhPtLeadingIsolated[p][i] = 0;
72  for(Int_t j = 0; j < 2; j++)
73  {
74  fhZHard[p][j][i] = 0;
75  fhZHardIsolated[p][j][i] = 0;
76  fhZParton[p][j][i] = 0;
77  fhZPartonIsolated[p][j][i] = 0;
78  fhZJet[p][j][i] = 0;
79  fhZJetIsolated[p][j][i] = 0;
80  fhXE[p][j][i] = 0;
81  fhXEIsolated[p][j][i] = 0;
82  fhXEUE[p][j][i] = 0;
83  fhXEUEIsolated[p][j][i] = 0;
84 
85  fhPtPartonTypeNear[p][j][i] = 0;
86  fhPtPartonTypeNearIsolated[p][j][i] = 0;
87 
88  fhPtPartonTypeAway[p][j][i] = 0;
89  fhPtPartonTypeAwayIsolated[p][j][i] = 0;
90 
91  if( p == 0 ) fhPtAcceptedGammaJet[j][i] = 0;
92  }
93  }
94  }
95 }
96 
97 //___________________________________________________________________________
99 //___________________________________________________________________________
101  Int_t partType,
102  Bool_t leading [fgkNIso],
103  Bool_t isolated[fgkNIso],
104  Int_t & iparton )
105 {
106  AliDebug(1,"Start");
107 
108  if( fNPrimaries < 7 )
109  {
110  AliDebug(1,Form("End, not enough partons, n primaries %d",fNPrimaries));
111  return kFALSE ;
112  }
113 
114  //
115  // Get the index of the mother
116  //
117  if(fStack) // ESD
118  {
119  iparton = (fStack->Particle(indexTrig))->GetFirstMother();
120  TParticle * mother = fStack->Particle(iparton);
121  while (iparton > 7)
122  {
123  iparton = mother->GetFirstMother();
124  if(iparton < 0)
125  {
126  AliDebug(1,"Negative index, skip ESD event");
127  return kFALSE;
128  }
129  mother = fStack->Particle(iparton);
130  }
131  }
132  else // AOD
133  {
134  iparton = ((AliAODMCParticle*) fAODMCparticles->At(indexTrig))->GetMother();
135  AliAODMCParticle * mother = (AliAODMCParticle*) fAODMCparticles->At(iparton);
136  while (iparton > 7)
137  {
138  iparton = mother->GetMother();
139  if(iparton < 0)
140  {
141  AliDebug(1,"Negative index, skip AOD event");
142  return kFALSE;
143  }
144  mother = (AliAODMCParticle*) fAODMCparticles->At(iparton);
145  }
146  }
147 
148  //printf("Mother is parton %d with pdg %d\n",imom,mother->GetPdgCode());
149 
150  if(iparton < 6)
151  {
152  AliDebug(1,Form("This particle is not from hard process - pdg %d, parton index %d\n",partType, iparton));
153  return kFALSE;
154  }
155 
156  //
157  // Get the kinematics
158  //
159  Float_t ptTrig = fTrigger.Pt();
160  Float_t phiTrig = fTrigger.Phi();
161  Float_t etaTrig = fTrigger.Eta();
162 
163  AliDebug(1,Form("Trigger pdg %d pT %2.2f, phi %2.2f, eta %2.2f", partType,ptTrig,phiTrig*TMath::RadToDeg(),etaTrig));
164 
165  Float_t jetPt = fJet6.Pt();
166  Float_t jetPhi = fJet6.Phi();
167  Float_t jetEta = fJet6.Eta();
168 
169  AliDebug(1,Form("Jet 6 pT %2.2f, phi %2.2f, eta %2.2f",jetPt,jetPhi*TMath::RadToDeg(),jetEta));
170 
171  Float_t awayJetPt = fJet7.Pt();
172  Float_t awayJetEta = fJet7.Eta();
173  Float_t awayJetPhi = fJet7.Phi();
174 
175  AliDebug(1,Form("Jet 7 pT %2.2f, phi %2.2f, eta %2.2f",awayJetPt,awayJetPhi*TMath::RadToDeg(),awayJetEta));
176 
177  Float_t partonPt = fParton6.Pt();
178 
179  Int_t nearPDG = fParton6PDG;
180  Int_t awayPDG = fParton7PDG;
181 
182  AliDebug(1,Form("Parton6 pT pT %2.2f, pdg %d",fParton6.Pt(),fParton6PDG));
183  AliDebug(1,Form("Parton7 pT pT %2.2f, pdg %d",fParton7.Pt(),fParton7PDG));
184 
185  if ( iparton == 7 )
186  {
187  partonPt = fParton7.Pt();
188 
189  jetPt = fJet7.Pt();
190  jetPhi = fJet7.Phi();
191  jetEta = fJet7.Eta();
192 
193  awayJetPt = fJet6.Pt();
194  awayJetEta = fJet6.Eta();
195  awayJetPhi = fJet6.Phi();
196 
197  nearPDG = fParton7PDG;
198  awayPDG = fParton6PDG;
199  }
200 
201  Float_t deltaPhi = TMath::Abs(phiTrig-awayJetPhi) *TMath::RadToDeg();
202  AliDebug(1,Form("Trigger Away jet phi %2.2f\n",deltaPhi));
203 
204  //
205  // Get id of parton in near and away side
206  //
207  Int_t away = -1;
208  Int_t near = -1;
209  if (nearPDG == 22) near = 0;
210  else if(nearPDG == 21) near = 1;
211  else near = 2;
212 
213  if (awayPDG == 22) away = 0;
214  else if(awayPDG == 21) away = 1;
215  else away = 2;
216 
217  // RATIOS
218 
219  Float_t zHard = -1;
220  Float_t zPart = -1;
221  Float_t zJet = -1;
222 
223  if( fPtHard > 0 ) zHard = ptTrig / fPtHard ;
224  if( partonPt > 0 ) zPart = ptTrig / partonPt;
225  if( jetPt > 0 ) zJet = ptTrig / jetPt ;
226 
227  //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard);
228 
229  //printf("Z: hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
230 
231  // conditions loop
232  for( Int_t i = 0; i < fgkNIso ; i++ )
233  {
234  fhPtPartonTypeNear[partType][leading[i]][i]->Fill(ptTrig, near, GetEventWeight());
235  fhPtPartonTypeAway[partType][leading[i]][i]->Fill(ptTrig, away, GetEventWeight());
236 
237  fhZHard [partType][leading[i]][i]->Fill(ptTrig, zHard, GetEventWeight());
238  fhZParton[partType][leading[i]][i]->Fill(ptTrig, zPart, GetEventWeight());
239  fhZJet [partType][leading[i]][i]->Fill(ptTrig, zJet , GetEventWeight());
240 
241  if(isolated[i])
242  {
243  fhPtPartonTypeNearIsolated[partType][leading[i]][i]->Fill(ptTrig, near, GetEventWeight());
244  fhPtPartonTypeAwayIsolated[partType][leading[i]][i]->Fill(ptTrig, away, GetEventWeight());
245 
246  fhZHardIsolated [partType][leading[i]][i]->Fill(ptTrig, zHard, GetEventWeight());
247  fhZPartonIsolated[partType][leading[i]][i]->Fill(ptTrig, zPart, GetEventWeight());
248  fhZJetIsolated [partType][leading[i]][i]->Fill(ptTrig, zJet , GetEventWeight());
249  }
250 
251  if(partType == kmcPrimPhoton && deltaPhi < 220 && deltaPhi > 140 && TMath::Abs(awayJetEta) < 0.6)
252  {
253  //printf("Accept jet\n");
254  fhPtAcceptedGammaJet[leading[i]][i]->Fill(ptTrig, away, GetEventWeight());
255  }
256  //else printf("Reject jet!!!\n");
257 
258  } // conditions loop
259 
260  AliDebug(1,"End TRUE");
261 
262  return kTRUE;
263 }
264 
265 //____________________________________________________
267 //____________________________________________________
269 {
270  TList * outputContainer = new TList() ;
271  outputContainer->SetName("GenKineHistos") ;
272 
274  Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
275  Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
276 
277  Int_t nptsumbins = GetHistogramRanges()->GetHistoNPtSumBins();
278  Float_t ptsummax = GetHistogramRanges()->GetHistoPtSumMax();
279  Float_t ptsummin = GetHistogramRanges()->GetHistoPtSumMin();
280 
281  fhPtHard = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax);
282  fhPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
283  outputContainer->Add(fhPtHard);
284 
285  fhPtParton = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax);
286  fhPtParton->SetXTitle("#it{p}_{T}^{parton} (GeV/#it{c})");
287  outputContainer->Add(fhPtParton);
288 
289  fhPtJet = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax);
290  fhPtJet->SetXTitle("#it{p}_{T}^{jet} (GeV/#it{c})");
291  outputContainer->Add(fhPtJet);
292 
293  fhPtPartonPtHard = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
294  fhPtPartonPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
295  fhPtPartonPtHard->SetYTitle("#it{p}_{T}^{parton}/#it{p}_{T}^{hard}");
296  outputContainer->Add(fhPtPartonPtHard);
297 
298  fhPtJetPtHard = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
299  fhPtJetPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
300  fhPtJetPtHard->SetYTitle("#it{p}_{T}^{jet}/#it{p}_{T}^{hard}");
301  outputContainer->Add(fhPtJetPtHard);
302 
303  fhPtJetPtParton = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
304  fhPtJetPtParton->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
305  fhPtJetPtParton->SetYTitle("#it{p}_{T}^{jet}/#it{p}_{T}^{parton}");
306  outputContainer->Add(fhPtJetPtParton);
307 
308  fhPtPi0Not2Gamma = new TH1F("hPtPi0Not2Gamma","#pi^{0} decay other than 2 #gamma",nptbins,ptmin,ptmax);
309  fhPtPi0Not2Gamma->SetXTitle("#it{p}_{T} (GeV/#it{c})");
310  outputContainer->Add(fhPtPi0Not2Gamma);
311 
312  fhPtEtaNot2Gamma = new TH1F("hPtEtaNot2Gamma","#eta decay other than 2 #gamma",nptbins,ptmin,ptmax);
313  fhPtEtaNot2Gamma->SetXTitle("#it{p}_{T} (GeV/#it{c})");
314  outputContainer->Add(fhPtEtaNot2Gamma);
315 
316  fhPtGammaFromPi0Not2Gamma = new TH1F("hPtGammaFromPi0Not2Gamma","#gamma from #pi^{0} decay other than 2 #gamma",nptbins,ptmin,ptmax);
317  fhPtGammaFromPi0Not2Gamma->SetXTitle("#it{p}_{T} (GeV/#it{c})");
318  outputContainer->Add(fhPtGammaFromPi0Not2Gamma);
319 
320  fhPtGammaFromEtaNot2Gamma = new TH1F("hPtGammaFromEtaNot2Gamma","#gamma from #eta decay other than 2 #gamma",nptbins,ptmin,ptmax);
321  fhPtGammaFromEtaNot2Gamma->SetXTitle("#it{p}_{T} (GeV/#it{c})");
322  outputContainer->Add(fhPtGammaFromEtaNot2Gamma);
323 
324  fhPtPi0Status = new TH2F("hPtPi0Status","#pi^{0} status",nptbins,ptmin,ptmax,101,-50,50);
325  fhPtPi0Status->SetXTitle("#it{p}_{T} (GeV/#it{c})");
326  fhPtPi0Status->SetYTitle("status");
327  outputContainer->Add(fhPtPi0Status);
328 
329  fhPtEtaStatus = new TH2F("hPtEtaStatus","#eta status",nptbins,ptmin,ptmax,101,-50,50);
330  fhPtEtaStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
331  fhPtEtaStatus->SetYTitle("status");
332  outputContainer->Add(fhPtEtaStatus);
333 
334  fhPtPi0DecayStatus = new TH2F("hPtPi0DecayStatus","#gamma from #pi^{0}, mother status",nptbins,ptmin,ptmax,101,-50,50);
335  fhPtPi0DecayStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
336  fhPtPi0DecayStatus->SetYTitle("status");
337  outputContainer->Add(fhPtPi0DecayStatus);
338 
339  fhPtEtaDecayStatus = new TH2F("hPtEtaStatus","#gamma from #eta, mother status",nptbins,ptmin,ptmax,101,-50,50);
340  fhPtEtaDecayStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
341  fhPtEtaDecayStatus->SetYTitle("status");
342  outputContainer->Add(fhPtEtaDecayStatus);
343 
344  fhPtOtherDecayStatus = new TH2F("hPtOtherDecayStatus","#gamma from other decay particle, mother status",nptbins,ptmin,ptmax,101,-50,50);
345  fhPtOtherDecayStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
346  fhPtOtherDecayStatus->SetYTitle("status");
347  outputContainer->Add(fhPtOtherDecayStatus);
348 
349  fhPtOtherDecayMesonId = new TH2F("hPtOtherDecayMesonId","Particle decaying into #gamma, Id",nptbins,ptmin,ptmax,8,0,8) ;
350  fhPtOtherDecayMesonId->SetXTitle("#it{p}_{T} (GeV/#it{c})");
351  fhPtOtherDecayMesonId->SetYTitle("Origin");
352  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(1 ,"Resonances");
353  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(2 ,"#rho");
354  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(3 ,"#omega");
355  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(4 ,"K");
356  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(5 ,"Other");
357  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(6 ,"#eta");
358  fhPtOtherDecayMesonId->GetYaxis()->SetBinLabel(7 ,"#eta prime");
359  outputContainer->Add(fhPtOtherDecayMesonId) ;
360 
361  TString name [] = {"","_EMC","_Photon","_EMC_Photon"};
362  TString title [] = {"",", neutral in EMCal",", neutral only #gamma-like",", neutral in EMCal and only #gamma-like"};
363  TString leading[] = {"NotLeading","Leading"};
364 
365  TString partTitl[] = {"#gamma_{direct}","#gamma_{decay}^{#pi}","#gamma_{decay}^{#eta}","#gamma_{decay}^{other}","#pi^{0}","#eta"};
366  TString particle[] = {"DirectPhoton" ,"Pi0DecayPhoton" ,"EtaDecayPhoton" ,"OtherDecayPhoton" ,"Pi0" ,"Eta"};
367 
368  for(Int_t p = 0; p < fgkNmcPrimTypes; p++)
369  {
370  fhPt[p] = new TH1F(Form("h%sPt",particle[p].Data()),Form("Input %s p_{T}",partTitl[p].Data()),nptbins,ptmin,ptmax);
371  fhPt[p]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
372  outputContainer->Add(fhPt[p]);
373 
374  fhPhi[p] = new TH1F(Form("h%sPhi",particle[p].Data()),Form("Input %s #phi with #it{p}_{T} > %f GeV/c",partTitl[p].Data(),GetMinPt()),180,0,TMath::TwoPi());
375  fhPhi[p]->SetXTitle("#phi (rad)");
376  outputContainer->Add(fhPhi[p]);
377 
378  fhEta[p] = new TH1F(Form("h%sEta",particle[p].Data()),Form("Input %s #eta with #it{p}_{T} > %f GeV/c",partTitl[p].Data(),GetMinPt()),200,-2,2);
379  fhEta[p]->SetXTitle("#eta");
380  outputContainer->Add(fhEta[p]);
381 
382  fhEtaPhi[p] = new TH2F(Form("h%sEtaPhi",particle[p].Data()),
383  Form("Input %s #eta vs #phi with #it{p}_{T} > %f GeV/c",partTitl[p].Data(),GetMinPt()),
384  200,-2,2,180,0,TMath::TwoPi());
385  fhEtaPhi[p]->SetXTitle("#eta");
386  fhEtaPhi[p]->SetYTitle("#phi (rad)");
387  outputContainer->Add(fhEtaPhi[p]);
388 
389  fhPhiStatus[p] = new TH2F(Form("h%sPhiStatus",particle[p].Data()),
390  Form("Input %s #phi vs status code with #it{p}_{T} > %f GeV/c",partTitl[p].Data(),GetMinPt()),
391  180,0,TMath::TwoPi(),101,-50,50);
392  fhPhiStatus[p]->SetXTitle("#phi (rad)");
393  fhPhiStatus[p]->SetYTitle("status code");
394  outputContainer->Add(fhPhiStatus[p]);
395 
396  fhEtaStatus[p] = new TH2F(Form("h%sEtaStatus",particle[p].Data()),
397  Form("Input %s #eta vs status code with #it{p}_{T} > %f GeV/c",partTitl[p].Data(),GetMinPt()),
398  200,-2,2,101,-50,50);
399  fhEtaStatus[p]->SetXTitle("#eta");
400  fhEtaStatus[p]->SetYTitle("status code");
401  outputContainer->Add(fhEtaStatus[p]);
402 
403 
404  fhPtOrigin[p] = new TH2F(Form("h%sPtOrigin",particle[p].Data()),Form("Input %s p_{T} vs origin",partTitl[p].Data()),nptbins,ptmin,ptmax,11,0,11) ;
405  fhPtOrigin[p]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
406  fhPtOrigin[p]->SetYTitle("Origin");
407  fhPtOrigin[p]->GetYaxis()->SetBinLabel(1 ,"Status 21");
408  fhPtOrigin[p]->GetYaxis()->SetBinLabel(2 ,"Quark");
409  fhPtOrigin[p]->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
410  fhPtOrigin[p]->GetYaxis()->SetBinLabel(4 ,"Resonances");
411  fhPtOrigin[p]->GetYaxis()->SetBinLabel(5 ,"#rho");
412  fhPtOrigin[p]->GetYaxis()->SetBinLabel(6 ,"#omega");
413  fhPtOrigin[p]->GetYaxis()->SetBinLabel(7 ,"K");
414  fhPtOrigin[p]->GetYaxis()->SetBinLabel(8 ,"Other");
415  fhPtOrigin[p]->GetYaxis()->SetBinLabel(9 ,"#eta");
416  fhPtOrigin[p]->GetYaxis()->SetBinLabel(10 ,"#eta prime");
417  outputContainer->Add(fhPtOrigin[p]) ;
418 
419  fhPtOriginNotFinal[p] = new TH2F(Form("h%sPtOriginNotFinal",particle[p].Data()),Form("Input %s p_{T} vs origin, status 0",partTitl[p].Data()),nptbins,ptmin,ptmax,11,0,11) ;
420  fhPtOriginNotFinal[p]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
421  fhPtOriginNotFinal[p]->SetYTitle("Origin");
422  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(1 ,"Status 21");
423  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(2 ,"Quark");
424  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
425  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(4 ,"Resonances");
426  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(5 ,"#rho");
427  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(6 ,"#omega");
428  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(7 ,"K");
429  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(8 ,"Other");
430  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(9 ,"#eta");
431  fhPtOriginNotFinal[p]->GetYaxis()->SetBinLabel(10 ,"#eta prime");
432  outputContainer->Add(fhPtOriginNotFinal[p]) ;
433 
434  for(Int_t i = 0; i < fgkNIso; i++)
435  {
436  // Pt
437 
438  fhPtLeading[p][i] = new TH1F(Form("h%sPtLeading%s", particle[p].Data(), name[i].Data()),
439  Form("%s: Leading of all particles%s",partTitl[p].Data(),title[i].Data()),
440  nptbins,ptmin,ptmax);
441  fhPtLeading[p][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
442  outputContainer->Add(fhPtLeading[p][i]);
443 
444  fhPtLeadingIsolated[p][i] = new TH1F(Form("h%sPtLeadingIsolated%s", particle[p].Data(), name[i].Data()),
445  Form("%s: Leading of all particles%s, isolated",partTitl[p].Data(),title[i].Data()),
446  nptbins,ptmin,ptmax);
447  fhPtLeadingIsolated[p][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
448  outputContainer->Add(fhPtLeadingIsolated[p][i]);
449 
450  fhPtLeadingSumPt[p][i] = new TH2F(Form("h%sPtLeadingSumPt%s", particle[p].Data(), name[i].Data()),
451  Form("%s: Leading of all particles%s",partTitl[p].Data(),title[i].Data()),
452  nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
453  fhPtLeadingSumPt[p][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
454  fhPtLeadingSumPt[p][i]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
455  outputContainer->Add(fhPtLeadingSumPt[p][i]);
456 
457 
458  // Leading or not loop
459  for(Int_t j = 0; j < fgkNLead; j++)
460  {
461  if(p==0)
462  {
463  fhPtAcceptedGammaJet[j][i] = new TH2F(Form("hPtAcceptedGammaJet%s%s", leading[j].Data(), name[i].Data()),
464  Form("#gamma-jet: %s of all particles%s", leading[j].Data(), title[i].Data()),
465  nptbins,ptmin,ptmax,3,0,3);
466  fhPtAcceptedGammaJet[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
467  fhPtAcceptedGammaJet[j][i]->SetYTitle("Parton type");
468  fhPtAcceptedGammaJet[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
469  fhPtAcceptedGammaJet[j][i]->GetYaxis()->SetBinLabel(2,"g");
470  fhPtAcceptedGammaJet[j][i]->GetYaxis()->SetBinLabel(3,"q");
471  outputContainer->Add(fhPtAcceptedGammaJet[j][i]);
472  }
473  // Near side parton
474 
475  fhPtPartonTypeNear[p][j][i] = new TH2F(Form("h%sPtPartonTypeNear%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
476  Form("%s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
477  nptbins,ptmin,ptmax,3,0,3);
478  fhPtPartonTypeNear[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
479  fhPtPartonTypeNear[p][j][i]->SetYTitle("Parton type");
480  fhPtPartonTypeNear[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
481  fhPtPartonTypeNear[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
482  fhPtPartonTypeNear[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
483  outputContainer->Add(fhPtPartonTypeNear[p][j][i]);
484 
485  fhPtPartonTypeNearIsolated[p][j][i] = new TH2F(Form("h%sPtPartonTypeNear%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
486  Form("%s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
487  nptbins,ptmin,ptmax,3,0,3);
488  fhPtPartonTypeNearIsolated[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
489  fhPtPartonTypeNearIsolated[p][j][i]->SetYTitle("Parton type");
490  fhPtPartonTypeNearIsolated[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
491  fhPtPartonTypeNearIsolated[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
492  fhPtPartonTypeNearIsolated[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
493  outputContainer->Add(fhPtPartonTypeNearIsolated[p][j][i]);
494 
495 
496  // Away side parton
497 
498  fhPtPartonTypeAway[p][j][i] = new TH2F(Form("h%sPtPartonTypeAway%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
499  Form("%s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
500  nptbins,ptmin,ptmax,3,0,3);
501  fhPtPartonTypeAway[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
502  fhPtPartonTypeAway[p][j][i]->SetYTitle("Parton type");
503  fhPtPartonTypeAway[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
504  fhPtPartonTypeAway[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
505  fhPtPartonTypeAway[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
506  outputContainer->Add(fhPtPartonTypeAway[p][j][i]);
507 
508  fhPtPartonTypeAwayIsolated[p][j][i] = new TH2F(Form("h%sPtPartonTypeAway%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
509  Form("%s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
510  nptbins,ptmin,ptmax,3,0,3);
511  fhPtPartonTypeAwayIsolated[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
512  fhPtPartonTypeAwayIsolated[p][j][i]->SetYTitle("Parton type");
513  fhPtPartonTypeAwayIsolated[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
514  fhPtPartonTypeAwayIsolated[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
515  fhPtPartonTypeAwayIsolated[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
516  outputContainer->Add(fhPtPartonTypeAwayIsolated[p][j][i]);
517 
518  // zHard
519 
520  fhZHard[p][j][i] = new TH2F(Form("h%sZHard%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
521  Form("#it{z}_{Hard} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
522  nptbins,ptmin,ptmax,200,0,2);
523  fhZHard[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
524  fhZHard[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
525  outputContainer->Add(fhZHard[p][j][i]);
526 
527  fhZHardIsolated[p][j][i] = new TH2F(Form("h%sZHard%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
528  Form("#it{z}_{Hard} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
529  nptbins,ptmin,ptmax,200,0,2);
530  fhZHardIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
531  fhZHardIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
532  outputContainer->Add(fhZHardIsolated[p][j][i]);
533 
534  // zHard
535 
536  fhZParton[p][j][i] = new TH2F(Form("h%sZParton%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
537  Form("#it{z}_{Parton} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
538  nptbins,ptmin,ptmax,200,0,2);
539  fhZParton[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
540  fhZParton[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
541  outputContainer->Add(fhZParton[p][j][i]);
542 
543  fhZPartonIsolated[p][j][i] = new TH2F(Form("h%sZParton%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
544  Form("#it{z}_{Parton} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
545  nptbins,ptmin,ptmax,200,0,2);
546  fhZPartonIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
547  fhZPartonIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
548  outputContainer->Add(fhZPartonIsolated[p][j][i]);
549 
550 
551  // zJet
552 
553  fhZJet[p][j][i] = new TH2F(Form("h%sZJet%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
554  Form("#it{z}_{Jet} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
555  nptbins,ptmin,ptmax,200,0,2);
556  fhZJet[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
557  fhZJet[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
558  outputContainer->Add(fhZJet[p][j][i]);
559 
560  fhZJetIsolated[p][j][i] = new TH2F(Form("h%sZJet%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
561  Form("#it{z}_{Jet} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
562  nptbins,ptmin,ptmax,200,0,2);
563  fhZJetIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
564  fhZJetIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
565  outputContainer->Add(fhZJetIsolated[p][j][i]);
566 
567 
568  // XE
569 
570  fhXE[p][j][i] = new TH2F(Form("h%sXE%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
571  Form("#it{z}_{Jet} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
572  nptbins,ptmin,ptmax,200,0,2);
573  fhXE[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
574  fhXE[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
575  outputContainer->Add(fhXE[p][j][i]);
576 
577  fhXEIsolated[p][j][i] = new TH2F(Form("h%sXE%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
578  Form("#it{z}_{Jet} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
579  nptbins,ptmin,ptmax,200,0,2);
580  fhXEIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
581  fhXEIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
582  outputContainer->Add(fhXEIsolated[p][j][i]);
583 
584 
585  // XE from UE
586 
587  fhXEUE[p][j][i] = new TH2F(Form("h%sXEUE%s%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
588  Form("#it{z}_{Jet} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
589  nptbins,ptmin,ptmax,200,0,2);
590  fhXEUE[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
591  fhXEUE[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
592  outputContainer->Add(fhXEUE[p][j][i]);
593 
594  fhXEUEIsolated[p][j][i] = new TH2F(Form("h%sXEUE%sIsolated%s", particle[p].Data(), leading[j].Data(), name[i].Data()),
595  Form("#it{z}_{Jet} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
596  nptbins,ptmin,ptmax,200,0,2);
597  fhXEUEIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
598  fhXEUEIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
599  outputContainer->Add(fhXEUEIsolated[p][j][i]);
600  }
601  }
602  }
603 
604  return outputContainer;
605 }
606 
607 //____________________________________________
609 //____________________________________________
611 {
612  AliDebug(1,"Start");
613 
614 // if( nPrimary > 2 ) fParton2 = fStack->Particle(2) ;
615 // if( nPrimary > 3 ) fParton3 = fStack->Particle(3) ;
616 
617  Float_t p6phi = -1 ;
618  Float_t p6eta = -10;
619  Float_t p6pt = 0 ;
620 
621  if( fNPrimaries > 6 )
622  {
623  p6pt = fParton6.Pt();
624  p6eta = fParton6.Eta();
625  p6phi = fParton6.Phi();
626  if(p6phi < 0) p6phi +=TMath::TwoPi();
627  }
628 
629  Float_t p7phi = -1 ;
630  Float_t p7eta = -10;
631  Float_t p7pt = 0 ;
632 
633  if( fNPrimaries > 7 )
634  {
635  p7pt = fParton7.Pt();
636  p7phi = fParton7.Eta();
637  p7phi = fParton7.Phi();
638  if(p7phi < 0) p7phi +=TMath::TwoPi();
639  }
640 
641  //printf("parton6: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",p6pt,p6eta,p6phi, fParton6PDG);
642  //printf("parton7: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",p7pt,p7eta,p7phi, fParton7PDG);
643 
644  // Get the jets, only for pythia
645  if(!strcmp(GetReader()->GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
646  {
647  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetReader()->GetGenEventHeader();
648 
649  fPtHard = pygeh->GetPtHard();
650 
651  //printf("pt Hard %2.2f\n",fPtHard);
652 
653  const Int_t nTriggerJets = pygeh->NTriggerJets();
654 
655  Float_t tmpjet[]={0,0,0,0};
656 
657  // select the closest jet to parton
658  Float_t jet7R = 100;
659  Float_t jet6R = 100;
660 
661  for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
662  {
663  pygeh->TriggerJet(ijet, tmpjet);
664 
665  fLVTmp.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
666  Float_t jphi = fLVTmp.Phi();
667  if(jphi < 0) jphi +=TMath::TwoPi();
668 
669  Double_t radius6 = GetIsolationCut()->Radius(p6eta, p6phi, fLVTmp.Eta() , jphi) ;
670  Double_t radius7 = GetIsolationCut()->Radius(p7eta, p7phi, fLVTmp.Eta() , jphi) ;
671 
672  //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);
673 
674  if (radius6 < jet6R)
675  {
676  jet6R = radius6;
677  fJet6 = fLVTmp;
678 
679  }
680 
681  if (radius7 < jet7R)
682  {
683  jet7R = radius7;
684  fJet7 = fLVTmp;
685  }
686  } // jet loop
687 
688  //printf("jet6: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet6.Pt(),fJet6.Eta(),fJet6.Phi());
689  //printf("jet7: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet7.Pt(),fJet7.Eta(),fJet6.Phi());
690  } // pythia header
691 
692  fhPtHard ->Fill(fPtHard , GetEventWeight());
693  fhPtJet ->Fill(fJet6.Pt(), GetEventWeight());
694  fhPtJet ->Fill(fJet7.Pt(), GetEventWeight());
695  fhPtParton ->Fill(p6pt , GetEventWeight());
696  fhPtParton ->Fill(p7pt , GetEventWeight());
697 
698  if( fPtHard > 0 )
699  {
704  }
705 
706  if( p6pt > 0 ) fhPtJetPtParton ->Fill(fPtHard, fJet6.Pt()/p6pt, GetEventWeight());
707  if( p7pt > 0 ) fhPtJetPtParton ->Fill(fPtHard, fJet7.Pt()/p7pt, GetEventWeight());
708 
709  AliDebug(1,"End");
710 }
711 
712 //_____________________________________________________
714 //_____________________________________________________
715 void AliAnaGeneratorKine::GetXE(Int_t indexTrig,
716  Int_t partType,
717  Bool_t leading [fgkNIso],
718  Bool_t isolated[fgkNIso],
719  Int_t iparton)
720 {
721  AliDebug(1,"Start");
722 
723  Float_t ptTrig = fTrigger.Pt();
724  Float_t phiTrig = fTrigger.Phi();
725  if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
726 
727  Int_t pdg = 0;
728  Int_t status = 0;
729  Int_t ipartonAway = 0;
730  Int_t charge = 0;
731  //Loop on primaries, start from position 8, no partons
732  for(Int_t ipr = 8; ipr < fNPrimaries; ipr ++ )
733  {
734  if(ipr==indexTrig) continue;
735 
736  // Get ESD particle kinematics
737  if(fStack)
738  {
739  TParticle * particle = fStack->Particle(ipr) ;
740 
741  pdg = particle->GetPdgCode();
742  status = particle->GetStatusCode();
743 
744  // Compare trigger with final state particles
745  if( status != 1) continue ; // do it here to avoid crashes
746 
747  particle->Momentum(fLVTmp);
748 
749  charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
750  }
751  else // AOD particle kinematics
752  {
753  AliAODMCParticle * particle = (AliAODMCParticle*) fAODMCparticles->At(ipr) ;
754 
755  pdg = particle->GetPdgCode();
756  status = particle->GetStatus();
757 
758  // Compare trigger with final state particles
759  if( status != 1) continue ; // do it here to avoid crashes
760 
761  fLVTmp.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->E());
762 
763  charge = particle->Charge();
764  }
765 
766  // construct xe only with charged
767  if( charge == 0 ) continue;
768 
769  Float_t pt = fLVTmp.Pt();
770  Float_t phi = fLVTmp.Phi();
771  if(phi < 0 ) phi += TMath::TwoPi();
772 
773  if( pt < fMinChargedPt) continue ;
774 
775  Bool_t inTPC = GetFiducialCut()->IsInFiducialCut(fLVTmp.Eta(),fLVTmp.Phi(),kCTS) ;
776 
777  if(!inTPC) continue;
778 
779  // ---------------------------------------------------
780  // Get the index of the mother, get from what parton
781  // ESD
782  if(fStack)
783  {
784  ipartonAway = fStack->Particle(ipr)->GetFirstMother();
785  if(ipartonAway < 0)
786  {
787  AliDebug(1,"End, no mother index");
788  return;
789  }
790 
791  TParticle * mother = fStack->Particle(ipartonAway);
792  while (ipartonAway > 7)
793  {
794  ipartonAway = mother->GetFirstMother();
795  if(ipartonAway < 0) break;
796  mother = fStack->Particle(ipartonAway);
797  }
798  }
799  else // AOD
800  {
801  ipartonAway = ((AliAODMCParticle*) fAODMCparticles->At(ipr))->GetMother();
802  if(ipartonAway < 0)
803  {
804  AliDebug(1,"End, no mother index");
805  return;
806  }
807 
808  AliAODMCParticle * mother = (AliAODMCParticle*) fAODMCparticles->At(ipartonAway);
809  while (ipartonAway > 7)
810  {
811  ipartonAway = mother->GetMother();
812  if(ipartonAway < 0) break;
813  mother = (AliAODMCParticle*) fAODMCparticles->At(ipartonAway);
814  }
815  }
816 
817  //-----------------------------------------
818  // Get XE of particles belonging to the jet
819  // on the opposite side of the trigger
820 
821  Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
822 
823  if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway)
824  {
825  for( Int_t i = 0; i < fgkNIso; i++ )
826  {
827  fhXE[partType][leading[i]][i] ->Fill(ptTrig, xe, GetEventWeight());
828 
829  if(isolated[i])
830  {
831  fhXEIsolated[partType][leading[i]][i] ->Fill(ptTrig, xe, GetEventWeight());
832  }
833  } // conditions loop
834  } // Away side
835 
836  //----------------------------------------------------------
837  // Get the XE from particles not attached to any of the jets
838  if(ipartonAway!=6 && ipartonAway!=7)
839  {
840  for( Int_t i = 0; i < fgkNIso; i++ )
841  {
842  fhXEUE[partType][leading[i]][i] ->Fill(ptTrig, xe, GetEventWeight());
843 
844  if(isolated[i])
845  {
846  fhXEUEIsolated[partType][leading[i]][i] ->Fill(ptTrig, xe, GetEventWeight());
847  }
848  } // conditions loop
849  } // Away side
850 
851  } // primary loop
852 
853  AliDebug(1,"End");
854 }
855 
856 //________________________________________
858 //________________________________________
860 {
861  AddToHistogramsName("AnaGenKine_");
862 
864 
865  fMinChargedPt = 0.2;
866  fMinNeutralPt = 0.3;
867 }
868 
869 //_____________________________________________________________________
872 //_____________________________________________________________________
874  Int_t partType,
875  Bool_t leading[fgkNIso],
876  Bool_t isolated[fgkNIso])
877 {
878  AliDebug(1,"Start");
879 
880  Float_t ptMaxCharged = 0; // all charged
881  Float_t ptMaxNeutral = 0; // all neutral
882  Float_t ptMaxNeutEMCAL = 0; // for neutral, select them in EMCAL acceptance
883  Float_t ptMaxNeutPhot = 0; // for neutral, take only photons
884  Float_t ptMaxNeutEMCALPhot = 0; // for neutral, take only photons in EMCAL acceptance
885 
886  leading[0] = 0;
887  leading[1] = 0;
888  leading[2] = 0;
889  leading[3] = 0;
890 
891  isolated[0] = 0;
892  isolated[1] = 0;
893  isolated[2] = 0;
894  isolated[3] = 0;
895 
896  Float_t ptTrig = fTrigger.Pt();
897  Float_t etaTrig = fTrigger.Eta();
898  Float_t phiTrig = fTrigger.Phi();
899  if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
900 
901  // Minimum track or cluster energy
902 
903  //Isolation cuts
904  Float_t ptThresIC = GetIsolationCut()->GetPtThreshold();
905  Float_t sumThresIC = GetIsolationCut()->GetPtThreshold();
906  Float_t rThresIC = GetIsolationCut()->GetConeSize();
907  Float_t isoMethod = GetIsolationCut()->GetICMethod();
908 
909  // Counters
910  Int_t nICTrack = 0;
911  Int_t nICNeutral = 0;
912  Int_t nICNeutEMCAL = 0;
913  Int_t nICNeutPhot = 0;
914  Int_t nICNeutEMCALPhot = 0;
915 
916  // Sum of pT
917  Float_t sumNePt = 0;
918  Float_t sumNePtPhot = 0;
919  Float_t sumNePtEMC = 0;
920  Float_t sumNePtEMCPhot = 0;
921  Float_t sumChPt = 0;
922 
923  // Loop on primaries, start from position 8, no partons
924 
925  Int_t imother = -1;
926  Int_t pdg = 0;
927  Int_t status = 0;
928  Int_t charge = 0;
929  for(Int_t ipr = 8; ipr < fNPrimaries; ipr ++ )
930  {
931  if(ipr == indexTrig) continue;
932 
933  if(fStack) // ESD
934  {
935  TParticle * particle = fStack->Particle(ipr) ;
936 
937  imother = particle->GetFirstMother();
938  pdg = particle->GetPdgCode();
939  status = particle->GetStatusCode();
940 
941  // Compare trigger with final state particles
942  if( status != 1) continue ; // do it here to avoid crashes
943 
944  charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
945  particle->Momentum(fLVTmp);
946  }
947  else // AOD
948  {
949  AliAODMCParticle * particle = (AliAODMCParticle*) fAODMCparticles->At(ipr) ;
950 
951  imother = particle->GetMother();
952  pdg = particle->GetPdgCode();
953  status = particle->GetStatus();
954 
955  // Compare trigger with final state particles
956  if( status != 1) continue ; // do it here to avoid crashes
957 
958  charge = particle->Charge();
959  fLVTmp.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->E());
960  }
961 
962  // Do not consider the photon decays from pi0 and eta
963  //printf("Leading ipr %d - mother %d - iTrig\n",ipr, imother,indexTrig);
964  if( imother == indexTrig) continue ;
965 
966  Float_t pt = fLVTmp.Pt();
967  Float_t eta = fLVTmp.Eta();
968  Float_t phi = fLVTmp.Phi();
969  if(phi < 0 ) phi += TMath::TwoPi();
970 
971  // Select all particles in at least the TPC acceptance
972  Bool_t inTPC = GetFiducialCut()->IsInFiducialCut(eta,phi,kCTS) ;
973  if(!inTPC) continue;
974 
975  //Isolation
976  Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
977 
978  if(charge==0)
979  {
980  if(pt < fMinNeutralPt) continue ;
981 
982  if( ptMaxNeutral < pt ) ptMaxNeutral = pt;
983 
984  if( radius < rThresIC )
985  {
986  if( pt > ptThresIC ) nICNeutral++ ;
987  sumNePt+= pt;
988  }
989 
990  Bool_t phPDG = kFALSE;
991  if(pdg==22 || pdg==111) phPDG = kTRUE;
992 
993 // 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());
994  if(phPDG)
995  {
996  if( ptMaxNeutPhot < pt) ptMaxNeutPhot = pt;
997 
998  if( radius < rThresIC )
999  {
1000  if(pt > ptThresIC) nICNeutPhot++ ;
1001  sumNePtPhot += pt;
1002  }
1003  }
1004 
1005  // Calorimeter acceptance
1006  Bool_t inCalo = GetFiducialCut()->IsInFiducialCut(eta,phi,GetCalorimeter()) ;
1007  if(!inCalo) continue;
1008 
1009  if( ptMaxNeutEMCAL < pt ) ptMaxNeutEMCAL = pt;
1010  if( radius < rThresIC )
1011  {
1012  if( pt > ptThresIC ) nICNeutEMCAL++ ;
1013  sumNePtEMC += pt;
1014  }
1015 
1016  if(phPDG)
1017  {
1018  if( ptMaxNeutEMCALPhot < pt ) ptMaxNeutEMCALPhot = pt;
1019  if( radius < rThresIC )
1020  {
1021  if (pt > ptThresIC) nICNeutEMCALPhot++ ;
1022  sumNePtEMCPhot += pt;
1023  }
1024  }
1025  }
1026  else
1027  {
1028  if( pt < fMinChargedPt) continue ;
1029 
1030  if( ptMaxCharged < pt ) ptMaxCharged = pt;
1031 
1032  if( radius < rThresIC )
1033  {
1034 // 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",
1035 // ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
1036  if( pt > ptThresIC ) nICTrack++ ;
1037  sumChPt += pt;
1038  }
1039  }
1040  } // particle loop
1041 
1042  // Leding decision
1043  if(ptTrig > ptMaxCharged)
1044  {
1045 // printf("pt charged %2.2f, pt neutral %2.2f, pt neutral emcal %2.2f, pt photon %2.2f, pt photon emcal %2.2f\n",
1046 // ptMaxCharged, ptMaxNeutral, ptMaxNeutEMCAL,ptMaxNeutPhot, ptMaxNeutEMCALPhot);
1047  if(ptTrig > ptMaxNeutral ) leading[0] = kTRUE ;
1048  if(ptTrig > ptMaxNeutEMCAL ) leading[1] = kTRUE ;
1049  if(ptTrig > ptMaxNeutPhot ) leading[2] = kTRUE ;
1050  if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
1051  }
1052 
1053 // printf("N in cone over threshold: tracks %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n",
1054 // nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
1055 
1056  //------------------
1057  // Isolation decision
1058  if( isoMethod == AliIsolationCut::kPtThresIC )
1059  {
1060  if( nICTrack == 0 )
1061  {
1062  if(nICNeutral == 0 ) isolated[0] = kTRUE ;
1063  if(nICNeutEMCAL == 0 ) isolated[1] = kTRUE ;
1064  if(nICNeutPhot == 0 ) isolated[2] = kTRUE ;
1065  if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
1066  }
1067  }
1068  else if( isoMethod == AliIsolationCut::kSumPtIC )
1069  {
1070  if(sumChPt + sumNePt < sumThresIC ) isolated[0] = kTRUE ;
1071  if(sumChPt + sumNePtEMC < sumThresIC ) isolated[1] = kTRUE ;
1072  if(sumChPt + sumNePtPhot < sumThresIC ) isolated[2] = kTRUE ;
1073  if(sumChPt + sumNePtEMCPhot < sumThresIC ) isolated[3] = kTRUE ;
1074  }
1075 
1076  //----------------------------------------------------
1077  // Fill histograms if conditions apply for all 4 cases
1078  for( Int_t i = 0; i < fgkNIso; i++ )
1079  {
1080  if(leading[i])
1081  {
1082  fhPtLeading[partType][i]->Fill(ptTrig, GetEventWeight());
1083 
1084  if (i == 0) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePt , GetEventWeight());
1085  else if(i == 1) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePtEMC , GetEventWeight());
1086  else if(i == 2) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePtPhot , GetEventWeight());
1087  else if(i == 3) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePtEMCPhot, GetEventWeight());
1088 
1089  if(isolated[i]) fhPtLeadingIsolated[partType][i]->Fill(ptTrig, GetEventWeight());
1090  }
1091  } // conditions loop
1092 
1093  AliDebug(1,"End");
1094 }
1095 
1096 //_____________________________________________________
1098 //_____________________________________________________
1100 {
1101  AliDebug(1,"Start");
1102 
1103  fParton6.SetPxPyPzE(0,0,0,0);
1104  fParton7.SetPxPyPzE(0,0,0,0);
1105  fParton6PDG = 0;
1106  fParton7PDG = 0;
1107 
1108  //
1109  // Get the ESD MC particles container
1110  fStack = 0;
1111  if( GetReader()->ReadStack() )
1112  {
1113  fStack = GetMCStack();
1114  if( !fStack )
1115  {
1116  AliFatal("Stack not available, is the MC handler called? STOP");
1117  return;
1118  }
1119 
1120  fNPrimaries = fStack->GetNprimary(); // GetNtrack();
1121 
1122  if(fNPrimaries > 6)
1123  {
1124  (fStack->Particle(6))->Momentum(fParton6) ;
1125  fParton6PDG = (fStack->Particle(6))->GetPdgCode();
1126  }
1127 
1128  if(fNPrimaries > 7)
1129  {
1130  (fStack->Particle(7))->Momentum(fParton7) ;
1131  fParton7PDG = (fStack->Particle(7))->GetPdgCode();
1132  }
1133  }
1134 
1135  // Get the AOD MC particles container
1136  fAODMCparticles = 0;
1137  if( GetReader()->ReadAODMCParticles() )
1138  {
1140  if( !fAODMCparticles )
1141  {
1142  AliFatal("Standard MCParticles not available!");
1143  return;
1144  }
1145 
1146  fNPrimaries = fAODMCparticles->GetEntriesFast();
1147  AliAODMCParticle * primAOD = 0;
1148  if(fNPrimaries > 6)
1149  {
1150  primAOD = (AliAODMCParticle *) fAODMCparticles->At(6);
1151  fParton6.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1152 
1153  fParton6PDG = primAOD->GetPdgCode();
1154  }
1155 
1156  if(fNPrimaries > 7)
1157  {
1158  primAOD = (AliAODMCParticle *) fAODMCparticles->At(7);
1159  fParton7.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1160 
1161  fParton7PDG = primAOD->GetPdgCode();
1162  }
1163  }
1164 
1166 
1167  // Main particle loop
1168  Int_t pdgTrig = 0;
1169  Int_t statusTrig = 0;
1170  Int_t imother = 0;
1171  Float_t ptTrig = 0;
1172  Int_t momStatus = 0;
1173  Int_t momPdg = 0;
1174  Int_t momNDaugh = 0;
1175  Int_t momImom = 0;
1176  Int_t pdg0 = 0;
1177  Int_t pdg1 = 0;
1178  Int_t id0 = 0;
1179  Int_t id1 = 0;
1180  Int_t nDaughters = 0;
1181 
1182  for(Int_t ipr = 0; ipr < fNPrimaries; ipr ++ )
1183  {
1184  if(fStack) // ESD
1185  {
1186  TParticle * particle = fStack->Particle(ipr) ;
1187 
1188  pdgTrig = particle->GetPdgCode();
1189  statusTrig = particle->GetStatusCode();
1190  imother = particle->GetFirstMother();
1191  ptTrig = particle->Pt();
1192  nDaughters = particle->GetNDaughters();
1193  id0 = particle->GetDaughter(0);
1194  id1 = particle->GetDaughter(1);
1195  // Recover the kinematics:
1196  particle->Momentum(fTrigger);
1197  }
1198  else // AOD
1199  {
1200  AliAODMCParticle* particle = (AliAODMCParticle*) fAODMCparticles->At(ipr) ;
1201 
1202  pdgTrig = particle->GetPdgCode();
1203  statusTrig = particle->GetStatus();
1204  imother = particle->GetMother();
1205  nDaughters = particle->GetNDaughters();
1206  id0 = particle->GetDaughter(0);
1207  id1 = particle->GetDaughter(1);
1208  // Recover the kinematics:
1209  fTrigger.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->E());
1210  }
1211 
1212  //
1213  // Select final state photons, pi0s or eta's (mesons not final state)
1214  //
1215  // Final state particles have status code 1
1216  // pi0 and etas can have status code 0, (not decayed by pythia but geant?)
1217  // In the current code photons are always final state with status code 1,
1218  // avoid complications with conversions although pi0 with status code 0
1219  // generate photons with status code 0
1220 
1221  if( pdgTrig == 22 && statusTrig != 1 ) continue ;
1222 
1223  if( pdgTrig != 111 && pdgTrig != 22 && pdgTrig !=221 ) continue ;
1224 
1225  //
1226  // Pt cut
1227  //
1228  Float_t ptTrig = fTrigger.Pt();
1229 
1230  if( ptTrig < GetMinPt() ) continue ;
1231 
1232  //
1233  // Select and tag the particles being, direct photon (prompt, fragmentation or isr)
1234  // decay photon from pi0, eta or other, and pi0 or eta
1235  //
1236  Int_t partType = -1;
1237 
1238  if (pdgTrig==22 )
1239  {
1240  if(imother > 0 )
1241  {
1242  if(fStack) // ESD
1243  {
1244  momStatus = (fStack->Particle(imother))->GetStatusCode();
1245  momPdg = (fStack->Particle(imother))->GetPdgCode();
1246  momNDaugh = (fStack->Particle(imother))->GetNDaughters();
1247  momImom = (fStack->Particle(imother))->GetFirstMother();
1248  }
1249  else // AOD
1250  {
1251  momStatus = ((AliAODMCParticle*) fAODMCparticles->At(imother))->GetStatus();
1252  momPdg = ((AliAODMCParticle*) fAODMCparticles->At(imother))->GetPdgCode();
1253  momNDaugh = ((AliAODMCParticle*) fAODMCparticles->At(imother))->GetNDaughters();
1254  momImom = ((AliAODMCParticle*) fAODMCparticles->At(imother))->GetMother();
1255  }
1256 
1257  if (imother < 8 && statusTrig == 1)
1258  {
1259  partType = kmcPrimPhoton ;
1260  }
1261  else if(momPdg == 111 )
1262  {
1263  partType = kmcPrimPi0Decay;
1264  }
1265  else if(momPdg == 221 )
1266  {
1267  partType = kmcPrimEtaDecay;
1268  }
1269  else if(TMath::Abs(momStatus) > 0 )
1270  {
1271  partType = kmcPrimOtherDecay ;
1272  }
1273  }
1274  }
1275  else if( (pdgTrig==111 || pdgTrig==221) && nDaughters == 2 )
1276  {
1277  if(fStack) // ESD
1278  {
1279  pdg0 = fStack->Particle(id0)->GetPdgCode();
1280  pdg1 = fStack->Particle(id1)->GetPdgCode();
1281  }
1282  else // AOD
1283  {
1284  pdg0 = ((AliAODMCParticle*) fAODMCparticles->At(id0))->GetPdgCode();
1285  pdg1 = ((AliAODMCParticle*) fAODMCparticles->At(id1))->GetPdgCode();
1286  }
1287 
1288  if( pdg0 == 22 && pdg1== 22 )
1289  {
1290  if ( pdgTrig==111 ) partType = kmcPrimPi0;
1291  else if( pdgTrig==221 ) partType = kmcPrimEta;
1292  }
1293  }
1294  else if( (pdgTrig==111 || pdgTrig==221) )
1295  {
1296  // Fill histogram to see how many pi0/eta decay other than 2 photons in trigger detector acceptance
1298  if(! in ) continue ;
1299 
1300  if(pdgTrig==111) fhPtPi0Not2Gamma->Fill(ptTrig, GetEventWeight());
1301  if(pdgTrig==221) fhPtEtaNot2Gamma->Fill(ptTrig, GetEventWeight());
1302  }
1303 
1304  if(partType < 0 ) continue ;
1305 
1306  //
1307  // Fill particle acceptance histograms
1308  //
1309  Float_t eta = fTrigger.Eta();
1310  Float_t phi = fTrigger.Phi();
1311  if(phi < 0) phi+=TMath::TwoPi();
1312 
1313  fhPhi [partType]->Fill(phi, GetEventWeight());
1314  fhEta [partType]->Fill(eta, GetEventWeight());
1315  fhEtaPhi[partType]->Fill(eta, phi, GetEventWeight());
1316 
1317  if(partType < 4 && partType!=0)
1318  {
1319  fhPhiStatus[partType]->Fill(phi, momStatus, GetEventWeight());
1320  fhEtaStatus[partType]->Fill(eta, momStatus, GetEventWeight());
1321  }
1322  else
1323  {
1324  fhPhiStatus[partType]->Fill(phi, statusTrig, GetEventWeight());
1325  fhEtaStatus[partType]->Fill(eta, statusTrig, GetEventWeight());
1326  }
1327 
1328  //
1329  // Select particles in trigger detector acceptance
1330  //
1332  if(! in ) continue ;
1333 
1334 
1335  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",
1336  ipr, pdgTrig, partType, statusTrig, imother, ptTrig, fTrigger.Eta(), fTrigger.Phi()*TMath::RadToDeg()));
1337 
1338  //
1339  // Fill particle pT histograms, check also status
1340  //
1341 
1342  fhPt[partType]->Fill(ptTrig, GetEventWeight());
1343 
1344  if (partType==kmcPrimPi0)
1345  {
1346  fhPtPi0Status->Fill(ptTrig, statusTrig, GetEventWeight());
1347  }
1348  else if(partType==kmcPrimEta)
1349  {
1350  fhPtEtaStatus->Fill(ptTrig, statusTrig, GetEventWeight());
1351  }
1352  else if(partType == kmcPrimPi0Decay )
1353  {
1354  fhPtPi0DecayStatus ->Fill(ptTrig, momStatus, GetEventWeight());
1355 
1356  if(momNDaugh!=2) fhPtGammaFromPi0Not2Gamma->Fill(ptTrig, GetEventWeight());
1357  }
1358  else if(partType == kmcPrimEtaDecay )
1359  {
1360  fhPtEtaDecayStatus ->Fill(ptTrig, momStatus, GetEventWeight());
1361 
1362  if(momNDaugh!=2) fhPtGammaFromEtaNot2Gamma->Fill(ptTrig, GetEventWeight());
1363  }
1364  else if(partType == kmcPrimOtherDecay)
1365  {
1366  fhPtOtherDecayStatus->Fill(ptTrig, momStatus, GetEventWeight());
1367 
1368  //Fill histogram with origin of this kind of photon
1369  if (momPdg > 2100 && momPdg < 2210)
1370  fhPtOtherDecayMesonId->Fill(ptTrig, 0.5, GetEventWeight());// resonances
1371  else if(momPdg == 221) fhPtOtherDecayMesonId->Fill(ptTrig, 5.5, GetEventWeight());//eta
1372  else if(momPdg == 331) fhPtOtherDecayMesonId->Fill(ptTrig, 6.5, GetEventWeight());//eta prime
1373  else if(momPdg == 213) fhPtOtherDecayMesonId->Fill(ptTrig, 1.5, GetEventWeight());//rho
1374  else if(momPdg == 223) fhPtOtherDecayMesonId->Fill(ptTrig, 2.5, GetEventWeight());//omega
1375  else if(momPdg >= 310 && momPdg <= 323)
1376  fhPtOtherDecayMesonId->Fill(ptTrig, 3.5, GetEventWeight());//k0S, k+-,k*
1377  else if(momPdg == 130) fhPtOtherDecayMesonId->Fill(ptTrig, 3.5, GetEventWeight());//k0L
1378  else fhPtOtherDecayMesonId->Fill(ptTrig, 4.5, GetEventWeight());//other?
1379  }
1380 
1381  //
1382  // Origin of particle, which mother of the pi0 or eta or
1383  // the eta or pi0 or other meson at the origin of the decay photon
1384  //
1385  Int_t momOrgPdg = -1;
1386  Int_t momOrgStatus = -1;
1387  if(fStack)
1388  {
1389  if(momImom >=0 || imother >=0)
1390  {
1391  TParticle* mother = 0;
1392  if(partType < 4 && partType!=0 )
1393  mother = fStack->Particle(momImom);
1394  else
1395  mother = fStack->Particle(imother);
1396 
1397  momOrgPdg = TMath::Abs(mother->GetPdgCode());
1398  momOrgStatus = mother->GetStatusCode();
1399  }
1400  }
1401  else
1402  {
1403  if(momImom >=0 || imother >=0)
1404  {
1405  AliAODMCParticle* mother = 0;
1406  if(partType < 4 && partType!=0 )
1407  mother = (AliAODMCParticle*) fAODMCparticles->At(momImom);
1408  else
1409  mother = (AliAODMCParticle*) fAODMCparticles->At(imother);
1410 
1411  momOrgPdg = TMath::Abs(mother->GetPdgCode());
1412  momOrgStatus = mother->GetStatus();
1413  }
1414  }
1415 
1416  if (momOrgStatus == 21) fhPtOrigin[partType]->Fill(ptTrig, 0.5, GetEventWeight());//parton
1417  else if(momOrgPdg < 22 ) fhPtOrigin[partType]->Fill(ptTrig, 1.5, GetEventWeight());//quark
1418  else if(momOrgPdg > 2100 && momOrgPdg < 2210)
1419  fhPtOrigin[partType]->Fill(ptTrig, 2.5, GetEventWeight());// resonances
1420  else if(momOrgPdg == 221) fhPtOrigin[partType]->Fill(ptTrig, 8.5, GetEventWeight());//eta
1421  else if(momOrgPdg == 331) fhPtOrigin[partType]->Fill(ptTrig, 9.5, GetEventWeight());//eta prime
1422  else if(momOrgPdg == 213) fhPtOrigin[partType]->Fill(ptTrig, 4.5, GetEventWeight());//rho
1423  else if(momOrgPdg == 223) fhPtOrigin[partType]->Fill(ptTrig, 5.5, GetEventWeight());//omega
1424  else if(momOrgPdg >= 310 && momOrgPdg <= 323)
1425  fhPtOrigin[partType]->Fill(ptTrig, 6.5, GetEventWeight());//k0S, k+-,k*
1426  else if(momOrgPdg == 130) fhPtOrigin[partType]->Fill(ptTrig, 6.5, GetEventWeight());//k0L
1427  else if(momOrgStatus == 11 || momOrgStatus == 12 )
1428  fhPtOrigin[partType]->Fill(ptTrig, 3.5, GetEventWeight());//resonances
1429  else fhPtOrigin[partType]->Fill(ptTrig, 7.5, GetEventWeight());//other?
1430 
1431  if(statusTrig == 0)
1432  {
1433  // Histogram will not be filled for photons, leave it like this for now
1434  // in case we leave not final photons in the future
1435  if (momOrgStatus == 21) fhPtOriginNotFinal[partType]->Fill(ptTrig, 0.5, GetEventWeight());//parton
1436  else if(momOrgPdg < 22 ) fhPtOriginNotFinal[partType]->Fill(ptTrig, 1.5, GetEventWeight());//quark
1437  else if(momOrgPdg > 2100 && momOrgPdg < 2210)
1438  fhPtOriginNotFinal[partType]->Fill(ptTrig, 2.5, GetEventWeight());// resonances
1439  else if(momOrgPdg == 221) fhPtOriginNotFinal[partType]->Fill(ptTrig, 8.5, GetEventWeight());//eta
1440  else if(momOrgPdg == 331) fhPtOriginNotFinal[partType]->Fill(ptTrig, 9.5, GetEventWeight());//eta prime
1441  else if(momOrgPdg == 213) fhPtOriginNotFinal[partType]->Fill(ptTrig, 4.5, GetEventWeight());//rho
1442  else if(momOrgPdg == 223) fhPtOriginNotFinal[partType]->Fill(ptTrig, 5.5, GetEventWeight());//omega
1443  else if(momOrgPdg >= 310 && momOrgPdg <= 323)
1444  fhPtOriginNotFinal[partType]->Fill(ptTrig, 6.5, GetEventWeight());//k0S, k+-,k*
1445  else if(momOrgPdg == 130) fhPtOriginNotFinal[partType]->Fill(ptTrig, 6.5, GetEventWeight());//k0L
1446  else if(momOrgStatus == 11 || momOrgStatus == 12 )
1447  fhPtOriginNotFinal[partType]->Fill(ptTrig, 3.5, GetEventWeight());//resonances
1448  else fhPtOriginNotFinal[partType]->Fill(ptTrig, 7.5, GetEventWeight());//other?
1449  }
1450 
1451  //
1452  // Check if it is leading or isolated
1453  //
1454  Bool_t leading [fgkNIso] ;
1455  Bool_t isolated[fgkNIso] ;
1456 
1457  IsLeadingAndIsolated(ipr, partType, leading, isolated);
1458 
1459  //
1460  // Correlate trigger particle with partons or jets
1461  //
1462  Int_t iparton = -1;
1463  Int_t ok = CorrelateWithPartonOrJet(ipr, partType, leading, isolated, iparton);
1464  if(!ok) continue;
1465 
1466  //
1467  // Correlate trigger particle with hadrons
1468  //
1469  GetXE(ipr,partType,leading,isolated,iparton) ;
1470 
1471  }
1472 
1473  AliDebug(1,"End fill histograms");
1474 }
1475 
1476 //_________________________________________________________
1478 //_________________________________________________________
1480 {
1481  fTriggerDetectorString = det;
1482 
1483  if (det=="EMCAL") fTriggerDetector = kEMCAL;
1484  else if(det=="PHOS" ) fTriggerDetector = kPHOS;
1485  else if(det=="CTS") fTriggerDetector = kCTS;
1486  else if(det=="DCAL") fTriggerDetector = kDCAL;
1487  else if(det.Contains("DCAL") && det.Contains("PHOS")) fTriggerDetector = kDCALPHOS;
1488  else AliFatal(Form("Detector < %s > not known!", det.Data()));
1489 }
1490 
1491 //_____________________________________________________
1493 //_____________________________________________________
1495 {
1496  fTriggerDetector = det;
1497 
1498  if (det==kEMCAL) fTriggerDetectorString = "EMCAL";
1499  else if(det==kPHOS ) fTriggerDetectorString = "PHOS";
1500  else if(det==kCTS) fTriggerDetectorString = "CTS";
1501  else if(det==kDCAL) fTriggerDetectorString = "DCAL";
1502  else if(det==kDCALPHOS) fTriggerDetectorString = "DCAL_PHOS";
1503  else AliFatal(Form("Detector < %d > not known!", det));
1504 }
1505 
Int_t charge
Float_t GetHistoPtMax() const
Int_t pdg
TH1F * fhPtHard
! pt of parton
TH2F * fhPtEtaDecayStatus
! Input eta decay meson status
Float_t GetHistoPtMin() const
Int_t GetICMethod() const
TLorentzVector fParton6
! Final state Parton
virtual void AddToHistogramsName(TString add)
TH1F * fhPtJet
! pt of jet
const char * title
Definition: MakeQAPdf.C:26
TH2F * fhXEIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! xE away side
Bool_t ReadAODMCParticles() const
TH2F * fhPtPi0DecayStatus
! Input pi0 decay meson status
Float_t fPtHard
! Generated pT hard
TH2F * fhPtPi0Status
! Input pi0 status
virtual AliIsolationCut * GetIsolationCut()
TH2F * fhPtOriginNotFinal[fgkNmcPrimTypes]
! Input particle pt vs particle type originating it (if meson decay, its mother) if trigger is not fi...
TH2F * fhEtaPhi[fgkNmcPrimTypes]
! Input particle eta vs phi
TH2F * fhPtPartonPtHard
! pt of parton divided to pt hard, trigger is photon
TH2F * fhZJet[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! zHard
TH1F * fhPhi[fgkNmcPrimTypes]
! Input particle phi
TH1F * fhPtPi0Not2Gamma
! Input pi0 not 2 gamma decay
TH1F * fhPtGammaFromPi0Not2Gamma
! Input gamma from pi0 not 2 gamma decay
Get trigger particles/partons/jets and correlations at generator level.
virtual AliGenEventHeader * GetGenEventHeader() const
Int_t fParton7PDG
! Final state Parton PDG
TH2F * fhZPartonIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! isolated, zHard
TH2F * fhZHardIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! isolated, zHard
TH2F * fhPtPartonTypeAway[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! pt versus away side parton type
TH1F * fhPt[fgkNmcPrimTypes]
! Input particle pt
static const Int_t fgkNmcPrimTypes
Number of MC indeces for histograms arrays.
TLorentzVector fJet7
! Pythia jet close to parton in position 7
TLorentzVector fParton7
! Final state Parton
Float_t GetPtThreshold() const
TLorentzVector fJet6
! Pythia jet close to parton in position 6
Base class for CaloTrackCorr analysis algorithms.
TString fTriggerDetectorString
Detector : EMCAL, PHOS, CTS.
virtual AliFiducialCut * GetFiducialCut()
virtual TClonesArray * GetAODMCParticles() const
virtual AliHistogramRanges * GetHistogramRanges()
TH2F * fhPtPartonTypeNearIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! pt versus originating parton type
Float_t fMinChargedPt
! Minimum energy for charged particles in correlation
const Double_t ptmax
void GetPartonsAndJets()
Fill data members with partons,jets and generated pt hard.
TH2F * fhXEUEIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! xE away side
TH2F * fhXE[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! xE away side
TList * GetCreateOutputObjects()
Create histograms to be saved in output file.
AliStack * fStack
! access ESD stack
TH1F * fhPtLeadingIsolated[fgkNmcPrimTypes][fgkNIso]
! isolated
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
TH2F * fhZHard[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! zHard
TH2F * fhPtLeadingSumPt[fgkNmcPrimTypes][fgkNIso]
! pT vs sum in cone
const Double_t ptmin
TH2F * fhPhiStatus[fgkNmcPrimTypes]
! Input particle phi vs status
Int_t GetHistoNPtSumBins() const
static const Int_t fgkNIso
Number of isolation cases.
virtual Double_t GetEventWeight() const
AliAnaGeneratorKine()
Default Constructor. Initialize parameters with default values.
Int_t fTriggerDetector
Detector : EMCAL, PHOS, CTS.
AliFiducialCut * GetFiducialCutForTrigger()
TH2F * fhZParton[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! zHard
Int_t fNPrimaries
! N primaries
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Float_t GetHistoPtSumMin() const
TH2F * fhXEUE[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! xE away side
void MakeAnalysisFillHistograms()
Particle-Parton/Jet/Hadron Correlation Analysis, main method.
TH2F * fhEtaStatus[fgkNmcPrimTypes]
! Input particle eta vs status
TLorentzVector fTrigger
! Trigger momentum, avoid generating TLorentzVectors per event
Int_t GetHistoPtBins() const
TH2F * fhPtOrigin[fgkNmcPrimTypes]
! Input particle pt vs particle type originating it (if meson decay, its mother)
TH1F * fhPtEtaNot2Gamma
! Input eta not 2 gamma decay
TH2F * fhPtJetPtHard
! pt of jet divided to pt hard, trigger is photon
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TH2F * fhPtPartonTypeAwayIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! isolated, particle pt versus away side parton type
TH2F * fhPtAcceptedGammaJet[fgkNLead][fgkNIso]
! gamma-jet pair in acceptance (jet in good eta window)
TClonesArray * fAODMCparticles
! access AOD stack
void SetTriggerDetector(TString &det)
Set the calorimeter for the analysis.
TH2F * fhPtJetPtParton
! pt of parton divided to pt parton, trigger is photon
Int_t fParton6PDG
! Final state Parton PDG
Bool_t CorrelateWithPartonOrJet(Int_t indexTrig, Int_t pdgTrig, Bool_t leading[fgkNIso], Bool_t isolated[fgkNIso], Int_t &iparton)
Correlate trigger with partons or jets, get z.
virtual AliCaloTrackReader * GetReader() const
TH1F * fhPtParton
! pt of parton
TH1F * fhEta[fgkNmcPrimTypes]
! Input particle eta
Float_t Radius(Float_t etaCandidate, Float_t phiCandidate, Float_t eta, Float_t phi) const
TH2F * fhPtOtherDecayMesonId
! Decay photons, not originating from pi0 or eta, ID of the particle
DCal, not used so far, just in case.
void IsLeadingAndIsolated(Int_t indexTrig, Int_t pdgTrig, Bool_t leading[fgkNIso], Bool_t isolated[fgkNIso])
TH1F * fhPtGammaFromEtaNot2Gamma
! Input gamma from eta not 2 gamma decay
void InitParameters()
Initialize the parameters of the analysis.
Float_t GetHistoPtSumMax() const
void GetXE(Int_t indexTrig, Int_t pdgTrig, Bool_t leading[fgkNIso], Bool_t isolated[fgkNIso], Int_t iparton)
Calculate the real XE and the UE XE.
Int_t nptbins
TH1F * fhPtLeading[fgkNmcPrimTypes][fgkNIso]
! pT
static const Int_t fgkNLead
Number of leadingness cases.
TH2F * fhPtPartonTypeNear[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! particle pt versus originating parton type
Float_t fMinNeutralPt
! Minimum energy for neutral particles in correlation
Float_t GetConeSize() const
TH2F * fhPtEtaStatus
! Input eta status
TLorentzVector fLVTmp
! Momentum, avoid generating TLorentzVectors per event
TH2F * fhPtOtherDecayStatus
! Input other decay particle status
TH2F * fhZJetIsolated[fgkNmcPrimTypes][fgkNLead][fgkNIso]
! isolated, zHard