Add histograms to study shape of bad clusters
[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 **************************************************************************/
e09cf5ef 15
16//_________________________________________________________________________
17// Class for the analysis of particle - hadron correlations
18// Particle (for example direct gamma) must be found in a previous analysis
19//-- Author: Gustavo Conesa (LNF-INFN)
20
21// Modified by Yaxian Mao:
22// 1. add the UE subtraction for corrlation study
23// 2. change the correlation variable
24// 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)
25// 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)
26// 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06)
045396c8 27// 6. Add the possibility for event selection analysis based on vertex and multiplicity bins (10/10/2010)
28// 7. change the way of delta phi cut for UE study due to memory issue (reduce histograms)
29// 8. Add the possibility to request the absolute leading particle at the near side or not, set trigger bins, general clean-up (08/2011)
e09cf5ef 30//////////////////////////////////////////////////////////////////////////////
31
32
33// --- ROOT system ---
34//#include "TClonesArray.h"
35#include "TClass.h"
36#include "TMath.h"
37#include "TH3D.h"
045396c8 38#include "TDatabasePDG.h"
e09cf5ef 39
40//---- ANALYSIS system ----
41#include "AliNeutralMesonSelection.h"
42#include "AliAnaParticleHadronCorrelation.h"
43#include "AliCaloTrackReader.h"
44#include "AliCaloPID.h"
45#include "AliAODPWG4ParticleCorrelation.h"
46#include "AliFiducialCut.h"
47#include "AliVTrack.h"
48#include "AliVCluster.h"
49#include "AliMCAnalysisUtils.h"
50#include "TParticle.h"
51#include "AliStack.h"
52#include "AliAODMCParticle.h"
53#include "AliMixedEvent.h"
54
55ClassImp(AliAnaParticleHadronCorrelation)
56
57
05d0d05d 58//___________________________________________________________________
e09cf5ef 59 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
745913ae 60 AliAnaCaloTrackCorrBaseClass(),
66e64043 61 fMinTriggerPt(0.),
62 fMaxAssocPt(1000.), fMinAssocPt(0.),
05d0d05d 63 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
64 fSelectIsolated(0), fMakeSeveralUE(0),
65 fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.),
66 fPi0AODBranchName(""), fNeutralCorr(0),
67 fPi0Trigger(0), fMakeAbsoluteLeading(0),
68 fLeadingTriggerIndex(-1),
69 fNAssocPtBins(0), fAssocPtBinLimit(),
70 //Histograms
71 fhPtLeading(0), fhPhiLeading(0),
72 fhEtaLeading(0), fhDeltaPhiDeltaEtaCharged(0),
73 fhPhiCharged(0), fhEtaCharged(0),
74 fhDeltaPhiCharged(0), fhDeltaEtaCharged(0),
75 fhDeltaPhiChargedPt(0), fhDeltaPhiUeChargedPt(0),
76 fhPtImbalanceCharged(0), fhPtImbalanceUeCharged(0),
77 fhPtImbalancePosCharged(0), fhPtImbalanceNegCharged(0),
78 fhPtHbpCharged(0), fhPtHbpUeCharged(0),
79 fhDeltaPhiUeLeftCharged(0), fhDeltaPhiUeRightCharged(0),
80 fhPtImbalanceUeLeftCharged(0), fhPtImbalanceUeRightCharged(0),
81 fhPtHbpUeLeftCharged(0), fhPtHbpUeRightCharged(0),
82 fhPtTrigPout(0), fhPtTrigCharged(0),
83 fhTrigDeltaPhiCharged(0x0), fhTrigDeltaEtaCharged(0x0),
84 fhTrigCorr(0x0), fhTrigUeCorr(0x0),
85 fhAssocPt(0), fhAssocPtBkg(0),
66e64043 86 fhDeltaPhiAssocPtBin(0),
87 fhDeltaPhiBradAssocPtBin(0), fhDeltaPhiBrad(0),
88 fhXEAssocPtBin(0), fhXE(0),
e09cf5ef 89 fhDeltaPhiDeltaEtaNeutral(0),
05d0d05d 90 fhPhiNeutral(0), fhEtaNeutral(0),
91 fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),
92 fhDeltaPhiNeutralPt(0), fhDeltaPhiUeNeutralPt(0),
93 fhPtImbalanceNeutral(0), fhPtImbalanceUeNeutral(0),
94 fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),
95 fhDeltaPhiUeLeftNeutral(0), fhDeltaPhiUeRightNeutral(0),
96 fhPtImbalanceUeLeftNeutral(0), fhPtImbalanceUeRightNeutral(0),
97 fhPtHbpUeLeftNeutral(0), fhPtHbpUeRightNeutral(0),
e09cf5ef 98 fhPtPi0DecayRatio(0),
05d0d05d 99 fhDeltaPhiDecayCharged(0), fhPtImbalanceDecayCharged(0),
100 fhDeltaPhiDecayNeutral(0), fhPtImbalanceDecayNeutral(0),
045396c8 101 fh2phiLeadingParticle(0x0),
102 fhMCLeadingCount(0),
05d0d05d 103 fhMCEtaCharged(0), fhMCPhiCharged(0),
104 fhMCDeltaEtaCharged(0), fhMCDeltaPhiCharged(0x0),
105 fhMCDeltaPhiDeltaEtaCharged(0), fhMCDeltaPhiChargedPt(0),
045396c8 106 fhMCPtImbalanceCharged(0),
107 fhMCPtHbpCharged(0),
108 fhMCPtTrigPout(0),
109 fhMCPtAssocDeltaPhi(0)
e09cf5ef 110{
111 //Default Ctor
112
113 //Initialize parameters
114 InitParameters();
115}
116
05d0d05d 117//____________________________________________________________
118TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
119{
120 //Save parameters used for analysis
121 TString parList ; //this will be list of parameters used for this analysis.
122 const Int_t buffersize = 560;
123 char onePar[buffersize] ;
124
125 snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
126 parList+=onePar ;
66e64043 127 snprintf(onePar,buffersize," Pt Trigger > %3.2f ", fMinTriggerPt) ;
05d0d05d 128 parList+=onePar ;
66e64043 129 snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f ", fMinAssocPt, fMaxAssocPt) ;
05d0d05d 130 parList+=onePar ;
66e64043 131 snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f ", fDeltaPhiMinCut, fDeltaPhiMaxCut) ;
05d0d05d 132 parList+=onePar ;
66e64043 133 snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron < %3.2f ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
05d0d05d 134 parList+=onePar ;
66e64043 135 snprintf(onePar,buffersize,"Isolated Trigger? %d\n", fSelectIsolated) ;
05d0d05d 136 parList+=onePar ;
66e64043 137 snprintf(onePar,buffersize,"Several UE? %d\n", fMakeSeveralUE) ;
05d0d05d 138 parList+=onePar ;
66e64043 139 snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ", fPi0AODBranchName.Data());
05d0d05d 140 parList+=onePar ;
141 snprintf(onePar,buffersize,"Do Decay-hadron correlation ? %d", fPi0Trigger) ;
142 parList+=onePar ;
143 snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d\n", fMakeAbsoluteLeading) ;
144 parList+=onePar ;
145 snprintf(onePar,buffersize,"Associated particle pt bins %d: ", fNAssocPtBins) ;
146 parList+=onePar ;
147 for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
148 snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
149 }
150 parList+=onePar ;
151
152 //Get parameters set in base class.
153 parList += GetBaseParametersList() ;
154
155 //Get parameters set in PID class.
156 //parList += GetCaloPID()->GetPIDParametersList() ;
157
158 //Get parameters set in FiducialCut class (not available yet)
159 //parlist += GetFidCut()->GetFidCutParametersList()
160
161 return new TObjString(parList) ;
162
163}
164
165//________________________________________________________________
e09cf5ef 166TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
167{
168
169 // Create histograms to be saved in output file and
170 // store them in fOutputContainer
171 TList * outputContainer = new TList() ;
172 outputContainer->SetName("CorrelationHistos") ;
173
745913ae 174 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
175 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
176 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
e09cf5ef 177
e09cf5ef 178 fhPtLeading = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
179 fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
180
181 fhPhiLeading = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
182 fhPhiLeading->SetYTitle("#phi (rad)");
183
184 fhEtaLeading = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
185 fhEtaLeading->SetYTitle("#eta ");
045396c8 186
e09cf5ef 187 outputContainer->Add(fhPtLeading);
188 outputContainer->Add(fhPhiLeading);
189 outputContainer->Add(fhEtaLeading);
190
191 //Correlation with charged hadrons
192 if(GetReader()->IsCTSSwitchedOn()) {
193 fhDeltaPhiDeltaEtaCharged = new TH2F
194 ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
195 140,-2.,5.,200,-2,2);
196 fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
197 fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
198
199 fhPhiCharged = new TH2F
200 ("PhiCharged","#phi_{h^{#pm}} vs p_{T #pm}",
201 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
202 fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
203 fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
204
205 fhEtaCharged = new TH2F
206 ("EtaCharged","#eta_{h^{#pm}} vs p_{T #pm}",
207 nptbins,ptmin,ptmax,netabins,etamin,etamax);
208 fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
209 fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
210
211 fhDeltaPhiCharged = new TH2F
212 ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
213 nptbins,ptmin,ptmax,140,-2.,5.);
214 fhDeltaPhiCharged->SetYTitle("#Delta #phi");
215 fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
216
217 fhDeltaPhiChargedPt = new TH2F
218 ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
219 nptbins,ptmin,ptmax,140,-2.,5.);
220 fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
221 fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
222
223 fhDeltaPhiUeChargedPt = new TH2F
224 ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
225 nptbins,ptmin,ptmax,140,-2.,5.);
226 fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
227 fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
228
229 fhDeltaEtaCharged = new TH2F
230 ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
231 nptbins,ptmin,ptmax,200,-2,2);
232 fhDeltaEtaCharged->SetYTitle("#Delta #eta");
233 fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
234
235 fhPtImbalanceCharged =
236 new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
237 nptbins,ptmin,ptmax,200,0.,2.);
238 fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
239 fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
240
241 fhPtImbalanceUeCharged =
242 new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
243 nptbins,ptmin,ptmax,200,0.,2.);
244 fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
245 fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
246
247 fhPtImbalancePosCharged =
248 new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",
249 nptbins,ptmin,ptmax,200,0.,2.);
250 fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");
251 fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");
252
253 fhPtImbalanceNegCharged =
254 new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",
255 nptbins,ptmin,ptmax,200,0.,2.);
256 fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");
257 fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");
258
259 fhPtHbpCharged =
260 new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
261 nptbins,ptmin,ptmax,200,0.,10.);
262 fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
263 fhPtHbpCharged->SetXTitle("p_{T trigger}");
264
265 fhPtHbpUeCharged =
266 new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
267 nptbins,ptmin,ptmax,200,0.,10.);
268 fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
269 fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
270
271 fhPtTrigPout =
272 new TH2F("PtTrigPout","Pout with triggers",
273 nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
274 fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");
275 fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
276
e09cf5ef 277 fhPtTrigCharged =
278 new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",
279 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
280 fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");
281 fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
282
283 outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
284 outputContainer->Add(fhPhiCharged) ;
285 outputContainer->Add(fhEtaCharged) ;
286 outputContainer->Add(fhDeltaPhiCharged) ;
287 outputContainer->Add(fhDeltaEtaCharged) ;
288 outputContainer->Add(fhDeltaPhiChargedPt) ;
289 outputContainer->Add(fhDeltaPhiUeChargedPt) ;
290 outputContainer->Add(fhPtImbalanceCharged) ;
291 outputContainer->Add(fhPtImbalancePosCharged) ;
292 outputContainer->Add(fhPtImbalanceNegCharged) ;
293 outputContainer->Add(fhPtImbalanceUeCharged) ;
294 outputContainer->Add(fhPtHbpCharged) ;
295 outputContainer->Add(fhPtHbpUeCharged) ;
296 outputContainer->Add(fhPtTrigPout) ;
e09cf5ef 297 outputContainer->Add(fhPtTrigCharged) ;
298
299 if(DoEventSelect()){
300 Int_t nMultiBins = GetMultiBin();
301 fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;
302 fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;
05d0d05d 303 fhTrigCorr = new TH2F*[nMultiBins];
304 fhTrigUeCorr = new TH2F*[nMultiBins];
e09cf5ef 305 for(Int_t im=0; im<nMultiBins; im++){
306 fhTrigDeltaPhiCharged[im] = new TH2F
307 (Form("fhTrigDeltaPhiCharged_%d",im),Form("fhTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.);
308 fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
309 fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
310 fhTrigDeltaEtaCharged[im] = new TH2F
311 (Form("fhTrigDeltaEtaCharged_%d",im),Form("fhTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2);
312 fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
313 fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
314 fhTrigCorr[im] = new TH2F
315 (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
316 fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
317 fhTrigCorr[im]->SetXTitle("p_{T trigger}");
318 fhTrigUeCorr[im] = new TH2F
319 (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
320 fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
321 fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");
045396c8 322
e09cf5ef 323 outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;
324 outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;
325 outputContainer->Add(fhTrigCorr[im]);
326 outputContainer->Add(fhTrigUeCorr[im]);
327
045396c8 328 }
329 }
330
05d0d05d 331 fhAssocPt = new TH2F("fhAssocPt", " Trigger p_{T} vs associated hadron p_{T}",
332 nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
333 fhAssocPt->SetXTitle("p_{T trigger}");
334 fhAssocPt->SetYTitle("p_{T associated}");
335 outputContainer->Add(fhAssocPt) ;
336
337 fhAssocPtBkg = new TH2F("fhAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
338 nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
339 fhAssocPtBkg->SetXTitle("p_{T trigger}");
340 fhAssocPtBkg->SetYTitle("p_{T associated}");
341 outputContainer->Add(fhAssocPtBkg) ;
342
66e64043 343 fhDeltaPhiBrad = new TH2F("fhDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ",
344 nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
345 fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
346 fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
347
348 outputContainer->Add(fhDeltaPhiBrad) ;
349
350 fhXE = new TH2F("fhXE", "x_{E} vs p_{T trigger}",
351 nptbins, ptmin, ptmax,50, 0.0, 2.0);
352 fhXE->SetXTitle("p_{T trigger}");
353 fhXE->SetYTitle("x_{E}");
354
355 outputContainer->Add(fhXE);
356
357
05d0d05d 358 fhDeltaPhiAssocPtBin = new TH2F*[fNAssocPtBins] ;
359 fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
360 fhXEAssocPtBin = new TH2F*[fNAssocPtBins] ;
361 for(Int_t i = 0 ; i < fNAssocPtBins ; i++){
362 fhDeltaPhiAssocPtBin[i] = new TH2F(Form("fhDeltaPhiPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
363 Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
364 nptbins, ptmin, ptmax,140,-2.,5.);
365 fhDeltaPhiAssocPtBin[i]->SetXTitle("p_{T trigger}");
366 fhDeltaPhiAssocPtBin[i]->SetYTitle("#Delta #phi");
367
368 fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("fhDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
369 Form("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
370 nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
371 fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
372 fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
373
374
375 fhXEAssocPtBin[i] = new TH2F(Form("fhXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
66e64043 376 Form("x_{E} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
05d0d05d 377 nptbins, ptmin, ptmax,50, 0.0, 2.0);
378 fhXEAssocPtBin[i]->SetXTitle("p_{T trigger}");
379 fhXEAssocPtBin[i]->SetYTitle("x_{E}");
380
381 outputContainer->Add(fhDeltaPhiAssocPtBin[i]) ;
382 outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
383 outputContainer->Add(fhXEAssocPtBin[i]);
e09cf5ef 384 }
045396c8 385
e09cf5ef 386 if(fPi0Trigger){
387 fhPtPi0DecayRatio = new TH2F
388 ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
389 nptbins,ptmin,ptmax, 100,0.,2.);
390 fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
391 fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
392
393 fhDeltaPhiDecayCharged = new TH2F
394 ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
395 nptbins,ptmin,ptmax,140,-2.,5.);
396 fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
397 fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
398
399 fhPtImbalanceDecayCharged =
400 new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
401 nptbins,ptmin,ptmax,200,0.,2.);
402 fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");
403 fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");
404
405 outputContainer->Add(fhPtPi0DecayRatio) ;
406 outputContainer->Add(fhDeltaPhiDecayCharged) ;
407 outputContainer->Add(fhPtImbalanceDecayCharged) ;
408 }
409
410
411 if(fMakeSeveralUE){
412 fhDeltaPhiUeLeftCharged = new TH2F
413 ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
414 nptbins,ptmin,ptmax,140,-2.,5.);
415 fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
416 fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
417 outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
418
419 fhDeltaPhiUeRightCharged = new TH2F
420 ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
421 nptbins,ptmin,ptmax,140,-2.,5.);
422 fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
423 fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
424 outputContainer->Add(fhDeltaPhiUeRightCharged) ;
425
426 fhPtImbalanceUeLeftCharged =
427 new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
428 nptbins,ptmin,ptmax,200,0.,2.);
429 fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
430 fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
431 outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
432
433 fhPtImbalanceUeRightCharged =
434 new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
435 nptbins,ptmin,ptmax,200,0.,2.);
436 fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
437 fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
438 outputContainer->Add(fhPtImbalanceUeRightCharged) ;
439
440 fhPtHbpUeLeftCharged =
441 new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
442 nptbins,ptmin,ptmax,200,0.,10.);
443 fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
444 fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
445 outputContainer->Add(fhPtHbpUeLeftCharged) ;
446
447 fhPtHbpUeRightCharged =
448 new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
449 nptbins,ptmin,ptmax,200,0.,10.);
450 fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
451 fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
452 outputContainer->Add(fhPtHbpUeRightCharged) ;
453
454 }
455 } //Correlation with charged hadrons
456
457 //Correlation with neutral hadrons
458 if(fNeutralCorr){
459
460 fhDeltaPhiDeltaEtaNeutral = new TH2F
461 ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
462 140,-2.,5.,200,-2,2);
463 fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
464 fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
465
466 fhPhiNeutral = new TH2F
467 ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #pi^{0}}",
468 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
469 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
470 fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
471
472 fhEtaNeutral = new TH2F
473 ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #pi^{0}}",
474 nptbins,ptmin,ptmax,netabins,etamin,etamax);
475 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
476 fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
477
478 fhDeltaPhiNeutral = new TH2F
479 ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
480 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
481 fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
482 fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
483
484 fhDeltaPhiNeutralPt = new TH2F
485 ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
486 nptbins,ptmin,ptmax,140,-2.,5.);
487 fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
488 fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
489
490 fhDeltaPhiUeNeutralPt = new TH2F
491 ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
492 nptbins,ptmin,ptmax,140,-2.,5.);
493 fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
494 fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
495
496 fhDeltaEtaNeutral = new TH2F
497 ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
498 nptbins,ptmin,ptmax,200,-2,2);
499 fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
500 fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
501
502 fhPtImbalanceNeutral =
503 new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
504 nptbins,ptmin,ptmax,200,0.,2.);
505 fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
506 fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
507
508 fhPtImbalanceUeNeutral =
509 new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
510 nptbins,ptmin,ptmax,200,0.,2.);
511 fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
512 fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
513
514 fhPtHbpNeutral =
515 new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
516 nptbins,ptmin,ptmax,200,0.,10.);
517 fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
518 fhPtHbpNeutral->SetXTitle("p_{T trigger}");
519
520 fhPtHbpUeNeutral =
521 new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
522 nptbins,ptmin,ptmax,200,0.,10.);
523 fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
524 fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
525
526
527 outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
528 outputContainer->Add(fhPhiNeutral) ;
529 outputContainer->Add(fhEtaNeutral) ;
530 outputContainer->Add(fhDeltaPhiNeutral) ;
531 outputContainer->Add(fhDeltaPhiNeutralPt) ;
532 outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
533 outputContainer->Add(fhDeltaEtaNeutral) ;
534 outputContainer->Add(fhPtImbalanceNeutral) ;
535 outputContainer->Add(fhPtImbalanceUeNeutral) ;
536 outputContainer->Add(fhPtHbpNeutral) ;
537 outputContainer->Add(fhPtHbpUeNeutral) ;
538
539 if(fPi0Trigger){
540 fhDeltaPhiDecayNeutral = new TH2F
541 ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
542 nptbins,ptmin,ptmax,140,-2.,5.);
543 fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
544 fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
545
546 fhPtImbalanceDecayNeutral =
547 new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
548 nptbins,ptmin,ptmax,200,0.,2.);
549 fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");
550 fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");
551
552 outputContainer->Add(fhDeltaPhiDecayNeutral) ;
553 outputContainer->Add(fhPtImbalanceDecayNeutral) ;
554 }
555
556
557 if(fMakeSeveralUE){
558 fhDeltaPhiUeLeftNeutral = new TH2F
559 ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
560 nptbins,ptmin,ptmax,140,-2.,5.);
561 fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
562 fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
563 outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
564
565 fhDeltaPhiUeRightNeutral = new TH2F
566 ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
567 nptbins,ptmin,ptmax,140,-2.,5.);
568 fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
569 fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
570 outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
571
572 fhPtImbalanceUeLeftNeutral =
573 new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
574 nptbins,ptmin,ptmax,140,0.,2.);
575 fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
576 fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
577 outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
578
579 fhPtImbalanceUeRightNeutral =
580 new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
581 nptbins,ptmin,ptmax,200,0.,2.);
582 fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
583 fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
584 outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
585
586 fhPtHbpUeLeftNeutral =
587 new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
588 nptbins,ptmin,ptmax,200,0.,10.);
589 fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
590 fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
591 outputContainer->Add(fhPtHbpUeLeftNeutral) ;
592
593 fhPtHbpUeRightNeutral =
594 new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
595 nptbins,ptmin,ptmax,200,0.,10.);
596 fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
597 fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
598 outputContainer->Add(fhPtHbpUeRightNeutral) ;
599
600 }
fe871b21 601
602 }//Correlation with neutral hadrons
603
604 //if data is MC, fill more histograms
605 if(IsDataMC()){
e09cf5ef 606
fe871b21 607 fh2phiLeadingParticle=new TH2F("fh2phiLeadingParticle","#phi resolustion for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
608 fh2phiLeadingParticle->GetXaxis()->SetTitle("p_{T gen Leading} (GeV/c)");
609 fh2phiLeadingParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
e09cf5ef 610
fe871b21 611 fhMCLeadingCount=new TH1F("MCLeadingTriggerCount","MCLeadingTriggerCount",nptbins,ptmin,ptmax);
612 fhMCLeadingCount->SetXTitle("p_{T trig}");
613
614 fhMCEtaCharged = new TH2F
615 ("MCEtaCharged","MC #eta_{h^{#pm}} vs p_{T #pm}",
616 nptbins,ptmin,ptmax,netabins,etamin,etamax);
617 fhMCEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
618 fhMCEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
619
620 fhMCPhiCharged = new TH2F
621 ("MCPhiCharged","#MC phi_{h^{#pm}} vs p_{T #pm}",
622 200,ptmin,ptmax,nphibins,phimin,phimax);
623 fhMCPhiCharged->SetYTitle("MC #phi_{h^{#pm}} (rad)");
624 fhMCPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
625
626 fhMCDeltaPhiDeltaEtaCharged = new TH2F
627 ("MCDeltaPhiDeltaEtaCharged","#MC phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
628 140,-2.,5.,200,-2,2);
629 fhMCDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
630 fhMCDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
631
632 fhMCDeltaEtaCharged = new TH2F
633 ("MCDeltaEtaCharged","MC #eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger} and p_{T assoc}",
634 nptbins,ptmin,ptmax,200,-2,2);
635 fhMCDeltaEtaCharged->SetYTitle("#Delta #eta");
636 fhMCDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
637
638 fhMCDeltaPhiCharged = new TH2F
639 ("MCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
640 nptbins,ptmin,ptmax,140,-2.,5.);
641 fhMCDeltaPhiCharged->SetYTitle("#Delta #phi");
642 fhMCDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
643
644 fhMCDeltaPhiChargedPt = new TH2F
645 ("MCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
646 nptbins,ptmin,ptmax,140,-2.,5.);
647 fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi");
648 fhMCDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
649
650 fhMCPtImbalanceCharged =
651 new TH2F("MCCorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
652 nptbins,ptmin,ptmax,200,0.,2.);
653 fhMCPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
654 fhMCPtImbalanceCharged->SetXTitle("p_{T trigger}");
655
656 fhMCPtHbpCharged =
657 new TH2F("MCHbpCharged","MC #xi = ln(1/x_{E}) with charged hadrons",
658 nptbins,ptmin,ptmax,200,0.,10.);
659 fhMCPtHbpCharged->SetYTitle("ln(1/x_{E})");
660 fhMCPtHbpCharged->SetXTitle("p_{T trigger}");
661
662 fhMCPtTrigPout =
663 new TH2F("MCPtTrigPout","AOD MC Pout with triggers",
664 nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
665 fhMCPtTrigPout->SetYTitle("p_{out} (GeV/c)");
666 fhMCPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
667
668 fhMCPtAssocDeltaPhi =
669 new TH2F("fhMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
670 nptbins,ptmin,ptmax,140,-2.,5.);
671 fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi");
672 fhMCPtAssocDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
673
674 outputContainer->Add(fh2phiLeadingParticle);
675 outputContainer->Add(fhMCLeadingCount);
676 outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
677 outputContainer->Add(fhMCPhiCharged) ;
678 outputContainer->Add(fhMCEtaCharged) ;
679 outputContainer->Add(fhMCDeltaEtaCharged) ;
680 outputContainer->Add(fhMCDeltaPhiCharged) ;
681
682 outputContainer->Add(fhMCDeltaPhiChargedPt) ;
683 outputContainer->Add(fhMCPtImbalanceCharged) ;
684 outputContainer->Add(fhMCPtHbpCharged) ;
685 outputContainer->Add(fhMCPtTrigPout) ;
686 outputContainer->Add(fhMCPtAssocDeltaPhi) ;
687 } //for MC histogram
688
e09cf5ef 689
690 return outputContainer;
691
692}
693
05d0d05d 694//____________________________________________________
e09cf5ef 695void AliAnaParticleHadronCorrelation::InitParameters()
696{
697
698 //Initialize the parameters of the analysis.
699 SetInputAODName("PWG4Particle");
700 SetAODObjArrayName("Hadrons");
701 AddToHistogramsName("AnaHadronCorr_");
045396c8 702
e09cf5ef 703 SetPtCutRange(0.,300);
045396c8 704 fDeltaPhiMinCut = 1.5 ;
705 fDeltaPhiMaxCut = 4.5 ;
706 fSelectIsolated = kFALSE;
707 fMakeSeveralUE = kFALSE;
708 fUeDeltaPhiMinCut = 1. ;
709 fUeDeltaPhiMaxCut = 1.5 ;
710 fNeutralCorr = kFALSE ;
711 fPi0Trigger = kFALSE ;
712 fMakeAbsoluteLeading = kTRUE;
713
66e64043 714 fNAssocPtBins = 7 ;
715 fAssocPtBinLimit[0] = 2. ;
716 fAssocPtBinLimit[1] = 4. ;
717 fAssocPtBinLimit[2] = 6. ;
718 fAssocPtBinLimit[3] = 8. ;
719 fAssocPtBinLimit[4] = 10. ;
720 fAssocPtBinLimit[5] = 13. ;
721 fAssocPtBinLimit[6] = 16. ;
722 fAssocPtBinLimit[7] = 20. ;
723 fAssocPtBinLimit[8] = 100.;
724 fAssocPtBinLimit[9] = 200.;
045396c8 725
e09cf5ef 726}
727
05d0d05d 728//__________________________________________________________
e09cf5ef 729void AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
730{
731 //Particle-Hadron Correlation Analysis, fill AODs
732
733 if(!GetInputAODBranch()){
734 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
735 abort();
736 }
737
738 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
739 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());
740 abort();
741 }
742
743 if(GetDebug() > 1){
744 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
745 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
66e64043 746 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks() ->GetEntriesFast());
e09cf5ef 747 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
66e64043 748 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetPHOSClusters() ->GetEntriesFast());
e09cf5ef 749 }
750
045396c8 751 //Get the vertex and check it is not too large in z
752 Double_t v[3] = {0,0,0}; //vertex ;
753 GetReader()->GetVertex(v);
754 if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
755
756 //Loop on stored AOD particles, find leading trigger
66e64043 757 Double_t ptTrig = fMinTriggerPt ;
045396c8 758 fLeadingTriggerIndex = -1 ;
759 Int_t naod = GetInputAODBranch()->GetEntriesFast() ;
e09cf5ef 760 for(Int_t iaod = 0; iaod < naod ; iaod++){
761 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
045396c8 762
763 //vertex cut in case of mixing
764 if (GetMixedEvent()) {
765 Int_t evt=-1;
766 Int_t id =-1;
767 if (particle->GetCaloLabel(0) >= 0 ){
768 id=particle->GetCaloLabel(0);
769 if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
770 }
771 else if(particle->GetTrackLabel(0) >= 0 ){
772 id=particle->GetTrackLabel(0);
773 if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
774 }
775 else continue;
776
777 if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
778 return ;
e09cf5ef 779 }
045396c8 780
781 // find the leading particles with highest momentum
782 if (particle->Pt() > ptTrig) {
783 ptTrig = particle->Pt() ;
784 fLeadingTriggerIndex = iaod ;
785 }
786 }// finish search of leading trigger particle
e09cf5ef 787
66e64043 788
e09cf5ef 789 //Do correlation with leading particle
045396c8 790 if(fLeadingTriggerIndex >= 0){
e09cf5ef 791
045396c8 792 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
793
794 //check if the particle is isolated or if we want to take the isolation into account
795 if(OnlyIsolated() && !particle->IsIsolated()) return;
796
e09cf5ef 797 //Make correlation with charged hadrons
045396c8 798 Bool_t okcharged = kTRUE;
799 Bool_t okneutral = kTRUE;
e09cf5ef 800 if(GetReader()->IsCTSSwitchedOn() )
045396c8 801 okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);
e09cf5ef 802
803 TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
804 if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
045396c8 805 okneutral = MakeNeutralCorrelation(particle, pi0list,kFALSE);
e09cf5ef 806
807 }//Correlate leading
808
809 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
810
811}
812
05d0d05d 813//_________________________________________________________________
e09cf5ef 814void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
815{
816 //Particle-Hadron Correlation Analysis, fill histograms
817
818 if(!GetInputAODBranch()){
819 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
820 abort();
821 }
822
823 if(GetDebug() > 1){
824 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
825 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
826 }
827
828 //Get the vertex and check it is not too large in z
829 Double_t v[3] = {0,0,0}; //vertex ;
830 GetReader()->GetVertex(v);
831 if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
832
833 //Loop on stored AOD particles, find leading
66e64043 834 Double_t ptTrig = fMinTriggerPt;
045396c8 835 if(fLeadingTriggerIndex < 0){//Search leading if not done before
836 Int_t naod = GetInputAODBranch()->GetEntriesFast() ;
837 for(Int_t iaod = 0; iaod < naod ; iaod++){ //loop on input trigger AOD file
838 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
839 //vertex cut in case of mixing
840 if (GetMixedEvent()) {
841 Int_t evt=-1;
842 Int_t id =-1;
843 if (particle->GetCaloLabel(0) >= 0 ){
844 id=particle->GetCaloLabel(0);
845 if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
846 }
847 else if(particle->GetTrackLabel(0) >= 0 ){
848 id=particle->GetTrackLabel(0);
849 if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
850 }
851 else continue;
852
853 if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
854 return ;
e09cf5ef 855 }
05d0d05d 856
045396c8 857 //check if the particle is isolated or if we want to take the isolation into account
858 if(OnlyIsolated() && !particle->IsIsolated()) continue;
e09cf5ef 859
045396c8 860 //check if inside the vertex cut
861 //find the leading particles with highest momentum
862 if (particle->Pt() > ptTrig) {
863 ptTrig = particle->Pt() ;
864 fLeadingTriggerIndex = iaod ;
e09cf5ef 865 }
045396c8 866 }//finish search of leading trigger particle
867 }//Search leading if not done before
868
869 if(fLeadingTriggerIndex >= 0 ){ //using trigger particle to do correlations
870 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
871
872 //check if the particle is isolated or if we want to take the isolation into account
873 if(OnlyIsolated() && !particle->IsIsolated()) return;
e09cf5ef 874
875 //Make correlation with charged hadrons
045396c8 876 Bool_t okcharged = kTRUE;
877 Bool_t okneutral = kTRUE;
878 if(GetReader()->IsCTSSwitchedOn() ){
879 okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);
a5f1e836 880 if(IsDataMC()){
045396c8 881 MakeMCChargedCorrelation(particle);
882 }
a5f1e836 883 }
884
e09cf5ef 885 TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
a5f1e836 886 if(fNeutralCorr && pi0list){
887 if(pi0list->GetEntriesFast() > 0)
888 okneutral = MakeNeutralCorrelation(particle, pi0list,kTRUE);
889 }
e09cf5ef 890
045396c8 891 // Fill leading particle histogram if correlation went well and
892 // no problem was found, like not absolute leading, or bad vertex in mixing.
893 if(okcharged && okneutral){
894 fhPtLeading->Fill(particle->Pt());
895 Float_t phi = particle->Phi();
896 if(phi<0)phi+=TMath::TwoPi();
897 fhPhiLeading->Fill(particle->Pt(), phi);
898 fhEtaLeading->Fill(particle->Pt(), particle->Eta());
899 }//ok charged && neutral
e09cf5ef 900 }//Aod branch loop
901
045396c8 902 //Reinit for next event
903 fLeadingTriggerIndex = -1;
e09cf5ef 904
045396c8 905 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
e09cf5ef 906}
907
c5693f62 908//___________________________________________________________________________________________________________
045396c8 909Bool_t AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle,
c5693f62 910 const TObjArray* pl, const Bool_t bFillHisto)
e09cf5ef 911{
912 // Charged Hadron Correlation Analysis
913 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
914
045396c8 915 Int_t evtIndex11 = -1 ; //cluster trigger or pi0 trigger
916 Int_t evtIndex12 = -1 ; // pi0 trigger
917 Int_t evtIndex13 = -1 ; // charged trigger
e09cf5ef 918 Int_t indexPhoton1 = -1 ;
919 Int_t indexPhoton2 = -1 ;
045396c8 920
921 Double_t v[3] = {0,0,0}; //vertex ;
e09cf5ef 922 GetReader()->GetVertex(v);
e09cf5ef 923
924 if (GetMixedEvent()) {
925 evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
926 evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
927 evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
928 }
929
e09cf5ef 930 Double_t phiTrig = aodParticle->Phi();
045396c8 931 Double_t etaTrig = aodParticle->Eta();
932 Double_t ptTrig = aodParticle->Pt();
045396c8 933
934 Double_t pt = -100. ;
935 Double_t px = -100. ;
936 Double_t py = -100. ;
937 Double_t rat = -100. ;
938 Double_t xE = -100. ;
939 Double_t cosi = -100. ;
940 Double_t phi = -100. ;
941 Double_t eta = -100. ;
942 Double_t pout = -100. ;
e09cf5ef 943
045396c8 944 Double_t ptDecay1 = 0. ;
945 Double_t pxDecay1 = 0. ;
946 Double_t pyDecay1 = 0. ;
947 Double_t phiDecay1 = 0. ;
948 Double_t ptDecay2 = 0. ;
949 Double_t pxDecay2 = 0. ;
950 Double_t pyDecay2 = 0. ;
951 Double_t phiDecay2 = 0. ;
952
953 Double_t ratDecay1 = -100. ;
954 Double_t ratDecay2 = -100. ;
955 Double_t deltaPhi = -100. ;
956 Double_t deltaPhiOrg = -100. ;
957 Double_t deltaPhiDecay1 = -100. ;
958 Double_t deltaPhiDecay2 = -100. ;
959
960 TVector3 p3;
e09cf5ef 961 TLorentzVector photonMom ;
045396c8 962 TObjArray * clusters = 0x0 ;
963 TObjArray * reftracks = 0x0;
964 Int_t nrefs = 0;
965 Int_t nTracks = GetCTSTracks()->GetEntriesFast() ;
e09cf5ef 966
967 if(fPi0Trigger){
968 indexPhoton1 = aodParticle->GetCaloLabel (0);
969 indexPhoton2 = aodParticle->GetCaloLabel (1);
970 if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
971
972 if(indexPhoton1!=-1 && indexPhoton2!=-1){
973 if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
045396c8 974 else clusters = GetPHOSClusters() ;
e09cf5ef 975 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
976 AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
977 photon->GetMomentum(photonMom,GetVertex(0)) ;
978 if(photon->GetID()==indexPhoton1) {
979 ptDecay1 = photonMom.Pt();
980 pxDecay1 = photonMom.Px();
981 pyDecay1 = photonMom.Py();
982 phiDecay1 = photonMom.Phi();
983 if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);
984 }
985 if(photon->GetID()==indexPhoton2) {
986 ptDecay2 = photonMom.Pt();
987 pxDecay2 = photonMom.Px();
988 pyDecay2 = photonMom.Py();
989 phiDecay2 = photonMom.Phi();
990 if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);
991 }
992 if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
993 } //cluster loop
994 } //index of decay photons found
e09cf5ef 995 } //make decay-hadron correlation
996
997 //Track loop, select tracks with good pt, phi and fill AODs or histograms
e09cf5ef 998 for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
999 AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
045396c8 1000
e09cf5ef 1001 Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
1002 p3.SetXYZ(mom[0],mom[1],mom[2]);
1003 pt = p3.Pt();
1004 px = p3.Px();
1005 py = p3.Py();
1006 eta = p3.Eta();
1007 phi = p3.Phi() ;
1008 if(phi < 0) phi+=TMath::TwoPi();
045396c8 1009
e09cf5ef 1010 //Select only hadrons in pt range
66e64043 1011 if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1012
e09cf5ef 1013 //remove trigger itself for correlation when use charged triggers
1014 if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
66e64043 1015 track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3) )
045396c8 1016 continue ;
1017
e09cf5ef 1018 if(IsFiducialCutOn()){
1019 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
1020 if(! in ) continue ;
1021 }
045396c8 1022
1023 //jump out this event if near side associated particle pt larger than trigger
1024 if (fMakeAbsoluteLeading){
1025 if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) return kFALSE;
1026 }
1027
1028 //Only for mixed event
1029 Int_t evtIndex2 = 0 ;
1030 if (GetMixedEvent()) {
1031 evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
1032 if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
1033 continue ;
1034 //vertex cut
1035 if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
1036 return kFALSE;
1037 }
1038
1039 if(fPi0Trigger){
e09cf5ef 1040 if(indexPhoton1!=-1 && indexPhoton2!=-1){
1041 if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
1042 if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
045396c8 1043 deltaPhiDecay1 = phiDecay1-phi;
1044 deltaPhiDecay2 = phiDecay2-phi;
1045 if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
1046 if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
1047 if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
1048 if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
e09cf5ef 1049 }
1050 } //do decay-hadron correlation
045396c8 1051
e09cf5ef 1052 //Selection within angular range
045396c8 1053 deltaPhi = phiTrig-phi;
1054 deltaPhiOrg = deltaPhi;
1055 if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1056 if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
1057
045396c8 1058 rat = pt/ptTrig ;
66e64043 1059 pout = pt*TMath::Sin(deltaPhi) ;
1060 xE =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1061 //if(xE <0.)xE =-xE;
1062 cosi = -100;
1063 if(xE > 0 ) cosi = TMath::Log(1./xE);
1064
e09cf5ef 1065 if(GetDebug() > 2)
66e64043 1066 printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f, pT Trig min %2.2f \n",
1067 pt,phi, phiTrig,fDeltaPhiMinCut, deltaPhi, fDeltaPhiMaxCut, fMinTriggerPt);
e09cf5ef 1068
045396c8 1069 // Fill Histograms
e09cf5ef 1070 if(bFillHisto){
045396c8 1071
05d0d05d 1072 Int_t assocBin = -1;
1073 for(Int_t i = 0 ; i < fNAssocPtBins ; i++){
66e64043 1074 if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
045396c8 1075 }
1076
05d0d05d 1077 fhAssocPt->Fill(ptTrig,pt);
1078
66e64043 1079 fhXE->Fill(ptTrig, xE) ;
1080
05d0d05d 1081 if(TMath::Cos(deltaPhi) < 0 && assocBin >= 0 )//away side
1082 fhXEAssocPtBin[assocBin]->Fill(ptTrig, xE) ;
1083
1084 //Hardcoded values, BAD, FIXME
1085 Double_t dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
1086 if(TMath::Abs(dphiBrad)>0.325 && TMath::Abs(dphiBrad)<0.475){
1087 fhAssocPtBkg->Fill(ptTrig, pt);
1088 }
1089
1090 if(dphiBrad<-1./3) dphiBrad += 2;
66e64043 1091 fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
1092
05d0d05d 1093 if(assocBin>=0){
1094 fhDeltaPhiBradAssocPtBin[assocBin]->Fill(ptTrig, dphiBrad);
1095 fhDeltaPhiAssocPtBin [assocBin]->Fill(ptTrig, deltaPhi);
1096 }
1097
1098 fhEtaCharged ->Fill(pt,eta);
1099 fhPhiCharged ->Fill(pt,phi);
e09cf5ef 1100 fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
045396c8 1101 fhDeltaPhiCharged->Fill(ptTrig, deltaPhi);
745913ae 1102 if(pt >2 ) fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi,aodParticle->Eta()-eta);
e09cf5ef 1103
1104 if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
1105 //fill different multiplicity histogram
1106 if(DoEventSelect()){
1107 for(Int_t im=0; im<GetMultiBin(); im++){
1108 if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1)){
045396c8 1109 fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaPhi);
e09cf5ef 1110 fhTrigDeltaEtaCharged[im]->Fill(ptTrig,aodParticle->Eta()-eta);
045396c8 1111 }
e09cf5ef 1112 }
1113 }
66e64043 1114
e09cf5ef 1115 //delta phi cut for correlation
66e64043 1116 if ( (deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) ) {
1117 fhDeltaPhiChargedPt ->Fill(pt,deltaPhi);
e09cf5ef 1118 fhPtImbalanceCharged->Fill(ptTrig,xE);
05d0d05d 1119 fhPtHbpCharged ->Fill(ptTrig, cosi);
1120 fhPtTrigPout ->Fill(ptTrig, pout) ;
e09cf5ef 1121 fhPtTrigCharged->Fill(ptTrig, pt) ;
05d0d05d 1122 if(track->Charge() > 0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;
1123 else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;
e09cf5ef 1124 //fill different multiplicity histogram
1125 if(DoEventSelect()){
1126 for(Int_t im=0; im<GetMultiBin(); im++){
1127 if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
1128 fhTrigCorr[im]->Fill(ptTrig,xE);
1129 }
1130 } //multiplicity events selection
1131 } //delta phi cut for correlation
66e64043 1132 else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) ) { //UE study
045396c8 1133 fhDeltaPhiUeChargedPt->Fill(pt,deltaPhi);
1134 Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
1135 Double_t uexE = -(pt/ptTrig)*TMath::Cos(randomphi);
1136 if(uexE < 0.) uexE = -uexE;
1137 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - xe = %f, uexE = %f \n", xE, uexE);
1138 fhPtImbalanceUeCharged->Fill(ptTrig,uexE);
c514f286 1139 if(uexE>0)fhPtHbpUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
e09cf5ef 1140 if(DoEventSelect()){
1141 for(Int_t im=0; im<GetMultiBin(); im++){
1142 if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
1143 fhTrigUeCorr[im]->Fill(ptTrig,xE);
1144 }
1145 } //multiplicity events selection
1146
1147 } //UE study
1148
1149 if(fPi0Trigger){
1150 if(indexPhoton1!=-1 && indexPhoton2!=-1){
045396c8 1151 fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
1152 fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
1153 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
1154 if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
e09cf5ef 1155 fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1);
045396c8 1156 if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
e09cf5ef 1157 fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);
045396c8 1158 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
e09cf5ef 1159 } //index of decay photons found
1160 } //make decay-hadron correlation
1161
1162 //several UE calculation
1163 if(fMakeSeveralUE){
045396c8 1164 if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut)){
1165 fhDeltaPhiUeLeftCharged->Fill(pt,deltaPhi);
e09cf5ef 1166 fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
1167 fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
1168 }
045396c8 1169 if((deltaPhi>fUeDeltaPhiMinCut) && (deltaPhi <fUeDeltaPhiMaxCut)){
1170 fhDeltaPhiUeRightCharged->Fill(pt,deltaPhi);
e09cf5ef 1171 fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
1172 fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
1173
1174 }
1175 } //several UE calculation
1176
045396c8 1177 //Fill leading particle histogram
1178 fhPtLeading->Fill(ptTrig);
1179 if(phiTrig<0)phiTrig+=TMath::TwoPi();
1180 fhPhiLeading->Fill(ptTrig, phiTrig);
1181 fhEtaLeading->Fill(ptTrig, etaTrig);
1182
e09cf5ef 1183 } //Fill histogram
1184 else{
1185 nrefs++;
1186 if(nrefs==1){
1187 reftracks = new TObjArray(0);
045396c8 1188 TString trackname = Form("%s+Tracks", GetAODObjArrayName().Data());
1189 reftracks->SetName(trackname.Data());
e09cf5ef 1190 reftracks->SetOwner(kFALSE);
1191 }
1192 reftracks->Add(track);
1193 }//aod particle loop
1194 }// track loop
1195
1196 //Fill AOD with reference tracks, if not filling histograms
1197 if(!bFillHisto && reftracks) {
1198 aodParticle->AddObjArray(reftracks);
1199 }
1200
045396c8 1201 return kTRUE;
e09cf5ef 1202
1203}
1204
05d0d05d 1205//________________________________________________________________________________________________________________
045396c8 1206Bool_t AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle,
c5693f62 1207 const TObjArray* pi0list, const Bool_t bFillHisto)
e09cf5ef 1208{
1209 // Neutral Pion Correlation Analysis
1210 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());
1211
1212 Int_t evtIndex11 = 0 ;
1213 Int_t evtIndex12 = 0 ;
1214 if (GetMixedEvent()) {
1215 evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
1216 evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
1217 }
1218
045396c8 1219 Double_t pt = -100. ;
1220 Double_t px = -100. ;
1221 Double_t py = -100. ;
1222 Double_t rat = -100. ;
1223 Double_t phi = -100. ;
1224 Double_t eta = -100. ;
1225 Double_t xE = -100. ;
1226 Double_t cosi = -100. ;
e09cf5ef 1227
1228 Double_t ptTrig = aodParticle->Pt();
1229 Double_t phiTrig = aodParticle->Phi();
1230 Double_t etaTrig = aodParticle->Eta();
66e64043 1231 //Double_t pxTrig = aodParticle->Px();
1232 //Double_t pyTrig = aodParticle->Py();
e09cf5ef 1233
045396c8 1234 Int_t indexPhoton1 =-1 ;
1235 Int_t indexPhoton2 =-1 ;
1236 Double_t ptDecay1 = 0. ;
e09cf5ef 1237 Double_t pxDecay1 = 0. ;
1238 Double_t pyDecay1 = 0. ;
045396c8 1239 Double_t phiDecay1 = 0. ;
1240 Double_t ptDecay2 = 0. ;
e09cf5ef 1241 Double_t pxDecay2 = 0. ;
1242 Double_t pyDecay2 = 0. ;
045396c8 1243 Double_t phiDecay2 = 0. ;
1244
1245 Double_t ratDecay1 = -100. ;
1246 Double_t ratDecay2 = -100. ;
1247 Double_t deltaPhi = -100. ;
1248 Double_t deltaPhiDecay1 = -100. ;
1249 Double_t deltaPhiDecay2 = -100. ;
e09cf5ef 1250
e09cf5ef 1251 TObjArray * clusters = 0x0 ;
045396c8 1252 TLorentzVector photonMom ;
1253
e09cf5ef 1254 if(fPi0Trigger){
1255 indexPhoton1 = aodParticle->GetCaloLabel (0);
1256 indexPhoton2 = aodParticle->GetCaloLabel (1);
045396c8 1257 if(GetDebug() > 1)
1258 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
e09cf5ef 1259
1260 if(indexPhoton1!=-1 && indexPhoton2!=-1){
1261 if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
1262 else clusters = GetPHOSClusters() ;
1263 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1264 AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
1265 photon->GetMomentum(photonMom,GetVertex(0)) ;
1266 if(photon->GetID()==indexPhoton1) {
1267 ptDecay1 = photonMom.Pt();
1268 pxDecay1 = photonMom.Px();
1269 pyDecay1 = photonMom.Py();
1270 phiDecay1 = photonMom.Phi();
1271 }
1272 if(photon->GetID()==indexPhoton2) {
1273 ptDecay2 = photonMom.Pt();
1274 pxDecay2 = photonMom.Px();
1275 pyDecay2 = photonMom.Py();
1276 phiDecay2 = photonMom.Phi();
1277 }
045396c8 1278 if(GetDebug() > 1)
1279 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
e09cf5ef 1280 } //photonAOD loop
1281 } //index of decay photons found
1282 if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);
1283 } //make decay-hadron correlation
1284
1285 TObjArray * refpi0 =0x0;
1286 Int_t nrefs = 0;
1287
1288 //Loop on stored AOD pi0
1289 Int_t naod = pi0list->GetEntriesFast();
045396c8 1290 if(GetDebug() > 0)
1291 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - aod branch entries %d\n", naod);
e09cf5ef 1292 for(Int_t iaod = 0; iaod < naod ; iaod++){
1293 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (pi0list->At(iaod));
1294
1295 Int_t evtIndex2 = 0 ;
1296 Int_t evtIndex3 = 0 ;
1297 if (GetMixedEvent()) {
1298 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
1299 evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
1300
045396c8 1301 if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
1302 evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
e09cf5ef 1303 continue ;
1304 }
1305
1306 //Int_t pdg = pi0->GetPdg();
1307 //if(pdg != AliCaloPID::kPi0) continue;
1308
1309 pt = pi0->Pt();
1310 px = pi0->Px();
1311 py = pi0->Py();
66e64043 1312
1313 if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1314
1315 //jumped out this event if near side associated particle pt larger than trigger
e09cf5ef 1316 if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;
045396c8 1317
e09cf5ef 1318 //Selection within angular range
1319 phi = pi0->Phi();
045396c8 1320 //Float_t deltaPhi = TMath::Abs(phiTrig-phi);
1321 //if( (deltaPhi < fDeltaPhiMinCut) || ( deltaPhi > fDeltaPhiMaxCut) ) continue ;
e09cf5ef 1322
1323 if(bFillHisto){
1324
045396c8 1325 deltaPhi = phiTrig-phi;
1326 if(deltaPhi<-TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1327 if(deltaPhi>3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
e09cf5ef 1328
1329 rat = pt/ptTrig ;
1330 phi = pi0->Phi() ;
1331 eta = pi0->Eta() ;
66e64043 1332
1333 rat = pt/ptTrig ;
1334 xE =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1335 //if(xE <0.)xE =-xE;
1336 cosi = -100;
1337 if(xE > 0 ) cosi = TMath::Log(1./xE);
e09cf5ef 1338
1339 if(fPi0Trigger){
1340 if(indexPhoton1!=-1 && indexPhoton2!=-1){
1341 if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
1342 if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
045396c8 1343 deltaPhiDecay1 = phiDecay1-phi;
1344 deltaPhiDecay2 = phiDecay2-phi;
1345 if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
1346 if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
1347 if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
1348 if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
1349 fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
1350 fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
1351 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
1352 if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
e09cf5ef 1353 fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1);
045396c8 1354 if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
e09cf5ef 1355 fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);
045396c8 1356 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
e09cf5ef 1357 }
1358 } //do decay-hadron correlation
1359
1360 fhEtaNeutral->Fill(pt,eta);
1361 fhPhiNeutral->Fill(pt,phi);
1362 fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
045396c8 1363 fhDeltaPhiNeutral->Fill(ptTrig,deltaPhi);
1364 fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi,etaTrig-eta);
e09cf5ef 1365
1366 //delta phi cut for correlation
045396c8 1367 if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) ) {
1368 fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
e09cf5ef 1369 fhPtImbalanceNeutral->Fill(ptTrig,rat);
1370 fhPtHbpNeutral->Fill(ptTrig,cosi);
1371 }
1372 else {
045396c8 1373 fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
e09cf5ef 1374 fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
1375 fhPtHbpUeNeutral->Fill(ptTrig,cosi);
1376 }
1377 //several UE calculation
1378 if(fMakeSeveralUE){
045396c8 1379 if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut)){
1380 fhDeltaPhiUeLeftNeutral->Fill(pt,deltaPhi);
e09cf5ef 1381 fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
1382 fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
1383 }
045396c8 1384 if((deltaPhi>fUeDeltaPhiMinCut) && (deltaPhi <fUeDeltaPhiMaxCut)){
1385 fhDeltaPhiUeRightNeutral->Fill(pt,deltaPhi);
e09cf5ef 1386 fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
1387 fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
1388 }
1389 } //several UE calculation
1390 }
1391 else{
1392 nrefs++;
1393 if(nrefs==1){
1394 refpi0 = new TObjArray(0);
1395 refpi0->SetName(GetAODObjArrayName()+"Pi0s");
1396 refpi0->SetOwner(kFALSE);
1397 }
1398 refpi0->Add(pi0);
1399 }//put references in trigger AOD
045396c8 1400
1401 if(GetDebug() > 2 )
1402 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
1403
1404 }//loop
1405
1406 return kTRUE;
1407}
1408
05d0d05d 1409//_________________________________________________________________________________________________________
045396c8 1410void AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
1411{
1412 // Charged Hadron Correlation Analysis with MC information
1413 if(GetDebug()>1)
1414 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Make trigger particle - charged hadron correlation in AOD MC level\n");
1415
1416 AliStack * stack = 0x0 ;
1417 TParticle * primary = 0x0 ;
1418 TClonesArray * mcparticles0 = 0x0 ;
1419 TClonesArray * mcparticles = 0x0 ;
1420 AliAODMCParticle * aodprimary = 0x0 ;
1421
1422 Double_t eprim = 0 ;
1423 Double_t ptprim = 0 ;
1424 Double_t phiprim = 0 ;
1425 Double_t etaprim = 0 ;
1426 Double_t pxprim = 0 ;
1427 Double_t pyprim = 0 ;
1428 Double_t pzprim = 0 ;
1429 Int_t nTracks = 0 ;
1430 Int_t iParticle = 0 ;
1431 Double_t charge = 0.;
1432
1433 Double_t mcrat =-100 ;
1434 Double_t mcxE =-100 ;
1435 Double_t mccosi =-100 ;
66e64043 1436 Double_t mcdeltaPhi = -100. ;
1437
045396c8 1438 //Track loop, select tracks with good pt, phi and fill AODs or histograms
1439 //Int_t currentIndex = -1 ;
1440 Double_t mcTrackPt = 0 ;
1441 Double_t mcTrackPhi = 0 ;
1442 Double_t mcTrackEta = 0 ;
1443 Double_t mcTrackPx = 0 ;
1444 Double_t mcTrackPy = 0 ;
1445 Double_t mcTrackPz = 0 ;
66e64043 1446
045396c8 1447 if(GetReader()->ReadStack()){
1448 nTracks = GetMCStack()->GetNtrack() ;
1449 }
1450 else nTracks = GetReader()->GetAODMCParticles()->GetEntriesFast() ;
1451 //Int_t trackIndex[nTracks];
1452
1453 Int_t label= aodParticle->GetLabel();
1454 if(label<0){
1455 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***: label %d \n", label);
1456 return;
1457 }
1458
1459 if(GetReader()->ReadStack()){
1460 stack = GetMCStack() ;
1461 if(!stack) {
1462 printf(" AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation- Stack not available, is the MC handler called? STOP\n");
1463 abort();
1464 }
1465
1466 nTracks=stack->GetNprimary();
1467 if(label >= stack->GetNtrack()) {
1468 if(GetDebug() > 2) printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1469 return ;
1470 }
1471 primary = stack->Particle(label);
1472 if(!primary){
1473 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no primary ***: label %d \n", label);
1474 return;
1475 }
1476
1477 eprim = primary->Energy();
1478 ptprim = primary->Pt();
1479 phiprim = primary->Phi();
1480 etaprim = primary->Eta();
1481 pxprim = primary->Px();
1482 pyprim = primary->Py();
1483 pzprim = primary->Pz();
1484
1485 if(primary){
e09cf5ef 1486
045396c8 1487 for (iParticle = 0 ; iParticle < nTracks ; iParticle++) {
1488 TParticle * particle = stack->Particle(iParticle);
1489 TLorentzVector momentum;
1490 //keep only final state particles
1491 if(particle->GetStatusCode()!=1) continue ;
1492 Int_t pdg = particle->GetPdgCode();
1493 charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1494 particle->Momentum(momentum);
1495
1496 //---------- Charged particles ----------------------
1497 if(charge != 0){
1498 //Particles in CTS acceptance
1499 Bool_t inCTS = GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
1500 if(TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
66e64043 1501 if(inCTS)
045396c8 1502 {
1503 mcTrackPt = particle->Pt();
1504 mcTrackPhi = particle->Phi();
1505 mcTrackEta = particle->Eta();
1506 mcTrackPx = particle->Px();
1507 mcTrackPy = particle->Py();
1508 mcTrackPz = particle->Pz();
66e64043 1509 if(mcTrackPhi < 0) mcTrackPhi+=TMath::TwoPi();
1510
045396c8 1511 //Select only hadrons in pt range
66e64043 1512 if(mcTrackPt < fMinAssocPt || mcTrackPt > fMaxAssocPt) continue ;
1513
045396c8 1514 //remove trigger itself for correlation when use charged triggers
c5693f62 1515 if(label==iParticle &&
1516 TMath::Abs(mcTrackPt -ptprim ) < 1e-6 &&
1517 TMath::Abs(mcTrackPhi-phiprim) < 1e-6 &&
1518 TMath::Abs(mcTrackEta-etaprim) < 1e-6)
66e64043 1519 continue ;
1520
045396c8 1521 //jumped out this event if near side associated partile pt larger than trigger
66e64043 1522 if( mcTrackPt > ptprim && TMath::Abs(mcTrackPhi-phiprim) < TMath::PiOver2())
045396c8 1523 return ;
1524
66e64043 1525 mcdeltaPhi= phiprim-mcTrackPhi;
1526 if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
1527 if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
1528
1529 mcrat = mcTrackPt/ptprim ;
1530 mcdeltaPhi= phiprim-mcTrackPhi;
1531 mcxE =-mcTrackPt/ptprim*TMath::Cos(mcdeltaPhi);// -(mcTrackPx*pxprim+mcTrackPy*pyprim)/(ptprim*ptprim);
1532 mccosi =-100 ;
1533 if(mcxE > 0 ) mccosi = TMath::Log(1./mcxE);
1534
045396c8 1535 // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
1536 // printf("phi = %f \n", phi);
1537
1538 //Selection within angular range
045396c8 1539 if( mcdeltaPhi< -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
1540 if( mcdeltaPhi>3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
1541 Double_t mcpout = mcTrackPt*TMath::Sin(mcdeltaPhi) ;
1542 if(GetDebug()>0 )
1543 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
1544 mcTrackPt,mcTrackPhi, phiprim,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());
1545 // Fill Histograms
1546 fhMCEtaCharged->Fill(mcTrackPt,mcTrackEta);
1547 fhMCPhiCharged->Fill(mcTrackPt,mcTrackPhi);
1548 fhMCDeltaEtaCharged->Fill(ptprim,etaprim-mcTrackEta);
1549 fhMCDeltaPhiCharged->Fill(ptprim,mcdeltaPhi);
1550 fhMCPtAssocDeltaPhi->Fill(mcTrackPt, mcdeltaPhi);
045396c8 1551 fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,etaprim-mcTrackEta);
1552
1553 //delta phi cut for correlation
1554 if( (mcdeltaPhi > fDeltaPhiMinCut) && ( mcdeltaPhi < fDeltaPhiMaxCut) ) {
1555 fhMCDeltaPhiChargedPt->Fill(mcTrackPt,mcdeltaPhi);
1556 fhMCPtImbalanceCharged->Fill(ptprim,mcxE);
1557 fhMCPtHbpCharged->Fill(ptprim,mccosi);
1558 fhMCPtTrigPout->Fill(ptprim, mcpout) ;
1559 }//delta phi cut for correlation
1560 } //tracks after cuts
1561 }//Charged
1562 } //track loop
1563 } //when the leading particles could trace back to MC
1564 } //ESD MC
1565 else if(GetReader()->ReadAODMCParticles()){
1566 //Get the list of MC particles
1567 mcparticles0 = GetReader()->GetAODMCParticles(0);
1568 if(!mcparticles0) return;
1569 if(label >=mcparticles0->GetEntriesFast()){
1570 if(GetDebug() > 2)
1571 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***: label %d, n tracks %d \n", label,mcparticles0->GetEntriesFast());
1572 return;
1573 }
1574 //Get the particle
1575 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1576 if(!aodprimary) {
1577 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no AOD primary ***: label %d \n", label);
1578 return;
1579 }
1580
1581 ptprim = aodprimary->Pt();
66e64043 1582 phiprim = aodprimary->Phi();
1583 etaprim = aodprimary->Eta();
1584 pxprim = aodprimary->Px();
1585 pyprim = aodprimary->Py();
1586 pzprim = aodprimary->Pz();
045396c8 1587 if(aodprimary){
1588 mcparticles= GetReader()->GetAODMCParticles();
1589 for (Int_t i=0; i<nTracks;i++) {
1590 AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(i);
1591 if (!part->IsPhysicalPrimary()) continue;
1592 Int_t pdg = part->GetPdgCode();
1593 charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1594 TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
1595 if(charge != 0){
1596 if(part->Pt()> GetReader()->GetCTSPtMin()){
1597 //Particles in CTS acceptance
1598 Bool_t inCTS = GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
1599 Int_t indexmother=part->GetMother();
1600 if(indexmother>-1)
1601 {
1602 Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
1603 if(TMath::Abs(pdg) == 11 && mPdg == 22) continue;
1604 }
66e64043 1605 if(inCTS)
045396c8 1606 {
66e64043 1607 mcTrackPt =part->Pt();
1608 mcTrackPhi =part->Phi();
1609 mcTrackEta =part->Eta();
1610 mcTrackPx =part->Px();
1611 mcTrackPy =part->Py();
1612 mcTrackPz =part->Pz();
1613 if(mcTrackPhi < 0) mcTrackPhi+=TMath::TwoPi();
1614
045396c8 1615 //Select only hadrons in pt range
1616 if(mcTrackPt < GetMinPt() || mcTrackPt > GetMaxPt()) continue ;
c5693f62 1617
045396c8 1618 //remove trigger itself for correlation when use charged triggers
c5693f62 1619 if(label==i &&
1620 TMath::Abs(mcTrackPt -ptprim ) < 1e-6 &&
1621 TMath::Abs(mcTrackPhi-phiprim) < 1e-6 &&
1622 TMath::Abs(mcTrackEta-etaprim) < 1e-6) continue ;
1623
045396c8 1624 //jumped out this event if near side associated partile pt larger than trigger
1625 if( mcTrackPt> ptprim && TMath::Abs(mcTrackPhi-phiprim)<TMath::PiOver2())
1626 return ;
1627
66e64043 1628 mcdeltaPhi= phiprim-mcTrackPhi;
1629 if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
1630 if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
1631
1632 mcrat = mcTrackPt/ptprim ;
1633 mcxE =-mcTrackPt/ptprim*TMath::Cos(mcdeltaPhi);// -(mcTrackPx*pxprim+mcTrackPy*pyprim)/(ptprim*ptprim);
1634 mccosi =-100 ;
1635 if(mcxE > 0 ) mccosi = TMath::Log(1./mcxE);
045396c8 1636
1637 //Selection within angular range
045396c8 1638 if( mcdeltaPhi< -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
1639 if( mcdeltaPhi>3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
1640 Double_t mcpout = mcTrackPt*TMath::Sin(mcdeltaPhi) ;
1641 if(GetDebug()>0)
1642 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
1643 mcTrackPt,mcTrackPhi, phiprim,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());
fe871b21 1644
045396c8 1645 // Fill Histograms
fe871b21 1646
66e64043 1647 fhMCEtaCharged ->Fill(mcTrackPt,mcTrackEta);
1648 fhMCPhiCharged ->Fill(mcTrackPt,mcTrackPhi);
045396c8 1649 fhMCDeltaEtaCharged->Fill(ptprim,etaprim-mcTrackEta);
1650 fhMCDeltaPhiCharged->Fill(ptprim,mcdeltaPhi);
1651 fhMCPtAssocDeltaPhi->Fill(mcTrackPt, mcdeltaPhi);
045396c8 1652 fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,etaprim-mcTrackEta);
1653
1654 //delta phi cut for correlation
1655 if( (mcdeltaPhi > fDeltaPhiMinCut) && ( mcdeltaPhi < fDeltaPhiMaxCut) ) {
66e64043 1656 fhMCDeltaPhiChargedPt ->Fill(mcTrackPt,mcdeltaPhi);
045396c8 1657 fhMCPtImbalanceCharged->Fill(ptprim,mcxE);
66e64043 1658 fhMCPtHbpCharged ->Fill(ptprim,mccosi);
1659 fhMCPtTrigPout ->Fill(ptprim, mcpout) ;
045396c8 1660 }//delta phi cut for correlation
1661
1662 } //tracks after cuts
1663
1664 } //with minimum pt cut
1665 } //only charged particles
1666 } //MC particle loop
1667 } //when the leading particles could trace back to MC
1668 }// AOD MC
e09cf5ef 1669}
045396c8 1670
05d0d05d 1671//_____________________________________________________________________
1672void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
1673{
1674
1675 //Print some relevant parameters set for the analysis
1676 if(! opt)
1677 return;
1678
1679 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 1680 AliAnaCaloTrackCorrBaseClass::Print(" ");
66e64043 1681 printf("Pt trigger > %3.2f\n", fMinTriggerPt) ;
1682 printf("Pt associated hadron < %3.2f\n", fMaxAssocPt) ;
1683 printf("Pt associated hadron > %3.2f\n", fMinAssocPt) ;
1684 printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ;
1685 printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;
1686 printf("Phi trigger particle-UeHadron < %3.2f\n", fUeDeltaPhiMaxCut) ;
1687 printf("Phi trigger particle-UeHadron > %3.2f\n", fUeDeltaPhiMinCut) ;
1688 printf("Isolated Trigger? %d\n" , fSelectIsolated) ;
1689 printf("Several UE? %d\n" , fMakeSeveralUE) ;
1690 printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
05d0d05d 1691 printf("Do Decay-hadron correlation ? %d\n", fPi0Trigger) ;
1692 printf("Select absolute leading for cluster triggers ? %d\n", fMakeAbsoluteLeading) ;
1693 printf("Trigger pt bins %d\n", fNAssocPtBins) ;
1694 for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
1695 printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
1696 }
1697
1698}
1699
1700//____________________________________________________________
1701void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
045396c8 1702{
1703 // Set number of bins
e09cf5ef 1704
05d0d05d 1705 fNAssocPtBins = n ;
1706
1707
045396c8 1708 if(n < 10 && n > 0)
1709 {
05d0d05d 1710 fNAssocPtBins = n ;
045396c8 1711 }
1712 else
1713 {
1714 printf("n = larger than 9 or too small, set to 9 \n");
05d0d05d 1715 fNAssocPtBins = 9;
045396c8 1716 }
1717}
e09cf5ef 1718
05d0d05d 1719//______________________________________________________________________________
1720void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
045396c8 1721{
1722 // Set the list of limits for the trigger pt bins
1723
05d0d05d 1724 if(ibin <= fNAssocPtBins || ibin >= 0)
045396c8 1725 {
05d0d05d 1726 fAssocPtBinLimit[ibin] = pt ;
045396c8 1727 }
1728 else {
05d0d05d 1729 printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ;
045396c8 1730
1731 }
1732}
1733