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