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