o add histograms for TPC + TOF
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleHadronCorrelation.cxx
CommitLineData
e09cf5ef 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/* $Id: $ */
16
17//_________________________________________________________________________
18// Class for the analysis of particle - hadron correlations
19// Particle (for example direct gamma) must be found in a previous analysis
20//-- Author: Gustavo Conesa (LNF-INFN)
21
22// Modified by Yaxian Mao:
23// 1. add the UE subtraction for corrlation study
24// 2. change the correlation variable
25// 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)
26// 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)
27// 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06)
28//////////////////////////////////////////////////////////////////////////////
29
30
31// --- ROOT system ---
32//#include "TClonesArray.h"
33#include "TClass.h"
34#include "TMath.h"
35#include "TH3D.h"
36
37//---- ANALYSIS system ----
38#include "AliNeutralMesonSelection.h"
39#include "AliAnaParticleHadronCorrelation.h"
40#include "AliCaloTrackReader.h"
41#include "AliCaloPID.h"
42#include "AliAODPWG4ParticleCorrelation.h"
43#include "AliFiducialCut.h"
44#include "AliVTrack.h"
45#include "AliVCluster.h"
46#include "AliMCAnalysisUtils.h"
47#include "TParticle.h"
48#include "AliStack.h"
49#include "AliAODMCParticle.h"
50#include "AliMixedEvent.h"
51
52ClassImp(AliAnaParticleHadronCorrelation)
53
54
55//____________________________________________________________________________
56 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
57 AliAnaPartCorrBaseClass(),
58 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
59 fMakeSeveralUE(0),
60 fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.),
61 fPi0AODBranchName(""),fNeutralCorr(0), fPi0Trigger(0),
62 // fMultiBin(0),fNZvertBin(0),fNrpBin(0),fZvtxCut(0.),
63 // fUseSelectEvent(kFALSE),
64 // fhNclustersNtracks(0), //fhVertex(0),
65 fhPtLeading(0),fhPhiLeading(0),fhEtaLeading(0),
66 fhDeltaPhiDeltaEtaCharged(0),
67 fhPhiCharged(0), fhEtaCharged(0),
68 fhDeltaPhiCharged(0),
69 fhDeltaEtaCharged(0),
70 fhDeltaPhiChargedPt(0),
71 fhDeltaPhiUeChargedPt(0),
72 fhPtImbalanceCharged(0),
73 fhPtImbalanceUeCharged(0),
74 fhPtImbalancePosCharged(0),fhPtImbalanceNegCharged(0),
75 fhPtHbpCharged(0), fhPtHbpUeCharged(0),
76 fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),
77 fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),
78 fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0),
79 fhPtTrigPout(0), fhPtAssocDeltaPhi(0),
80 fhPtTrigCharged(0),
81 fhTrigDeltaPhiCharged(0x0), fhTrigDeltaEtaCharged(0x0),fhTrigCorr(0x0),fhTrigUeCorr(0x0),
82 fhDeltaPhiDeltaEtaNeutral(0),
83 fhPhiNeutral(0), fhEtaNeutral(0),
84 fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),
85 fhDeltaPhiNeutralPt(0),fhDeltaPhiUeNeutralPt(0),
86 fhPtImbalanceNeutral(0),fhPtImbalanceUeNeutral(0),
87 fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),
88 fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),
89 fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),
90 fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0),
91 fhPtPi0DecayRatio(0),
92 fhDeltaPhiDecayCharged(0),
93 fhPtImbalanceDecayCharged(0),
94 fhDeltaPhiDecayNeutral(0),
95 fhPtImbalanceDecayNeutral(0)
96{
97 //Default Ctor
98
99 //Initialize parameters
100 InitParameters();
101}
102
103//________________________________________________________________________
104TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
105{
106
107 // Create histograms to be saved in output file and
108 // store them in fOutputContainer
109 TList * outputContainer = new TList() ;
110 outputContainer->SetName("CorrelationHistos") ;
111
112 Int_t nptbins = GetHistoPtBins();
113 Int_t nphibins = GetHistoPhiBins();
114 Int_t netabins = GetHistoEtaBins();
115 Float_t ptmax = GetHistoPtMax();
116 Float_t phimax = GetHistoPhiMax();
117 Float_t etamax = GetHistoEtaMax();
118 Float_t ptmin = GetHistoPtMin();
119 Float_t phimin = GetHistoPhiMin();
120 Float_t etamin = GetHistoEtaMin();
121
122
123// fhNclustersNtracks = new TH2F ("Multiplicity","Neutral cluster and charged track multiplicity",1000, 0., 1000.,1000, 0., 1000.);
124// fhNclustersNtracks->SetYTitle("# of tracks");
125// fhNclustersNtracks->SetXTitle("# of clusters");
126
127 fhPtLeading = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
128 fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
129
130 fhPhiLeading = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
131 fhPhiLeading->SetYTitle("#phi (rad)");
132
133 fhEtaLeading = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
134 fhEtaLeading->SetYTitle("#eta ");
135 // outputContainer->Add(fhNclustersNtracks);
136 outputContainer->Add(fhPtLeading);
137 outputContainer->Add(fhPhiLeading);
138 outputContainer->Add(fhEtaLeading);
139
140 //Correlation with charged hadrons
141 if(GetReader()->IsCTSSwitchedOn()) {
142 fhDeltaPhiDeltaEtaCharged = new TH2F
143 ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
144 140,-2.,5.,200,-2,2);
145 fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
146 fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
147
148 fhPhiCharged = new TH2F
149 ("PhiCharged","#phi_{h^{#pm}} vs p_{T #pm}",
150 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
151 fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
152 fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
153
154 fhEtaCharged = new TH2F
155 ("EtaCharged","#eta_{h^{#pm}} vs p_{T #pm}",
156 nptbins,ptmin,ptmax,netabins,etamin,etamax);
157 fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
158 fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
159
160 fhDeltaPhiCharged = new TH2F
161 ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
162 nptbins,ptmin,ptmax,140,-2.,5.);
163 fhDeltaPhiCharged->SetYTitle("#Delta #phi");
164 fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
165
166 fhDeltaPhiChargedPt = new TH2F
167 ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
168 nptbins,ptmin,ptmax,140,-2.,5.);
169 fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
170 fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
171
172 fhDeltaPhiUeChargedPt = new TH2F
173 ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
174 nptbins,ptmin,ptmax,140,-2.,5.);
175 fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
176 fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
177
178 fhDeltaEtaCharged = new TH2F
179 ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
180 nptbins,ptmin,ptmax,200,-2,2);
181 fhDeltaEtaCharged->SetYTitle("#Delta #eta");
182 fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
183
184 fhPtImbalanceCharged =
185 new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
186 nptbins,ptmin,ptmax,200,0.,2.);
187 fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
188 fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
189
190 fhPtImbalanceUeCharged =
191 new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
192 nptbins,ptmin,ptmax,200,0.,2.);
193 fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
194 fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
195
196 fhPtImbalancePosCharged =
197 new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",
198 nptbins,ptmin,ptmax,200,0.,2.);
199 fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");
200 fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");
201
202 fhPtImbalanceNegCharged =
203 new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",
204 nptbins,ptmin,ptmax,200,0.,2.);
205 fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");
206 fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");
207
208 fhPtHbpCharged =
209 new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
210 nptbins,ptmin,ptmax,200,0.,10.);
211 fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
212 fhPtHbpCharged->SetXTitle("p_{T trigger}");
213
214 fhPtHbpUeCharged =
215 new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
216 nptbins,ptmin,ptmax,200,0.,10.);
217 fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
218 fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
219
220 fhPtTrigPout =
221 new TH2F("PtTrigPout","Pout with triggers",
222 nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
223 fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");
224 fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
225
226 fhPtAssocDeltaPhi =
227 new TH2F("fhPtAssocDeltaPhi"," charged hadrons vs. delta phi",
228 nptbins,ptmin,ptmax,140,-2.,5.);
229 fhPtAssocDeltaPhi->SetXTitle("p_{T h^{#pm}} (GeV/c)");
230 fhPtAssocDeltaPhi->SetYTitle("#Delta #phi (GeV/c)");
231
232 fhPtTrigCharged =
233 new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",
234 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
235 fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");
236 fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
237
238 outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
239 outputContainer->Add(fhPhiCharged) ;
240 outputContainer->Add(fhEtaCharged) ;
241 outputContainer->Add(fhDeltaPhiCharged) ;
242 outputContainer->Add(fhDeltaEtaCharged) ;
243 outputContainer->Add(fhDeltaPhiChargedPt) ;
244 outputContainer->Add(fhDeltaPhiUeChargedPt) ;
245 outputContainer->Add(fhPtImbalanceCharged) ;
246 outputContainer->Add(fhPtImbalancePosCharged) ;
247 outputContainer->Add(fhPtImbalanceNegCharged) ;
248 outputContainer->Add(fhPtImbalanceUeCharged) ;
249 outputContainer->Add(fhPtHbpCharged) ;
250 outputContainer->Add(fhPtHbpUeCharged) ;
251 outputContainer->Add(fhPtTrigPout) ;
252 outputContainer->Add(fhPtAssocDeltaPhi) ;
253 outputContainer->Add(fhPtTrigCharged) ;
254
255 if(DoEventSelect()){
256 Int_t nMultiBins = GetMultiBin();
257 fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;
258 fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;
259 fhTrigCorr = new TH2F*[nMultiBins];
260 fhTrigUeCorr = new TH2F*[nMultiBins];
261 for(Int_t im=0; im<nMultiBins; im++){
262 fhTrigDeltaPhiCharged[im] = new TH2F
263 (Form("fhTrigDeltaPhiCharged_%d",im),Form("fhTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.);
264 fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
265 fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
266 fhTrigDeltaEtaCharged[im] = new TH2F
267 (Form("fhTrigDeltaEtaCharged_%d",im),Form("fhTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2);
268 fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
269 fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
270 fhTrigCorr[im] = new TH2F
271 (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
272 fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
273 fhTrigCorr[im]->SetXTitle("p_{T trigger}");
274 fhTrigUeCorr[im] = new TH2F
275 (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
276 fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
277 fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");
278
279 outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;
280 outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;
281 outputContainer->Add(fhTrigCorr[im]);
282 outputContainer->Add(fhTrigUeCorr[im]);
283
284 }
285 }
286
287 if(fPi0Trigger){
288 fhPtPi0DecayRatio = new TH2F
289 ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
290 nptbins,ptmin,ptmax, 100,0.,2.);
291 fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
292 fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
293
294 fhDeltaPhiDecayCharged = new TH2F
295 ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
296 nptbins,ptmin,ptmax,140,-2.,5.);
297 fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
298 fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
299
300 fhPtImbalanceDecayCharged =
301 new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
302 nptbins,ptmin,ptmax,200,0.,2.);
303 fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");
304 fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");
305
306 outputContainer->Add(fhPtPi0DecayRatio) ;
307 outputContainer->Add(fhDeltaPhiDecayCharged) ;
308 outputContainer->Add(fhPtImbalanceDecayCharged) ;
309 }
310
311
312 if(fMakeSeveralUE){
313 fhDeltaPhiUeLeftCharged = new TH2F
314 ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
315 nptbins,ptmin,ptmax,140,-2.,5.);
316 fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
317 fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
318 outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
319
320 fhDeltaPhiUeRightCharged = new TH2F
321 ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
322 nptbins,ptmin,ptmax,140,-2.,5.);
323 fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
324 fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
325 outputContainer->Add(fhDeltaPhiUeRightCharged) ;
326
327 fhPtImbalanceUeLeftCharged =
328 new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
329 nptbins,ptmin,ptmax,200,0.,2.);
330 fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
331 fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
332 outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
333
334 fhPtImbalanceUeRightCharged =
335 new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
336 nptbins,ptmin,ptmax,200,0.,2.);
337 fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
338 fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
339 outputContainer->Add(fhPtImbalanceUeRightCharged) ;
340
341 fhPtHbpUeLeftCharged =
342 new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
343 nptbins,ptmin,ptmax,200,0.,10.);
344 fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
345 fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
346 outputContainer->Add(fhPtHbpUeLeftCharged) ;
347
348 fhPtHbpUeRightCharged =
349 new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
350 nptbins,ptmin,ptmax,200,0.,10.);
351 fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
352 fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
353 outputContainer->Add(fhPtHbpUeRightCharged) ;
354
355 }
356 } //Correlation with charged hadrons
357
358 //Correlation with neutral hadrons
359 if(fNeutralCorr){
360
361 fhDeltaPhiDeltaEtaNeutral = new TH2F
362 ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
363 140,-2.,5.,200,-2,2);
364 fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
365 fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
366
367 fhPhiNeutral = new TH2F
368 ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #pi^{0}}",
369 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
370 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
371 fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
372
373 fhEtaNeutral = new TH2F
374 ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #pi^{0}}",
375 nptbins,ptmin,ptmax,netabins,etamin,etamax);
376 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
377 fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
378
379 fhDeltaPhiNeutral = new TH2F
380 ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
381 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
382 fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
383 fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
384
385 fhDeltaPhiNeutralPt = new TH2F
386 ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
387 nptbins,ptmin,ptmax,140,-2.,5.);
388 fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
389 fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
390
391 fhDeltaPhiUeNeutralPt = new TH2F
392 ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
393 nptbins,ptmin,ptmax,140,-2.,5.);
394 fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
395 fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
396
397 fhDeltaEtaNeutral = new TH2F
398 ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
399 nptbins,ptmin,ptmax,200,-2,2);
400 fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
401 fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
402
403 fhPtImbalanceNeutral =
404 new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
405 nptbins,ptmin,ptmax,200,0.,2.);
406 fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
407 fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
408
409 fhPtImbalanceUeNeutral =
410 new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
411 nptbins,ptmin,ptmax,200,0.,2.);
412 fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
413 fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
414
415 fhPtHbpNeutral =
416 new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
417 nptbins,ptmin,ptmax,200,0.,10.);
418 fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
419 fhPtHbpNeutral->SetXTitle("p_{T trigger}");
420
421 fhPtHbpUeNeutral =
422 new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
423 nptbins,ptmin,ptmax,200,0.,10.);
424 fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
425 fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
426
427
428 outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
429 outputContainer->Add(fhPhiNeutral) ;
430 outputContainer->Add(fhEtaNeutral) ;
431 outputContainer->Add(fhDeltaPhiNeutral) ;
432 outputContainer->Add(fhDeltaPhiNeutralPt) ;
433 outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
434 outputContainer->Add(fhDeltaEtaNeutral) ;
435 outputContainer->Add(fhPtImbalanceNeutral) ;
436 outputContainer->Add(fhPtImbalanceUeNeutral) ;
437 outputContainer->Add(fhPtHbpNeutral) ;
438 outputContainer->Add(fhPtHbpUeNeutral) ;
439
440 if(fPi0Trigger){
441 fhDeltaPhiDecayNeutral = new TH2F
442 ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
443 nptbins,ptmin,ptmax,140,-2.,5.);
444 fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
445 fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
446
447 fhPtImbalanceDecayNeutral =
448 new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
449 nptbins,ptmin,ptmax,200,0.,2.);
450 fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");
451 fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");
452
453 outputContainer->Add(fhDeltaPhiDecayNeutral) ;
454 outputContainer->Add(fhPtImbalanceDecayNeutral) ;
455 }
456
457
458 if(fMakeSeveralUE){
459 fhDeltaPhiUeLeftNeutral = new TH2F
460 ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
461 nptbins,ptmin,ptmax,140,-2.,5.);
462 fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
463 fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
464 outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
465
466 fhDeltaPhiUeRightNeutral = new TH2F
467 ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
468 nptbins,ptmin,ptmax,140,-2.,5.);
469 fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
470 fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
471 outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
472
473 fhPtImbalanceUeLeftNeutral =
474 new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
475 nptbins,ptmin,ptmax,140,0.,2.);
476 fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
477 fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
478 outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
479
480 fhPtImbalanceUeRightNeutral =
481 new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
482 nptbins,ptmin,ptmax,200,0.,2.);
483 fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
484 fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
485 outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
486
487 fhPtHbpUeLeftNeutral =
488 new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
489 nptbins,ptmin,ptmax,200,0.,10.);
490 fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
491 fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
492 outputContainer->Add(fhPtHbpUeLeftNeutral) ;
493
494 fhPtHbpUeRightNeutral =
495 new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
496 nptbins,ptmin,ptmax,200,0.,10.);
497 fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
498 fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
499 outputContainer->Add(fhPtHbpUeRightNeutral) ;
500
501 }
502
503 //Keep neutral meson selection histograms if requiered
504 //Setting done in AliNeutralMesonSelection
505 if(GetNeutralMesonSelection()){
506 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
507 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
508 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
509 delete nmsHistos;
510 }
511
512 }//Correlation with neutral hadrons
513
514 return outputContainer;
515
516}
517
518//____________________________________________________________________________
519void AliAnaParticleHadronCorrelation::InitParameters()
520{
521
522 //Initialize the parameters of the analysis.
523 SetInputAODName("PWG4Particle");
524 SetAODObjArrayName("Hadrons");
525 AddToHistogramsName("AnaHadronCorr_");
526
527 SetPtCutRange(0.,300);
528 fDeltaPhiMinCut = 1.5 ;
529 fDeltaPhiMaxCut = 4.5 ;
530 fSelectIsolated = kFALSE;
531 fMakeSeveralUE = kFALSE;
532 fUeDeltaPhiMinCut = 1. ;
533 fUeDeltaPhiMaxCut = 1.5 ;
534 fNeutralCorr = kFALSE ;
535 fPi0Trigger = kFALSE ;
536
537}
538
539//__________________________________________________________________
540void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
541{
542
543 //Print some relevant parameters set for the analysis
544 if(! opt)
545 return;
546
547 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
548 AliAnaPartCorrBaseClass::Print(" ");
549
550 printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ;
551 printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;
552 printf("Isolated Trigger? %d\n", fSelectIsolated) ;
553 printf("Phi trigger particle-UeHadron < %3.2f\n", fUeDeltaPhiMaxCut) ;
554 printf("Phi trigger particle-UeHadron > %3.2f\n", fUeDeltaPhiMinCut) ;
555 printf("Several UE? %d\n", fMakeSeveralUE) ;
556 printf("Name of AOD Pi0 Branch %s \n",fPi0AODBranchName.Data());
557 printf("Do Decay-hadron correlation ? %d\n", fPi0Trigger) ;
558
559
560}
561
562//__________________________________________________________________
563TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
564{
565 //Save parameters used for analysis
566 TString parList ; //this will be list of parameters used for this analysis.
567 const Int_t buffersize = 255;
568 char onePar[buffersize] ;
569
570 snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
571 parList+=onePar ;
572 snprintf(onePar,buffersize,"Phi trigger particle-Hadron < %3.2f ", fDeltaPhiMaxCut) ;
573 parList+=onePar ;
574 snprintf(onePar,buffersize,"Phi trigger particle-Hadron > %3.2f ", fDeltaPhiMinCut) ;
575 parList+=onePar ;
576 snprintf(onePar,buffersize,"Isolated Trigger? %d\n", fSelectIsolated) ;
577 parList+=onePar ;
578 snprintf(onePar,buffersize,"Phi trigger particle-UeHadron < %3.2f ", fUeDeltaPhiMaxCut) ;
579 parList+=onePar ;
580 snprintf(onePar,buffersize,"Phi trigger particle-UeHadron > %3.2f ", fUeDeltaPhiMinCut) ;
581 parList+=onePar ;
582 snprintf(onePar,buffersize,"Several UE? %d\n", fMakeSeveralUE) ;
583 parList+=onePar ;
584 snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ",fPi0AODBranchName.Data());
585 parList+=onePar ;
586 snprintf(onePar,buffersize,"Do Decay-hadron correlation ? %d", fPi0Trigger) ;
587 parList+=onePar ;
588
589 //Get parameters set in base class.
590 parList += GetBaseParametersList() ;
591
592 //Get parameters set in PID class.
593 //parList += GetCaloPID()->GetPIDParametersList() ;
594
595 //Get parameters set in FiducialCut class (not available yet)
596 //parlist += GetFidCut()->GetFidCutParametersList()
597
598 return new TObjString(parList) ;
599
600}
601
602
603//____________________________________________________________________________
604void AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
605{
606 //Particle-Hadron Correlation Analysis, fill AODs
607
608 if(!GetInputAODBranch()){
609 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
610 abort();
611 }
612
613 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
614 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
615 abort();
616 }
617
618 if(GetDebug() > 1){
619 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
620 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
621 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());
622 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
623 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetPHOSClusters()->GetEntriesFast());
624 }
625
626 //Loop on stored AOD particles, trigger
627 Double_t ptTrig = 0.;
628 Int_t trigIndex = -1;
629 Int_t naod = GetInputAODBranch()->GetEntriesFast();
630 //fhNclustersNtracks->Fill(naod, GetCTSTracks()->GetEntriesFast());
631 for(Int_t iaod = 0; iaod < naod ; iaod++){
632 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
633 //find the leading particles with highest momentum
634 if (particle->Pt()>ptTrig) {
635 ptTrig = particle->Pt() ;
636 trigIndex = iaod ;
637 }
638 }//Aod branch loop
639
640 //Do correlation with leading particle
641 if(trigIndex!=-1){
642
643 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
644 //Make correlation with charged hadrons
645 if(GetReader()->IsCTSSwitchedOn() )
646 MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);
647
648 TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
649 if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
650 MakeNeutralCorrelation(particle, pi0list,kFALSE);
651
652 }//Correlate leading
653
654 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
655
656}
657
658//____________________________________________________________________________
659void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
660{
661 //Particle-Hadron Correlation Analysis, fill histograms
662
663 if(!GetInputAODBranch()){
664 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
665 abort();
666 }
667
668 if(GetDebug() > 1){
669 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
670 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
671 }
672
673 //Get the vertex and check it is not too large in z
674 Double_t v[3] = {0,0,0}; //vertex ;
675 GetReader()->GetVertex(v);
676 if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
677
678 //Loop on stored AOD particles, find leading
679 Int_t naod = GetInputAODBranch()->GetEntriesFast();
680 // if(naod!=0)fhVertex->Fill(v[0],v[1],v[2]);
681 Double_t ptTrig = 0.;
682 Int_t trigIndex = -1 ;
683 for(Int_t iaod = 0; iaod < naod ; iaod++){ //loop on input trigger AOD file
684 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
685 //vertex cut in case of mixing
686 if (GetMixedEvent()) {
687 Int_t evt=-1;
688 Int_t id =-1;
689 if (particle->GetCaloLabel(0) >= 0 ){
690 id=particle->GetCaloLabel(0);
691 if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
692 }
693 else if(particle->GetTrackLabel(0) >= 0 ){
694 id=particle->GetTrackLabel(0);
695 if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
696 }
697 else continue;
698
699 if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
700 return ;
701 }
702
703 //check if the particle is isolated or if we want to take the isolation into account
704 if(OnlyIsolated() && !particle->IsIsolated()) continue;
705 //check if inside the vertex cut
706 //find the leading particles with highest momentum
707 if (particle->Pt()>ptTrig) {
708 ptTrig = particle->Pt() ;
709 trigIndex = iaod ;
710 }
711 }//finish searching for leading trigger particle
712 if(trigIndex!=-1){ //using trigger partilce to do correlations
713 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
714
715 if (GetMixedEvent()) {
716 Int_t evt=-1;
717 Int_t id = 0;
718 if (particle->GetCaloLabel(0) >= 0 ){
719 id=particle->GetCaloLabel(0);
720 if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
721 }
722 else if(particle->GetTrackLabel(0) >= 0 ){
723 id=particle->GetTrackLabel(0);
724 if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
725 }
726 else return;
727
728 if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
729 return ;
730 }
731
732 //Fill leading particle histogram
733 fhPtLeading->Fill(particle->Pt());
734 Float_t phi = particle->Phi();
735 if(phi<0)phi+=TMath::TwoPi();
736 fhPhiLeading->Fill(particle->Pt(), phi);
737 fhEtaLeading->Fill(particle->Pt(), particle->Eta());
738
739 //Make correlation with charged hadrons
740 if(GetReader()->IsCTSSwitchedOn() )
741 MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);
742
743 TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
744 if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
745 MakeNeutralCorrelation(particle, pi0list,kTRUE);
746
747 }//Aod branch loop
748
749 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
750
751}
752
753//____________________________________________________________________________
754void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)
755{
756 // Charged Hadron Correlation Analysis
757 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
758
759 Int_t evtIndex11 = -1 ; //cluster trigger or pi0 trigger
760 Int_t evtIndex12 = -1 ; // pi0 trigger
761 Int_t evtIndex13 = -1 ; // charged trigger
762 Int_t indexPhoton1 = -1 ;
763 Int_t indexPhoton2 = -1 ;
764 Double_t v[3] = {0,0,0}; //vertex ;
765 GetReader()->GetVertex(v);
766 if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
767
768 Int_t nTracks = GetCTSTracks()->GetEntriesFast() ;
769
770 if (GetMixedEvent()) {
771 evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
772 evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
773 evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
774 }
775
776 Double_t ptTrig = aodParticle->Pt();
777 Double_t pxTrig = aodParticle->Px();
778 Double_t pyTrig = aodParticle->Py();
779
780 Double_t phiTrig = aodParticle->Phi();
781 //Double_t etaTrig = aodParticle->Eta();
782
783 Double_t pt = -100.;
784 Double_t px = -100.;
785 Double_t py = -100.;
786 Double_t rat = -100.;
787 Double_t xE = -100.;
788 Double_t cosi = -100.;
789 Double_t phi = -100. ;
790 Double_t eta = -100. ;
791 TVector3 p3;
792
793 TObjArray * reftracks = 0x0;
794 Int_t nrefs = 0;
795
796 Double_t ptDecay1 = 0. ;
797 Double_t pxDecay1 = 0. ;
798 Double_t pyDecay1 = 0. ;
799 Double_t phiDecay1 = 0. ;
800 Double_t ptDecay2 = 0. ;
801 Double_t pxDecay2 = 0. ;
802 Double_t pyDecay2 = 0. ;
803 Double_t phiDecay2 = 0. ;
804
805 Double_t ratDecay1 = -100.;
806 Double_t ratDecay2 = -100.;
807 Float_t deltaphi = -100. ;
808 Float_t deltaphiDecay1 = -100. ;
809 Float_t deltaphiDecay2 = -100. ;
810 TObjArray * clusters = 0x0 ;
811 TLorentzVector photonMom ;
812
813 if(fPi0Trigger){
814 indexPhoton1 = aodParticle->GetCaloLabel (0);
815 indexPhoton2 = aodParticle->GetCaloLabel (1);
816 if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
817
818 if(indexPhoton1!=-1 && indexPhoton2!=-1){
819 if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
820 else clusters = GetPHOSClusters() ;
821 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
822 AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
823 photon->GetMomentum(photonMom,GetVertex(0)) ;
824 if(photon->GetID()==indexPhoton1) {
825 ptDecay1 = photonMom.Pt();
826 pxDecay1 = photonMom.Px();
827 pyDecay1 = photonMom.Py();
828 phiDecay1 = photonMom.Phi();
829 if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);
830 }
831 if(photon->GetID()==indexPhoton2) {
832 ptDecay2 = photonMom.Pt();
833 pxDecay2 = photonMom.Px();
834 pyDecay2 = photonMom.Py();
835 phiDecay2 = photonMom.Phi();
836 if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);
837 }
838 if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
839 } //cluster loop
840 } //index of decay photons found
841
842 } //make decay-hadron correlation
843
844 //Track loop, select tracks with good pt, phi and fill AODs or histograms
845 //Int_t currentIndex = -1 ;
846 for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
847 AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
848
849 //check if inside the vertex cut
850 //printf("charge = %d\n", track->Charge());
851 Int_t evtIndex2 = 0 ;
852 if (GetMixedEvent()) {
853 evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
854 if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
855 continue ;
856 //vertex cut
857 if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
858 return ;
859 // if(currentIndex == evtIndex2) // tracks from different event
860 // continue ;
861 // currentIndex = evtIndex2 ;
862 }
863 Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
864 p3.SetXYZ(mom[0],mom[1],mom[2]);
865 pt = p3.Pt();
866 px = p3.Px();
867 py = p3.Py();
868 eta = p3.Eta();
869 phi = p3.Phi() ;
870 if(phi < 0) phi+=TMath::TwoPi();
871
872 //Select only hadrons in pt range
873 if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
874 //remove trigger itself for correlation when use charged triggers
875 if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
876 track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3) ) continue ;
877 //&&pt==ptTrig && phi==phiTrig && eta==etaTrig) continue ;
878
879 if(IsFiducialCutOn()){
880 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
881 if(! in ) continue ;
882 }
883 //jumped out this event if near side associated partile pt larger than trigger
884 if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;
885 rat = pt/ptTrig ;
886 xE = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
887 if(xE <0.)xE =-xE;
888 cosi = TMath::Log(1/xE);
889 // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
890 // printf("phi = %f \n", phi);
891
892 if(fPi0Trigger){
893 if(indexPhoton1!=-1 && indexPhoton2!=-1){
894 if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
895 if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
896 deltaphiDecay1 = phiDecay1-phi;
897 deltaphiDecay2 = phiDecay2-phi;
898 if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();
899 if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();
900 if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();
901 if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();
902 }
903 } //do decay-hadron correlation
904
905 //Selection within angular range
906 deltaphi = phiTrig-phi;
907 if(deltaphi< -TMath::PiOver2()) deltaphi+=TMath::TwoPi();
908 if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
909
910 Double_t pout = pt*TMath::Sin(deltaphi) ;
911
912 if(GetDebug() > 2)
913 printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
914 pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());
915
916 if(bFillHisto){
917 // Fill Histograms
918 fhEtaCharged->Fill(pt,eta);
919 fhPhiCharged->Fill(pt,phi);
920 fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
921 fhDeltaPhiCharged->Fill(ptTrig, deltaphi);
922 fhDeltaPhiDeltaEtaCharged->Fill(deltaphi,aodParticle->Eta()-eta);
923 fhPtAssocDeltaPhi->Fill(pt, deltaphi);
924
925 if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
926 //fill different multiplicity histogram
927 if(DoEventSelect()){
928 for(Int_t im=0; im<GetMultiBin(); im++){
929 if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1)){
930 fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaphi);
931 fhTrigDeltaEtaCharged[im]->Fill(ptTrig,aodParticle->Eta()-eta);
932 }
933 }
934 }
935 //delta phi cut for correlation
936 if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
937 fhDeltaPhiChargedPt->Fill(pt,deltaphi);
938 fhPtImbalanceCharged->Fill(ptTrig,xE);
939 fhPtHbpCharged->Fill(ptTrig,cosi);
940 fhPtTrigPout->Fill(ptTrig, pout) ;
941 fhPtTrigCharged->Fill(ptTrig, pt) ;
942 if(track->Charge()>0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;
943 else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;
944 //fill different multiplicity histogram
945 if(DoEventSelect()){
946 for(Int_t im=0; im<GetMultiBin(); im++){
947 if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
948 fhTrigCorr[im]->Fill(ptTrig,xE);
949 }
950 } //multiplicity events selection
951 } //delta phi cut for correlation
952 else if ((deltaphi > fUeDeltaPhiMinCut) && ( deltaphi < fUeDeltaPhiMaxCut)) { //UE study
953 fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);
954 //fhUePoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
955 fhPtImbalanceUeCharged->Fill(ptTrig,xE);
956 fhPtHbpUeCharged->Fill(ptTrig,cosi);
957 if(DoEventSelect()){
958 for(Int_t im=0; im<GetMultiBin(); im++){
959 if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
960 fhTrigUeCorr[im]->Fill(ptTrig,xE);
961 }
962 } //multiplicity events selection
963
964 } //UE study
965
966 if(fPi0Trigger){
967 if(indexPhoton1!=-1 && indexPhoton2!=-1){
968 fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaphiDecay1);
969 fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaphiDecay2);
970 if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
971 if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
972 fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1);
973 if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
974 fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);
975 if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
976 } //index of decay photons found
977 } //make decay-hadron correlation
978
979 //several UE calculation
980 if(fMakeSeveralUE){
981 if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
982 fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);
983 fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
984 fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
985 }
986 if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
987 fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);
988 fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
989 fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
990
991 }
992 } //several UE calculation
993
994 } //Fill histogram
995 else{
996 nrefs++;
997 if(nrefs==1){
998 reftracks = new TObjArray(0);
999 reftracks->SetName(GetAODObjArrayName()+"Tracks");
1000 reftracks->SetOwner(kFALSE);
1001 }
1002 reftracks->Add(track);
1003 }//aod particle loop
1004 }// track loop
1005
1006 //Fill AOD with reference tracks, if not filling histograms
1007 if(!bFillHisto && reftracks) {
1008 aodParticle->AddObjArray(reftracks);
1009 }
1010
1011 //delete reftracks;
1012
1013}
1014
1015
1016
1017//____________________________________________________________________________
1018//void AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector)
1019//{
1020// // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed
1021// if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");
1022//
1023// if(!NewOutputAOD()){
1024// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");
1025// abort();
1026// }
1027//
1028// Double_t phiTrig = aodParticle->Phi();
1029// Int_t tag = 0;
1030// TLorentzVector gammai;
1031// TLorentzVector gammaj;
1032//
1033// //Get vertex for photon momentum calculation
1034//
1035// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1036// {
1037// for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {
1038// if (!GetMixedEvent())
1039// GetReader()->GetVertex(GetVertex(iev));
1040// else
1041// GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev));
1042// }
1043// }
1044// Double_t vertex2[] = {0.0,0.0,0.0} ; //vertex of second input aod
1045// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1046// {
1047// if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1048// }
1049//
1050// //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
1051// //Int_t iEvent= GetReader()->GetEventNumber() ;
1052// Int_t nclus = pl->GetEntriesFast();
1053// for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){
1054// AliVCluster * calo = (AliVCluster *) (pl->At(iclus)) ;
1055//
1056// Int_t evtIndex1 = 0 ;
1057// if (GetMixedEvent()) {
1058// evtIndex1=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1059// }
1060//
1061//
1062// //Input from second AOD?
1063// Int_t inputi = 0;
1064// if (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= iclus)
1065// inputi = 1 ;
1066// else if(aodParticle->GetDetector() == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= iclus)
1067// inputi = 1;
1068//
1069// //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
1070// //FIXME
1071// Int_t pdg=0;
1072// //if (inputi == 0 && !SelectCluster(calo, GetVertex(evtIndex1), gammai, pdg))
1073// continue ;
1074// //MEFIX
1075// else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg))
1076// continue ;
1077//
1078// if(GetDebug() > 2)
1079// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral cluster in %s: pt %f, phi %f, phi trigger %f. Cuts: delta phi min %2.2f, max %2.2f, pT min %2.2f \n",
1080// detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
1081//
1082// //2 gamma overlapped, found with PID
1083// if(pdg == AliCaloPID::kPi0){
1084//
1085// //Select only hadrons in pt range
1086// if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt())
1087// continue ;
1088//
1089// //Selection within angular range
1090// Float_t phi = gammai.Phi();
1091// if(phi < 0) phi+=TMath::TwoPi();
1092// //Float_t deltaphi = TMath::Abs(phiTrig-phi);
1093// //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
1094//
1095// AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);
1096// //pi0.SetLabel(calo->GetLabel());
1097// pi0.SetPdg(AliCaloPID::kPi0);
1098// pi0.SetDetector(detector);
1099//
1100// if(IsDataMC()){
1101// pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(),inputi));
1102// if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());
1103// }//Work with stack also
1104// //Set the indeces of the original caloclusters
1105// pi0.SetCaloLabel(calo->GetID(),-1);
1106// AddAODParticle(pi0);
1107//
1108// if(GetDebug() > 2)
1109// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());
1110//
1111// }// pdg = 111
1112//
1113// //Make invariant mass analysis
1114// else if(pdg == AliCaloPID::kPhoton){
1115// //Search the photon companion in case it comes from a Pi0 decay
1116// //Apply several cuts to select the good pair;
1117// for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
1118// AliVCluster * calo2 = (AliVCluster *) (pl->At(jclus)) ;
1119// Int_t evtIndex2 = 0 ;
1120// if (GetMixedEvent()) {
1121// evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1122// }
1123// if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1124// continue ;
1125//
1126// //Input from second AOD?
1127// Int_t inputj = 0;
1128// if (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= jclus)
1129// inputj = 1;
1130// else if(aodParticle->GetDetector() == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= jclus)
1131// inputj = 1;
1132//
1133// //Cluster selection, not charged with photon or pi0 id and in fiducial cut
1134// Int_t pdgj=0;
1135// //FIXME
1136// //if (inputj == 0 && !SelectCluster(calo2, GetVertex(evtIndex2), gammaj, pdgj))
1137// continue ;
1138// //MEFIX
1139//
1140// else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))
1141// continue ;
1142// //FIXME
1143// //if(!SelectCluster(calo2,GetVertex(evtIndex2), gammaj, pdgj))
1144// //MEFIX
1145// continue ;
1146//
1147// if(pdgj == AliCaloPID::kPhoton ){
1148//
1149// if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt())
1150// continue ;
1151//
1152// //Selection within angular range
1153// Float_t phi = (gammai+gammaj).Phi();
1154// if(phi < 0) phi+=TMath::TwoPi();
1155// //Float_t deltaphi = TMath::Abs(phiTrig-phi);
1156// //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
1157//
1158// //Select good pair (aperture and invariant mass)
1159// if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
1160//
1161// if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
1162// (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
1163//
1164// TLorentzVector pi0mom = gammai+gammaj;
1165// AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);
1166// //pi0.SetLabel(calo->GetLabel());
1167// pi0.SetPdg(AliCaloPID::kPi0);
1168// pi0.SetDetector(detector);
1169// if(IsDataMC()){
1170// //Check origin of the candidates
1171//
1172// Int_t label1 = calo->GetLabel();
1173// Int_t label2 = calo2->GetLabel();
1174// Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);
1175// Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);
1176//
1177// if(GetDebug() > 0)
1178// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1179// if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
1180//
1181// //Check if pi0 mother is the same
1182// if(GetReader()->ReadStack()){
1183// TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1184// label1 = mother1->GetFirstMother();
1185// //mother1 = GetMCStack()->Particle(label1);//pi0
1186//
1187// TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1188// label2 = mother2->GetFirstMother();
1189// //mother2 = GetMCStack()->Particle(label2);//pi0
1190// }
1191// else if(GetReader()->ReadAODMCParticles()){
1192// AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree
1193// label1 = mother1->GetMother();
1194// //mother1 = GetMCStack()->Particle(label1);//pi0
1195// AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree
1196// label2 = mother2->GetMother();
1197// //mother2 = GetMCStack()->Particle(label2);//pi0
1198// }
1199//
1200// //printf("mother1 %d, mother2 %d\n",label1,label2);
1201// if(label1 == label2)
1202// GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1203// }
1204// }//Work with mc information also
1205// pi0.SetTag(tag);
1206// //Set the indeces of the original caloclusters
1207// pi0.SetCaloLabel(calo->GetID(), calo2->GetID());
1208// AddAODParticle(pi0);
1209//
1210//
1211// }//Pair selected
1212// }//if pair of gammas
1213// }//2nd loop
1214// }// if pdg = 22
1215// }//1st loop
1216//
1217// if(GetDebug() > 1)
1218// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());
1219//}
1220
1221//____________________________________________________________________________
1222void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, TObjArray* pi0list, const Bool_t bFillHisto)
1223{
1224 // Neutral Pion Correlation Analysis
1225 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());
1226
1227 Int_t evtIndex11 = 0 ;
1228 Int_t evtIndex12 = 0 ;
1229 if (GetMixedEvent()) {
1230 evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
1231 evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
1232 }
1233
1234 Double_t pt = -100.;
1235 Double_t px = -100.;
1236 Double_t py = -100.;
1237 Double_t rat = -100.;
1238 Double_t phi = -100.;
1239 Double_t eta = -100.;
1240 Double_t xE = -100.;
1241 Double_t cosi = -100.;
1242
1243 Double_t ptTrig = aodParticle->Pt();
1244 Double_t phiTrig = aodParticle->Phi();
1245 Double_t etaTrig = aodParticle->Eta();
1246 Double_t pxTrig = aodParticle->Px();
1247 Double_t pyTrig = aodParticle->Py();
1248
1249
1250 Int_t indexPhoton1 = -1 ;
1251 Int_t indexPhoton2 = -1 ;
1252 Double_t ptDecay1 = 0. ;
1253 Double_t pxDecay1 = 0. ;
1254 Double_t pyDecay1 = 0. ;
1255 Double_t phiDecay1 = 0. ;
1256 Double_t ptDecay2 = 0. ;
1257 Double_t pxDecay2 = 0. ;
1258 Double_t pyDecay2 = 0. ;
1259 Double_t phiDecay2 = 0. ;
1260
1261 Double_t ratDecay1 = -100.;
1262 Double_t ratDecay2 = -100.;
1263 Float_t deltaphi = -100. ;
1264 Float_t deltaphiDecay1 = -100. ;
1265 Float_t deltaphiDecay2 = -100. ;
1266 TObjArray * clusters = 0x0 ;
1267 TLorentzVector photonMom ;
1268 if(fPi0Trigger){
1269 indexPhoton1 = aodParticle->GetCaloLabel (0);
1270 indexPhoton2 = aodParticle->GetCaloLabel (1);
1271 if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
1272
1273 if(indexPhoton1!=-1 && indexPhoton2!=-1){
1274 if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
1275 else clusters = GetPHOSClusters() ;
1276 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1277 AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
1278 photon->GetMomentum(photonMom,GetVertex(0)) ;
1279 if(photon->GetID()==indexPhoton1) {
1280 ptDecay1 = photonMom.Pt();
1281 pxDecay1 = photonMom.Px();
1282 pyDecay1 = photonMom.Py();
1283 phiDecay1 = photonMom.Phi();
1284 }
1285 if(photon->GetID()==indexPhoton2) {
1286 ptDecay2 = photonMom.Pt();
1287 pxDecay2 = photonMom.Px();
1288 pyDecay2 = photonMom.Py();
1289 phiDecay2 = photonMom.Phi();
1290 }
1291 if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
1292 } //photonAOD loop
1293 } //index of decay photons found
1294 if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);
1295 } //make decay-hadron correlation
1296
1297 TObjArray * refpi0 =0x0;
1298 Int_t nrefs = 0;
1299
1300 //Loop on stored AOD pi0
1301 Int_t naod = pi0list->GetEntriesFast();
1302 if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - aod branch entries %d\n", naod);
1303 for(Int_t iaod = 0; iaod < naod ; iaod++){
1304 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (pi0list->At(iaod));
1305
1306 Int_t evtIndex2 = 0 ;
1307 Int_t evtIndex3 = 0 ;
1308 if (GetMixedEvent()) {
1309 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
1310 evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
1311
1312 if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
1313 continue ;
1314 }
1315
1316 //Int_t pdg = pi0->GetPdg();
1317 //if(pdg != AliCaloPID::kPi0) continue;
1318
1319 pt = pi0->Pt();
1320 px = pi0->Px();
1321 py = pi0->Py();
1322 if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
1323 //jumped out this event if near side associated partile pt larger than trigger
1324 if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;
1325
1326 //Selection within angular range
1327 phi = pi0->Phi();
1328 //Float_t deltaphi = TMath::Abs(phiTrig-phi);
1329 //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
1330
1331 if(bFillHisto){
1332
1333 deltaphi = phiTrig-phi;
1334 if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();
1335 if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
1336
1337 rat = pt/ptTrig ;
1338 phi = pi0->Phi() ;
1339 eta = pi0->Eta() ;
1340 xE = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1341 if(xE <0.)xE =-xE;
1342 cosi = TMath::Log(1/xE);
1343
1344 if(fPi0Trigger){
1345 if(indexPhoton1!=-1 && indexPhoton2!=-1){
1346 if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
1347 if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
1348 deltaphiDecay1 = phiDecay1-phi;
1349 deltaphiDecay2 = phiDecay2-phi;
1350 if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();
1351 if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();
1352 if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();
1353 if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();
1354 fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaphiDecay1);
1355 fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaphiDecay2);
1356 if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
1357 if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
1358 fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1);
1359 if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
1360 fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);
1361 if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
1362 }
1363 } //do decay-hadron correlation
1364
1365 fhEtaNeutral->Fill(pt,eta);
1366 fhPhiNeutral->Fill(pt,phi);
1367 fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
1368 fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);
1369 fhDeltaPhiDeltaEtaNeutral->Fill(deltaphi,etaTrig-eta);
1370
1371 //delta phi cut for correlation
1372 if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
1373 fhDeltaPhiNeutralPt->Fill(pt,deltaphi);
1374 fhPtImbalanceNeutral->Fill(ptTrig,rat);
1375 fhPtHbpNeutral->Fill(ptTrig,cosi);
1376 }
1377 else {
1378 fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);
1379 fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
1380 fhPtHbpUeNeutral->Fill(ptTrig,cosi);
1381 }
1382 //several UE calculation
1383 if(fMakeSeveralUE){
1384 if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
1385 fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);
1386 fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
1387 fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
1388 }
1389 if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
1390 fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);
1391 fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
1392 fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
1393 }
1394 } //several UE calculation
1395 }
1396 else{
1397 nrefs++;
1398 if(nrefs==1){
1399 refpi0 = new TObjArray(0);
1400 refpi0->SetName(GetAODObjArrayName()+"Pi0s");
1401 refpi0->SetOwner(kFALSE);
1402 }
1403 refpi0->Add(pi0);
1404 }//put references in trigger AOD
1405
1406 //if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
1407
1408 }//loop
1409}
1410
1411
1412//____________________________________________________________________________
1413//Bool_t AliAnaParticleHadronCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {
1414// //Select cluster depending on its pid and acceptance selections
1415//
1416// //Skip matched clusters with tracks
1417// if(IsTrackMatched(calo)) return kFALSE;
1418//
1419// TString detector = "";
1420// if (calo->IsPHOS()) detector= "PHOS";
1421// else if(calo->IsEMCAL()) detector= "EMCAL";
1422//
1423// //Check PID
1424// calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1425// pdg = AliCaloPID::kPhoton;
1426// if(IsCaloPIDOn()){
1427// //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1428// //or redo PID, recommended option for EMCal.
1429//
1430// if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1431// pdg = GetCaloPID()->GetPdg(detector,calo->GetPID(),mom.E());//PID with weights
1432// else
1433// pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
1434//
1435// if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
1436//
1437// //If it does not pass pid, skip
1438// if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {
1439// return kFALSE ;
1440// }
1441// }//PID on
1442//
1443// //Check acceptance selection
1444// if(IsFiducialCutOn()){
1445// Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,detector) ;
1446// if(! in ) return kFALSE ;
1447// }
1448//
1449// if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1450//
1451// return kTRUE;
1452//
1453//}