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