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