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