6afc9bca363d010a920c7f3c11dbbed3414f0d3f
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / 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
16 //_________________________________________________________________________
17 // Class for the analysis of particle - hadron correlations
18 // Particle (for example direct gamma) must be found
19 // in a previous analysis
20 //
21 //-- Author: Gustavo Conesa (LNF-INFN) (LPSC-IN2P3-CNRS)
22 //           Yaxian Mao (LPSC-IN2P3-CNRS) and (CNWU) first usable implementation.
23 //           Xiangrong Zhu (CNWU), implementtion of own mixing.
24 //////////////////////////////////////////////////////////////////////////////
25
26
27 // --- ROOT system ---
28 //#include "TClonesArray.h"
29 #include <TClass.h>
30 #include <TMath.h>
31 #include <TH2F.h>
32 #include <TDatabasePDG.h>
33
34 //---- ANALYSIS system ----
35 #include "AliNeutralMesonSelection.h" 
36 #include "AliAnaParticleHadronCorrelation.h" 
37 #include "AliCaloTrackReader.h"
38 #include "AliAODPWG4ParticleCorrelation.h"
39 #include "AliFiducialCut.h"
40 #include "AliVTrack.h"
41 #include "AliVCluster.h"
42 #include "AliMCAnalysisUtils.h"
43 #include "TParticle.h"
44 #include "AliStack.h"
45 #include "AliAODMCParticle.h"
46 #include "AliMixedEvent.h"
47 #include "AliAnalysisManager.h"
48 #include "AliInputEventHandler.h"
49 #include "AliEventplane.h"
50
51 ClassImp(AliAnaParticleHadronCorrelation)
52
53
54 //___________________________________________________________________
55   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(): 
56     AliAnaCaloTrackCorrBaseClass(),
57     fFillAODWithReferences(0),      fCheckLeadingWithNeutralClusters(0),
58     fMaxAssocPt(1000.),             fMinAssocPt(0.),
59     fDeltaPhiMaxCut(0.),            fDeltaPhiMinCut(0.),   
60     fSelectIsolated(0),             fMakeSeveralUE(0),              
61     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.), 
62     fPi0AODBranchName(""),          fNeutralCorr(0),       
63     fPi0Trigger(0),                 fDecayTrigger(0),
64     fNDecayBits(0),                 fDecayBits(),
65     fMakeAbsoluteLeading(0),        fMakeNearSideLeading(0),      
66     fLeadingTriggerIndex(-1),       fHMPIDCorrelation(0),  fFillBradHisto(0),
67     fNAssocPtBins(0),               fAssocPtBinLimit(),
68     fCorrelVzBin(0),
69     fListMixTrackEvents(),          fListMixCaloEvents(),
70     fUseMixStoredInReader(0),       fFillNeutralEventMixPool(0),
71     fM02MaxCut(0),                  fM02MinCut(0),  
72     fSelectLeadingHadronAngle(0),   fFillLeadHadOppositeHisto(0),
73     fMinLeadHadPhi(0),              fMaxLeadHadPhi(0),
74     fMinLeadHadPt(0),               fMaxLeadHadPt(0),
75     fFillEtaGapsHisto(1),           fFillMomImbalancePtAssocBinsHisto(0),
76     fMCGenTypeMin(0),               fMCGenTypeMax(0),
77     fTrackVector(),                 fMomentum(),
78     fDecayMom1(),                   fDecayMom2(),
79     //Histograms
80     fhPtTriggerInput(0),            fhPtTriggerSSCut(0),
81     fhPtTriggerIsoCut(0),           fhPtTriggerFidCut(0),
82     fhPtTrigger(0),                 fhPtTriggerVtxBC0(0),
83     fhPtTriggerVzBin(0),            fhPtTriggerBin(0),                 
84     fhPhiTrigger(0),                fhEtaTrigger(0),   
85     fhPtTriggerMC(),
86     fhPtDecayTrigger(),             fhPtDecayTriggerMC(),
87     fhPtTriggerCentrality(0),       fhPtTriggerEventPlane(0), 
88     fhTriggerEventPlaneCentrality(0),
89     fhPtTriggerMixed(0),            fhPtTriggerMixedVzBin(0), fhPtTriggerMixedBin(0),              
90     fhPhiTriggerMixed(0),           fhEtaTriggerMixed(0),
91     fhPtLeadingOppositeHadron(0),   fhPtDiffPhiLeadingOppositeHadron(0), fhPtDiffEtaLeadingOppositeHadron(0),
92     fhPtNoLeadingOppositeHadron(0), fhEtaPhiNoLeadingOppositeHadron(0),
93     fhDeltaPhiDeltaEtaCharged(0),
94     fhPhiCharged(0),                fhEtaCharged(0), 
95     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0), 
96     fhDeltaPhiChargedPt(0),         fhDeltaPhiUeChargedPt(0), 
97     fhUePart(0),
98     fhXECharged(0),                 fhXECharged_Cone2(0),      fhXEUeCharged(0),
99     fhXEPosCharged(0),              fhXENegCharged(0),
100     fhPtHbpXECharged(0),            fhPtHbpXECharged_Cone2(0), fhPtHbpXEUeCharged(0),
101     fhZTCharged(0),                 fhZTUeCharged(0),
102     fhZTPosCharged(0),              fhZTNegCharged(0),
103     fhPtHbpZTCharged(0),            fhPtHbpZTUeCharged(0),
104     fhXEChargedMC(),                fhDeltaPhiChargedMC(),  
105     fhDeltaPhiDeltaEtaChargedPtA3GeV(0),
106     fhDeltaPhiChargedPtA3GeV(0),    fhDeltaEtaChargedPtA3GeV(0),
107     //Pile-Up
108     fhDeltaPhiChargedPileUp(),      fhDeltaEtaChargedPileUp(),
109     fhDeltaPhiChargedPtA3GeVPileUp(), fhDeltaEtaChargedPtA3GeVPileUp(),
110     fhXEChargedPileUp(),            fhXEUeChargedPileUp(),
111     fhZTChargedPileUp(),            fhZTUeChargedPileUp(), 
112     fhPtTrigChargedPileUp(),
113     fhDeltaPhiChargedOtherBC(),     fhDeltaPhiChargedPtA3GeVOtherBC(),
114     fhXEChargedOtherBC(),           fhXEUeChargedOtherBC(),
115     fhZTChargedOtherBC(),           fhZTUeChargedOtherBC(),
116     fhPtTrigChargedOtherBC(),
117     fhDeltaPhiChargedBC0(),         fhDeltaPhiChargedPtA3GeVBC0(),
118     fhXEChargedBC0(),               fhXEUeChargedBC0(),
119     fhZTChargedBC0(),               fhZTUeChargedBC0(),
120     fhPtTrigChargedBC0(),
121     fhDeltaPhiChargedVtxBC0(),      fhDeltaPhiChargedPtA3GeVVtxBC0(),
122     fhXEChargedVtxBC0(),            fhXEUeChargedVtxBC0(),
123     fhZTChargedVtxBC0(),            fhZTUeChargedVtxBC0(),
124     fhPtTrigChargedVtxBC0(),
125     fhDeltaPhiUeLeftCharged(0),
126     fhDeltaPhiUeLeftUpCharged(0),   fhDeltaPhiUeRightUpCharged(0),
127     fhDeltaPhiUeLeftDownCharged(0), fhDeltaPhiUeRightDownCharged(0),
128     fhXEUeLeftCharged(0),
129     fhXEUeLeftUpCharged(0),         fhXEUeRightUpCharged(0),
130     fhXEUeLeftDownCharged(0),       fhXEUeRightDownCharged(0),
131     fhPtHbpXEUeLeftCharged(0),      fhZTUeLeftCharged(0),
132     fhPtHbpZTUeLeftCharged(0),
133     fhPtTrigPout(0),                fhPtTrigCharged(0),
134     fhDeltaPhiChargedMult(0x0),     fhDeltaEtaChargedMult(0x0),
135     fhXEMult(0x0),                  fhXEUeMult(0x0),
136     fhZTMult(0x0),                  fhZTUeMult(0x0),
137     fhAssocPtBkg(0),                fhDeltaPhiDeltaEtaAssocPtBin(0),
138     fhDeltaPhiAssocPtBin(0),        
139     fhDeltaPhiAssocPtBinDEta08(0),  fhDeltaPhiAssocPtBinDEta0(0),
140     fhDeltaPhiAssocPtBinHMPID(0),   fhDeltaPhiAssocPtBinHMPIDAcc(0),
141     fhDeltaPhiBradAssocPtBin(0),    fhDeltaPhiBrad(0),
142     fhXEAssocPtBin(0),              fhZTAssocPtBin(0),
143     fhXEVZ(0),                      fhZTVZ(0),
144     fhDeltaPhiDeltaEtaNeutral(0),
145     fhPhiNeutral(0),                fhEtaNeutral(0), 
146     fhDeltaPhiNeutral(0),           fhDeltaEtaNeutral(0),
147     fhDeltaPhiNeutralPt(0),         fhDeltaPhiUeNeutralPt(0), 
148     fhXENeutral(0),                 fhXEUeNeutral(0),
149     fhPtHbpXENeutral(0),            fhPtHbpXEUeNeutral(0),
150     fhZTNeutral(0),                 fhZTUeNeutral(0),
151     fhPtHbpZTNeutral(0),            fhPtHbpZTUeNeutral(0),
152     fhDeltaPhiUeLeftNeutral(0),     fhXEUeLeftNeutral(0),
153     fhPtHbpXEUeLeftNeutral(0),      fhZTUeLeftNeutral(0),
154     fhPtHbpZTUeLeftNeutral(0),      fhPtPi0DecayRatio(0),
155     fhDeltaPhiPi0DecayCharged(0),   fhXEPi0DecayCharged(0),        fhZTPi0DecayCharged(0),
156     fhDeltaPhiPi0DecayNeutral(0),   fhXEPi0DecayNeutral(0),        fhZTPi0DecayNeutral(0),
157     fhDeltaPhiDecayCharged(),       fhXEDecayCharged(),            fhZTDecayCharged(),
158     fhDeltaPhiDecayChargedAssocPtBin(),
159     fhMCPtTrigger(),                fhMCPhiTrigger(),              fhMCEtaTrigger(),
160     fhMCPtTriggerNotLeading(),      fhMCPhiTriggerNotLeading(),    fhMCEtaTriggerNotLeading(),
161     fhMCEtaCharged(),               fhMCPhiCharged(),
162     fhMCDeltaEtaCharged(),          fhMCDeltaPhiCharged(),
163     fhMCDeltaPhiDeltaEtaCharged(),  fhMCDeltaPhiChargedPt(),
164     fhMCPtXECharged(),              fhMCPtXEUeCharged(),
165     fhMCPtXEUeLeftCharged(),
166     fhMCPtHbpXECharged(),           fhMCPtHbpXEUeCharged(),
167     fhMCPtHbpXEUeLeftCharged(),
168     fhMCUePart(),
169     fhMCPtZTCharged(),              fhMCPtZTUeCharged(),
170     fhMCPtZTUeLeftCharged(),
171     fhMCPtHbpZTCharged(),           fhMCPtHbpZTUeCharged(),
172     fhMCPtHbpZTUeLeftCharged(),
173     fhMCPtTrigPout(),               fhMCPtAssocDeltaPhi(),
174     //Mixing
175     fhNEventsTrigger(0),            fhNtracksMB(0),                 fhNclustersMB(0),
176     fhMixDeltaPhiCharged(0),        fhMixDeltaPhiDeltaEtaCharged(0),
177     fhMixXECharged(0),              fhMixXEUeCharged(0),            fhMixHbpXECharged(0),
178     fhMixDeltaPhiChargedAssocPtBin(),
179     fhMixDeltaPhiChargedAssocPtBinDEta08(),
180     fhMixDeltaPhiChargedAssocPtBinDEta0(),
181     fhMixDeltaPhiDeltaEtaChargedAssocPtBin(),
182     fhEventBin(0),                  fhEventMixBin(0),               fhEventMBBin(0)
183
184   //Default Ctor 
185     
186   //Initialize parameters
187   InitParameters();
188   
189   for(Int_t i = 0; i < fgkNmcTypes; i++)
190   { 
191     fhPtTriggerMC[i] = 0;
192     fhXEChargedMC[i] = 0;
193     fhDeltaPhiChargedMC[i] = 0;
194     for(Int_t ib = 0; ib < 4; ib++)  fhPtDecayTriggerMC[ib][i] = 0;
195   }
196
197   for(Int_t ib = 0; ib < 4; ib++)  fhPtDecayTrigger[ib] = 0;
198
199   for(Int_t i = 0; i < 7; i++)
200   {
201     fhPtTriggerPileUp             [i] = 0 ;
202     fhDeltaPhiChargedPileUp       [i] = 0 ; fhDeltaEtaChargedPileUp       [i] = 0 ;
203     fhXEChargedPileUp             [i] = 0 ; fhXEUeChargedPileUp           [i] = 0 ;
204     fhZTChargedPileUp             [i] = 0 ; fhZTUeChargedPileUp           [i] = 0 ;
205     fhPtTrigChargedPileUp         [i] = 0 ;
206     fhDeltaPhiChargedPtA3GeVPileUp[i] = 0 ; fhDeltaEtaChargedPtA3GeVPileUp[i] = 0 ;
207   }
208   
209 }
210
211 //_________________________________________________________________
212 AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation()
213 {
214   // Remove event containers
215   
216   if(DoOwnMix())
217   {
218     if(fListMixTrackEvents)
219     {      
220       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
221       {
222         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
223         {
224           for(Int_t irp=0; irp<GetNRPBin(); irp++)
225           {
226             Int_t bin = GetEventMixBin(ic, iz, irp);
227             fListMixTrackEvents[bin]->Delete() ;
228             delete fListMixTrackEvents[bin] ;
229           }
230         }
231       }
232     }
233
234     delete[] fListMixTrackEvents;
235
236     if(fListMixCaloEvents)
237     {      
238       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
239       {
240         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
241         {
242           for(Int_t irp=0; irp<GetNRPBin(); irp++)
243           {
244             Int_t bin = GetEventMixBin(ic, iz, irp);
245             fListMixCaloEvents[bin]->Delete() ;
246             delete fListMixCaloEvents[bin] ;
247           }
248         }
249       }
250     }
251   
252     delete[] fListMixCaloEvents;
253     
254   }
255 }
256
257 //__________________________________________________________________________________________________________________________________________
258 void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Float_t ptAssoc,  Float_t ptTrig,      Int_t   bin,
259                                                                               Float_t phiAssoc, Float_t phiTrig,     Float_t deltaPhi,
260                                                                               Float_t etaAssoc, Float_t etaTrig,  
261                                                                               Int_t   decayTag, Float_t hmpidSignal, Int_t  outTOF,
262                                                                               Int_t   cen,      Int_t   mcTag)
263 {
264   // Fill angular correlation related histograms
265   
266   Float_t deltaEta    = etaTrig-etaAssoc;
267   Float_t deltaPhiOrg = phiTrig-phiAssoc;
268   
269   fhEtaCharged       ->Fill(ptAssoc,etaAssoc);
270   fhPhiCharged       ->Fill(ptAssoc,phiAssoc);
271   fhDeltaEtaCharged  ->Fill(ptTrig ,deltaEta);
272   fhDeltaPhiCharged  ->Fill(ptTrig ,deltaPhi);
273   fhDeltaPhiChargedPt->Fill(ptAssoc, deltaPhi);
274   fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
275   
276   if(ptAssoc > 3 )
277   {
278     fhDeltaEtaChargedPtA3GeV        ->Fill(ptTrig  ,deltaEta);
279     fhDeltaPhiChargedPtA3GeV        ->Fill(ptTrig  ,deltaPhi);
280     fhDeltaPhiDeltaEtaChargedPtA3GeV->Fill(deltaPhi, deltaEta);    
281   }   
282   
283   // Pile up studies
284   
285   if(IsPileUpAnalysisOn())
286   {
287     if     (outTOF==1)
288     {
289       fhDeltaPhiChargedOtherBC->Fill(ptTrig ,deltaPhi) ;
290       if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVOtherBC->Fill(ptTrig ,deltaPhi) ;
291     }
292     else if(outTOF==0)
293     {
294         fhDeltaPhiChargedBC0->Fill(ptTrig ,deltaPhi) ;
295         if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVBC0->Fill(ptTrig ,deltaPhi) ;
296     }
297     
298     Int_t vtxBC = GetReader()->GetVertexBC();
299     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
300     {
301       fhDeltaPhiChargedVtxBC0->Fill(ptTrig ,deltaPhi) ;
302       if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVVtxBC0->Fill(ptTrig ,deltaPhi) ;
303     }
304     
305     if(GetReader()->IsPileUpFromSPD())               { fhDeltaEtaChargedPileUp[0]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[0]->Fill(ptTrig ,deltaPhi) ; }
306     if(GetReader()->IsPileUpFromEMCal())             { fhDeltaEtaChargedPileUp[1]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[1]->Fill(ptTrig ,deltaPhi) ; }
307     if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhDeltaEtaChargedPileUp[2]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[2]->Fill(ptTrig ,deltaPhi) ; }
308     if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhDeltaEtaChargedPileUp[3]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[3]->Fill(ptTrig ,deltaPhi) ; }
309     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhDeltaEtaChargedPileUp[4]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[4]->Fill(ptTrig ,deltaPhi) ; }
310     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhDeltaEtaChargedPileUp[5]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[5]->Fill(ptTrig ,deltaPhi) ; }
311     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhDeltaEtaChargedPileUp[6]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[6]->Fill(ptTrig ,deltaPhi) ; }
312     
313     if(ptAssoc > 3 )
314     {
315       if(GetReader()->IsPileUpFromSPD())               { fhDeltaEtaChargedPtA3GeVPileUp[0]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[0]->Fill(ptTrig ,deltaPhi) ; }
316       if(GetReader()->IsPileUpFromEMCal())             { fhDeltaEtaChargedPtA3GeVPileUp[1]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[1]->Fill(ptTrig ,deltaPhi) ; }
317       if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhDeltaEtaChargedPtA3GeVPileUp[2]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[2]->Fill(ptTrig ,deltaPhi) ; }
318       if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhDeltaEtaChargedPtA3GeVPileUp[3]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[3]->Fill(ptTrig ,deltaPhi) ; }
319       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhDeltaEtaChargedPtA3GeVPileUp[4]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[4]->Fill(ptTrig ,deltaPhi) ; }
320       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhDeltaEtaChargedPtA3GeVPileUp[5]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[5]->Fill(ptTrig ,deltaPhi) ; }
321       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhDeltaEtaChargedPtA3GeVPileUp[6]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[6]->Fill(ptTrig ,deltaPhi) ; }
322     }
323   }
324   
325   if(IsDataMC())
326   {
327     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
328     fhDeltaPhiChargedMC[mcIndex]->Fill(ptTrig , deltaPhi);
329     if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==2 )
330       fhDeltaPhiChargedMC[7]->Fill(ptTrig , deltaPhi);
331   }
332   
333   if(fDecayTrigger && decayTag > 0)
334   {
335     for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
336     {
337       if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit])) fhDeltaPhiDecayCharged[ibit]->Fill(ptTrig,deltaPhi);
338     }
339   }
340   
341   Double_t  dphiBrad = -100;
342   if(fFillBradHisto)
343   {
344     dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
345     if( TMath::Abs(dphiBrad) > 0.325 && TMath::Abs(dphiBrad) < 0.475 )  //Hardcoded values, BAD, FIXME
346     {
347       fhAssocPtBkg->Fill(ptTrig, ptAssoc);
348     }
349     
350     if( dphiBrad < -1./3 ) dphiBrad += 2;
351     fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
352   }
353   
354   // Fill histograms in bins of associated particle pT
355   if(bin>=0)
356   {
357       fhDeltaPhiDeltaEtaAssocPtBin    [bin]->Fill(deltaPhi,deltaEta);
358
359       fhDeltaPhiAssocPtBin            [bin]->Fill(ptTrig, deltaPhi);
360     
361     if(fFillEtaGapsHisto)
362     {
363       if(TMath::Abs(deltaEta)> 0.8)
364         fhDeltaPhiAssocPtBinDEta08      [bin]->Fill(ptTrig, deltaPhi);
365       
366       if(TMath::Abs(deltaEta)< 0.01)
367         fhDeltaPhiAssocPtBinDEta0       [bin]->Fill(ptTrig, deltaPhi);
368     }
369     
370     if (fFillBradHisto)
371       fhDeltaPhiBradAssocPtBin        [bin]->Fill(ptTrig, dphiBrad);
372     
373     if(fDecayTrigger && decayTag > 0 && fNDecayBits > 0 &&
374        GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
375       fhDeltaPhiDecayChargedAssocPtBin[bin]->Fill(ptTrig, deltaPhi);      
376     
377     if(fHMPIDCorrelation)
378     {
379       if( hmpidSignal > 0 )
380       {
381         //printf("Track pt %f with HMPID signal %f \n",pt,hmpidSignal);
382         fhDeltaPhiAssocPtBinHMPID[bin]->Fill(ptTrig, deltaPhi);        
383       }
384       
385       if(phiAssoc > 5*TMath::DegToRad() && phiAssoc < 20*TMath::DegToRad())
386       {
387         //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
388         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->Fill(ptTrig, deltaPhi);        
389       }
390     }
391   }
392   
393   //fill different multiplicity/centrality histogram
394   if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
395   {
396     fhDeltaPhiChargedMult[cen]->Fill(ptTrig,deltaPhi);
397     fhDeltaEtaChargedMult[cen]->Fill(ptTrig,deltaEta);
398   }  
399 }
400
401 //___________________________________________________________________________________________________________________________________
402 Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float_t mcAssocPt, Float_t mcAssocPhi, Float_t mcAssocEta,
403                                                                            Float_t mcTrigPt,  Float_t mcTrigPhi,  Float_t mcTrigEta,
404                                                                            Int_t histoIndex, Bool_t lostDecayPair)
405 {
406   // Fill MC histograms independently of AOD or ESD
407   
408   Bool_t lead = kTRUE;
409   
410   // In case we requested the trigger to be a leading particle,
411   // check if this is true at the MC level.
412   // Not sure if it is correct to skip or not skip this.
413   // Absolute leading?
414   if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt )     lead = kFALSE; // skip event
415   
416   // Skip this event if near side associated particle pt larger than trigger
417   if( mcAssocPhi < 0 ) mcAssocPhi+=TMath::TwoPi();
418   
419   Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi;
420   if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
421   if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
422
423   if( fMakeNearSideLeading)
424   {
425     if( mcAssocPt > mcTrigPt && mcdeltaPhi < TMath::PiOver2() ) lead = kFALSE; // skip event
426   }
427   
428   //
429   // Select only hadrons in pt range
430   if ( mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt ) return lead ; // exclude but continue
431   if ( mcAssocPt < GetReader()->GetCTSPtMin())              return lead ;
432
433   
434   //
435   // Remove trigger itself for correlation when use charged triggers
436   if(TMath::Abs(mcAssocPt -mcTrigPt ) < 1e-6 && 
437      mcdeltaPhi                       < 1e-6 &&
438      TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6)            return lead ; // exclude but continue
439   
440   Float_t mcxE    =-mcAssocPt/mcTrigPt*TMath::Cos(mcdeltaPhi);// -(mcAssocPx*pxprim+mcAssocPy*pyprim)/(mcTrigPt*mcTrigPt);  
441   Float_t mchbpXE =-100 ;
442   if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
443   
444   Float_t mczT    = mcAssocPt/mcTrigPt ;
445   Float_t mchbpZT =-100 ;
446   if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
447   
448   Double_t mcpout = mcAssocPt*TMath::Sin(mcdeltaPhi) ; 
449   
450   if(GetDebug() > 0 )
451   {
452     AliInfo(Form("Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f \n",
453                  mcAssocPt,mcAssocPhi, mcTrigPhi,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut));
454   }
455   
456   // Fill Histograms
457   fhMCEtaCharged     [histoIndex]->Fill(mcAssocPt, mcAssocEta);
458   fhMCPhiCharged     [histoIndex]->Fill(mcAssocPt, mcAssocPhi);
459   fhMCDeltaEtaCharged[histoIndex]->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
460   fhMCDeltaPhiCharged[histoIndex]->Fill(mcTrigPt , mcdeltaPhi);
461   fhMCPtAssocDeltaPhi[histoIndex]->Fill(mcAssocPt, mcdeltaPhi);
462   
463   fhMCDeltaPhiDeltaEtaCharged[histoIndex]->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
464   
465   //delta phi cut for correlation
466   if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) ) 
467   {
468     fhMCDeltaPhiChargedPt[histoIndex]->Fill(mcAssocPt,mcdeltaPhi);
469     fhMCPtXECharged      [histoIndex]->Fill(mcTrigPt, mcxE);
470     fhMCPtHbpXECharged   [histoIndex]->Fill(mcTrigPt, mchbpXE);
471     fhMCPtZTCharged      [histoIndex]->Fill(mcTrigPt, mczT);
472     fhMCPtHbpZTCharged   [histoIndex]->Fill(mcTrigPt, mchbpZT);
473     fhMCPtTrigPout       [histoIndex]->Fill(mcTrigPt, mcpout) ;
474   }
475
476   if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
477   {
478     fhMCEtaCharged     [7]->Fill(mcAssocPt, mcAssocEta);
479     fhMCPhiCharged     [7]->Fill(mcAssocPt, mcAssocPhi);
480     fhMCDeltaEtaCharged[7]->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
481     fhMCDeltaPhiCharged[7]->Fill(mcTrigPt , mcdeltaPhi);
482     fhMCPtAssocDeltaPhi[7]->Fill(mcAssocPt, mcdeltaPhi);
483     
484     fhMCDeltaPhiDeltaEtaCharged[7]->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
485     
486     //delta phi cut for correlation
487     if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
488     {
489       fhMCDeltaPhiChargedPt[7]->Fill(mcAssocPt,mcdeltaPhi);
490       fhMCPtXECharged      [7]->Fill(mcTrigPt, mcxE);
491       fhMCPtHbpXECharged   [7]->Fill(mcTrigPt, mchbpXE);
492       fhMCPtZTCharged      [7]->Fill(mcTrigPt, mczT);
493       fhMCPtHbpZTCharged   [7]->Fill(mcTrigPt, mchbpZT);
494       fhMCPtTrigPout       [7]->Fill(mcTrigPt, mcpout) ;
495     }
496   }
497   
498   // Underlying event
499   
500   // Right
501   if ( (mcdeltaPhi > fUeDeltaPhiMinCut) && (mcdeltaPhi < fUeDeltaPhiMaxCut) )
502   {
503     //Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
504     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
505     Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
506     Double_t mcUezT =   mcAssocPt/mcTrigPt;
507     
508     if(mcUexE < 0.)
509       printf("FillChargedMCCorrelationHistograms(): Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
510              mcUexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
511
512     fhMCPtXEUeCharged[histoIndex]->Fill(mcTrigPt,mcUexE);
513     if(mcUexE > 0) fhMCPtHbpXEUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
514     
515     fhMCPtZTUeCharged[histoIndex]->Fill(mcTrigPt,mcUezT);
516     if(mcUezT > 0) fhMCPtHbpZTUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
517     
518     fhMCUePart[histoIndex]->Fill(mcTrigPt);
519     
520     if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
521     {
522       fhMCPtXEUeCharged[7]->Fill(mcTrigPt,mcUexE);
523       if(mcUexE > 0) fhMCPtHbpXEUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
524       
525       fhMCPtZTUeCharged[7]->Fill(mcTrigPt,mcUezT);
526       if(mcUezT > 0) fhMCPtHbpZTUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
527       
528       fhMCUePart[7]->Fill(mcTrigPt);
529     }
530   }
531
532   if(fMakeSeveralUE)
533   {
534     // Left
535     if((mcdeltaPhi<-fUeDeltaPhiMinCut) || (mcdeltaPhi >2*fUeDeltaPhiMaxCut))
536     {
537       Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
538       Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
539       Double_t mcUezT =   mcAssocPt/mcTrigPt;
540       
541       if(mcUexE < 0.)
542         printf("FillChargedMCCorrelationHistograms(): Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
543                mcUexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
544       
545       fhMCPtXEUeLeftCharged[histoIndex]->Fill(mcTrigPt,mcUexE);
546       if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
547       
548       fhMCPtZTUeLeftCharged[histoIndex]->Fill(mcTrigPt,mcUezT);
549       if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
550       
551       if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
552       {
553         fhMCPtXEUeLeftCharged[7]->Fill(mcTrigPt,mcUexE);
554         if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[7]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
555         
556         fhMCPtZTUeLeftCharged[7]->Fill(mcTrigPt,mcUezT);
557         if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[7]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
558       }
559     }
560   }
561   
562   return lead;
563 }
564
565 //___________________________________________________________________________________________________________________
566 void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Float_t ptTrig,   Float_t ptAssoc, 
567                                                                              Float_t deltaPhi,
568                                                                              Int_t   cen,      Int_t   charge,
569                                                                              Int_t   bin,      Int_t   decayTag,
570                                                                              Int_t   outTOF,   Int_t   mcTag)
571
572 {
573   // Fill mostly momentum imbalance related histograms
574   
575   Float_t zT =   ptAssoc/ptTrig ;
576   Float_t xE   =-ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
577   Float_t pout = ptAssoc*TMath::Sin(deltaPhi) ;
578
579   if(xE < 0.)
580     printf("FillChargedMomentumImbalanceHistograms(): Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
581            xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
582
583   Float_t hbpXE = -100;
584   Float_t hbpZT = -100;
585
586   if(xE > 0 ) hbpXE = TMath::Log(1./xE);
587   if(zT > 0 ) hbpZT = TMath::Log(1./zT);
588   
589   fhXECharged         ->Fill(ptTrig , xE);
590   fhPtHbpXECharged    ->Fill(ptTrig , hbpXE);
591   fhZTCharged         ->Fill(ptTrig , zT); 
592   fhPtHbpZTCharged    ->Fill(ptTrig , hbpZT);
593   fhPtTrigPout        ->Fill(ptTrig , pout) ;
594   fhPtTrigCharged     ->Fill(ptTrig , ptAssoc) ;
595   if((deltaPhi > 5*TMath::Pi()/6.)   && (deltaPhi < 7*TMath::Pi()/6.))
596   {
597     fhXECharged_Cone2         ->Fill(ptTrig , xE);
598     fhPtHbpXECharged_Cone2    ->Fill(ptTrig , hbpXE);
599   }
600   
601   // MC
602   if(IsDataMC())
603   {
604     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
605     fhXEChargedMC[mcIndex]->Fill(ptTrig , xE);
606     if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==2 )
607       fhXEChargedMC[7]->Fill(ptTrig , xE);
608   }
609   
610   // Pile up studies
611   if(IsPileUpAnalysisOn())
612   {
613     if     (outTOF==1)
614     {
615       fhXEChargedOtherBC    ->Fill(ptTrig,xE);
616       fhZTChargedOtherBC    ->Fill(ptTrig,zT);
617       fhPtTrigChargedOtherBC->Fill(ptTrig,ptAssoc);
618     }
619     else if(outTOF==0)
620     {
621       fhXEChargedBC0    ->Fill(ptTrig,xE);
622       fhZTChargedBC0    ->Fill(ptTrig,zT);
623       fhPtTrigChargedBC0->Fill(ptTrig,ptAssoc);
624     }
625
626     Int_t vtxBC = GetReader()->GetVertexBC();
627     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
628     {
629       fhXEChargedVtxBC0    ->Fill(ptTrig,xE);
630       fhZTChargedVtxBC0    ->Fill(ptTrig,zT);
631       fhPtTrigChargedVtxBC0->Fill(ptTrig,ptAssoc);
632     }
633        
634     if(GetReader()->IsPileUpFromSPD())                { fhXEChargedPileUp[0]->Fill(ptTrig,xE); fhZTChargedPileUp[0]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[0]->Fill(ptTrig,ptAssoc); }
635     if(GetReader()->IsPileUpFromEMCal())              { fhXEChargedPileUp[1]->Fill(ptTrig,xE); fhZTChargedPileUp[1]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[1]->Fill(ptTrig,ptAssoc); }
636     if(GetReader()->IsPileUpFromSPDOrEMCal())         { fhXEChargedPileUp[2]->Fill(ptTrig,xE); fhZTChargedPileUp[2]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[2]->Fill(ptTrig,ptAssoc); }
637     if(GetReader()->IsPileUpFromSPDAndEMCal())        { fhXEChargedPileUp[3]->Fill(ptTrig,xE); fhZTChargedPileUp[3]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[3]->Fill(ptTrig,ptAssoc); }
638     if(GetReader()->IsPileUpFromSPDAndNotEMCal())     { fhXEChargedPileUp[4]->Fill(ptTrig,xE); fhZTChargedPileUp[4]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[4]->Fill(ptTrig,ptAssoc); }
639     if(GetReader()->IsPileUpFromEMCalAndNotSPD())     { fhXEChargedPileUp[5]->Fill(ptTrig,xE); fhZTChargedPileUp[5]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[5]->Fill(ptTrig,ptAssoc); }
640     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhXEChargedPileUp[6]->Fill(ptTrig,xE); fhZTChargedPileUp[6]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[6]->Fill(ptTrig,ptAssoc); }
641   }
642   
643   if(fDecayTrigger && decayTag > 0)
644   {
645     for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
646     {
647       if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
648       {
649         fhXEDecayCharged[ibit]->Fill(ptTrig,xE);
650         fhZTDecayCharged[ibit]->Fill(ptTrig,zT);
651       }
652     }
653   } // photon decay pi0/eta trigger
654   
655   if(bin >= 0 && fFillMomImbalancePtAssocBinsHisto)//away side
656   {
657     fhXEAssocPtBin[bin]->Fill(ptTrig, xE) ;
658     fhZTAssocPtBin[bin]->Fill(ptTrig, zT) ;
659   }        
660   
661   if(fCorrelVzBin)
662   {
663     Int_t vz = GetEventVzBin();
664     fhXEVZ[vz]->Fill(ptTrig, xE) ;
665     fhZTVZ[vz]->Fill(ptTrig, zT) ;
666   }
667   
668   if(charge > 0)
669   {
670     fhXEPosCharged->Fill(ptTrig,xE) ;
671     fhZTPosCharged->Fill(ptTrig,zT) ;
672   }
673   else  
674   {
675     fhXENegCharged->Fill(ptTrig,xE) ;
676     fhZTNegCharged->Fill(ptTrig,zT) ;
677   }
678   
679   //fill different multiplicity/centrality histogram
680   if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
681   {
682     fhXEMult[cen]->Fill(ptTrig,xE);
683     fhZTMult[cen]->Fill(ptTrig,zT);
684   } //multiplicity/centrality events selection
685
686
687 //_______________________________________________________________________________________________________________________
688 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(Float_t ptTrig,   Float_t ptAssoc,
689                                                                            Float_t deltaPhi, Int_t cen, Int_t outTOF)
690 {
691   // Fill underlying event histograms
692   
693   fhUePart->Fill(ptTrig);
694   
695   fhDeltaPhiUeChargedPt->Fill(ptAssoc,deltaPhi);
696   
697   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
698   Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
699   Double_t uezT =   ptAssoc/ptTrig;
700   
701   if(uexE < 0.)
702     printf("FillChargedUnderlyingEventHistograms(): Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
703            uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
704   
705   fhXEUeCharged->Fill(ptTrig,uexE);
706   if(uexE > 0) fhPtHbpXEUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
707   
708   fhZTUeCharged->Fill(ptTrig,uezT);
709   if(uezT > 0) fhPtHbpZTUeCharged->Fill(ptTrig,TMath::Log(1/uezT));
710   
711   // Pile up studies
712   
713   if(IsPileUpAnalysisOn())
714   {
715     if     (outTOF==1)
716     {
717       fhXEUeChargedOtherBC->Fill(ptTrig,uexE);
718       fhZTUeChargedOtherBC->Fill(ptTrig,uezT);
719     }
720     else if(outTOF==0)
721     {
722       fhXEUeChargedBC0->Fill(ptTrig,uexE);
723       fhZTUeChargedBC0->Fill(ptTrig,uezT);
724     }
725     
726     Int_t vtxBC = GetReader()->GetVertexBC();
727     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
728     {
729       fhXEUeChargedVtxBC0->Fill(ptTrig,uexE);
730       fhZTUeChargedVtxBC0->Fill(ptTrig,uezT);
731     }
732
733     if(GetReader()->IsPileUpFromSPD())               { fhXEUeChargedPileUp[0]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[0]->Fill(ptTrig,uezT);}
734     if(GetReader()->IsPileUpFromEMCal())             { fhXEUeChargedPileUp[1]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[1]->Fill(ptTrig,uezT);}
735     if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhXEUeChargedPileUp[2]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[2]->Fill(ptTrig,uezT);}
736     if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhXEUeChargedPileUp[3]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[3]->Fill(ptTrig,uezT);}
737     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhXEUeChargedPileUp[4]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[4]->Fill(ptTrig,uezT);}
738     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhXEUeChargedPileUp[5]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[5]->Fill(ptTrig,uezT);}
739     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhXEUeChargedPileUp[6]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[6]->Fill(ptTrig,uezT);}
740   }
741   
742   //fill different multiplicity/centrality histogram
743   if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
744   {
745     fhXEUeMult[cen]->Fill(ptTrig,uexE);
746     fhZTUeMult[cen]->Fill(ptTrig,uezT);
747   }
748 }
749
750 //_____________________________________________________________________________________________________
751 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(Float_t ptTrig, 
752                                                                                 Float_t ptAssoc, 
753                                                                                 Float_t deltaPhi)
754 {
755   // Fill underlying event histograms to the left and right of trigger
756   // Right cone is the default UE.
757   
758   if((deltaPhi<-fUeDeltaPhiMinCut) || (deltaPhi >2*fUeDeltaPhiMaxCut))
759   {
760     fhDeltaPhiUeLeftCharged->Fill(ptAssoc,deltaPhi);
761     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
762     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
763     Double_t uezT =   ptAssoc/ptTrig;
764     
765     if(uexE < 0.)
766       printf("FillChargedUnderlyingEventSidesHistograms(): Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
767               uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
768     
769     fhXEUeLeftCharged->Fill(ptTrig,uexE);
770     if(uexE > 0) fhPtHbpXEUeLeftCharged->Fill(ptTrig,TMath::Log(1/uexE));
771     
772     fhZTUeLeftCharged->Fill(ptTrig,uezT);
773     if(uezT > 0) fhPtHbpZTUeLeftCharged->Fill(ptTrig,TMath::Log(1/uezT));
774     
775     fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi);
776   }
777   
778   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-TMath::Pi()/2))
779   {
780     fhDeltaPhiUeLeftDownCharged->Fill(ptAssoc,deltaPhi);
781     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
782     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
783     
784     if(uexE < 0.)
785       printf("FillChargedUnderlyingEventSidesHistograms(): Careful!!, negative xE %2.2f for left-down UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
786              uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
787     
788     fhXEUeLeftDownCharged->Fill(ptTrig,uexE);
789   }
790   
791   if((deltaPhi>2*fUeDeltaPhiMaxCut) && (deltaPhi <3*TMath::Pi()/2))
792   {
793     fhDeltaPhiUeLeftUpCharged->Fill(ptAssoc,deltaPhi);
794     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
795     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
796     
797     if(uexE < 0.)
798       printf("FillChargedUnderlyingEventSidesHistograms(): Careful!!, negative xE %2.2f for left-up UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
799              uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
800     
801     fhXEUeLeftUpCharged->Fill(ptTrig,uexE);
802   }
803   
804   if((deltaPhi > TMath::Pi()/2) && (deltaPhi < fUeDeltaPhiMaxCut))
805   {
806     fhDeltaPhiUeRightUpCharged->Fill(ptAssoc,deltaPhi);
807     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
808     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
809     
810     if(uexE < 0.)
811       printf("FillChargedUnderlyingEventSidesHistograms(): Careful!!, negative xE %2.2f for right-up UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
812              uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
813     
814     fhXEUeRightUpCharged->Fill(ptTrig,uexE);
815   }
816   
817   if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < TMath::Pi()/2))
818   {
819     fhDeltaPhiUeRightDownCharged->Fill(ptAssoc,deltaPhi);
820     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
821     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
822     
823     if(uexE < 0.)
824       printf("FillChargedUnderlyingEventSidesHistograms(): Careful!!, negative xE %2.2f for right-down UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f\n",
825              uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg());
826     
827     fhXEUeRightDownCharged->Fill(ptTrig,uexE);
828   }
829
830
831 //_____________________________________________________________________________________________________________________________________
832 void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float_t ptAssoc, Float_t phiAssoc, Bool_t bChargedOrNeutral)
833 {
834   // Do correlation with decay photons of triggered pi0 or eta
835   
836   // Calculate the correlation parameters
837   Float_t ptDecay1 = fDecayMom1.Pt();
838   Float_t ptDecay2 = fDecayMom2.Pt();
839   
840   Float_t zTDecay1 = -100, zTDecay2 = -100;
841   if(ptDecay1 > 0) zTDecay1 = ptAssoc/ptDecay1 ;
842   if(ptDecay2 > 0) zTDecay2 = ptAssoc/ptDecay2 ;
843   
844   Float_t deltaPhiDecay1 = fDecayMom1.Phi()-phiAssoc;
845   if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
846   if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
847   
848   Float_t deltaPhiDecay2 = fDecayMom2.Phi()-phiAssoc;
849   if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
850   if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
851   
852   Float_t xEDecay1   =-zTDecay1*TMath::Cos(deltaPhiDecay1); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
853   Float_t xEDecay2   =-zTDecay2*TMath::Cos(deltaPhiDecay2); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
854   
855   if(bChargedOrNeutral) // correlate with charges
856   {
857     fhDeltaPhiPi0DecayCharged->Fill(ptDecay1, deltaPhiDecay1);
858     fhDeltaPhiPi0DecayCharged->Fill(ptDecay2, deltaPhiDecay2);
859     
860     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms( Charged corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
861     
862     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
863     {
864       fhZTPi0DecayCharged->Fill(ptDecay1,zTDecay1); 
865       fhXEPi0DecayCharged->Fill(ptDecay1,xEDecay1); 
866     }
867     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
868     {
869       fhZTPi0DecayCharged->Fill(ptDecay2,zTDecay2);
870       fhXEPi0DecayCharged->Fill(ptDecay2,xEDecay2);
871     }
872   }
873   else // correlate with neutrals
874   {
875     fhDeltaPhiPi0DecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
876     fhDeltaPhiPi0DecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
877     
878     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms(Neutral corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
879     
880     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
881     {
882       fhZTPi0DecayNeutral->Fill(ptDecay1,zTDecay1);
883       fhXEPi0DecayNeutral->Fill(ptDecay1,xEDecay1);
884     }
885     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
886     {
887       fhZTPi0DecayNeutral->Fill(ptDecay2,zTDecay2);
888       fhXEPi0DecayNeutral->Fill(ptDecay2,xEDecay2);
889     }    
890   }
891 }
892
893 //_____________________________________________________________________________________________________________________________
894 void AliAnaParticleHadronCorrelation::FillNeutralUnderlyingEventSidesHistograms(Float_t ptTrig,   Float_t ptAssoc, 
895                                                                                 Float_t zT,       Float_t hbpZT,
896                                                                                 Float_t deltaPhi)
897 {
898   // Fill underlying event histograms to the left of trigger
899   // Right is the default case
900   
901   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
902
903   Float_t xE  =-ptAssoc/ptTrig*TMath::Cos(randomphi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
904   Float_t hbpXE = -100;
905   if(xE > 0 ) hbpXE = TMath::Log(1./xE);
906   
907   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
908   {
909     fhDeltaPhiUeLeftNeutral->Fill(ptAssoc, deltaPhi);
910     fhXEUeLeftNeutral      ->Fill(ptTrig , xE);
911     fhPtHbpXEUeLeftNeutral ->Fill(ptTrig , hbpXE);
912     fhZTUeLeftNeutral      ->Fill(ptTrig , zT);
913     fhPtHbpZTUeLeftNeutral ->Fill(ptTrig , hbpZT);
914   }
915
916
917 //______________________________________________________
918 void AliAnaParticleHadronCorrelation::FillEventMixPool()
919 {
920   // Fill the pool with tracks or clusters if requested
921     
922   if ( !DoOwnMix() ) return;
923   
924   FillChargedEventMixPool();
925   
926   // Do the cluster pool filling only if requested
927   // or in case of isolation cut using clusters in the cone.
928   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
929   
930   if( !fFillNeutralEventMixPool && !isoCase) return;
931   
932   FillNeutralEventMixPool();
933 }
934
935 //_____________________________________________________________
936 void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
937 {
938   // Mixed event pool filling for tracks
939     
940   if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
941   {
942     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
943     return ; // pool filled previously for another trigger
944   }
945   
946   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
947   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
948   
949   if(!inputHandler) return ;
950     
951   // Do mixing only with MB event (or the chosen mask), if not skip
952   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
953   
954   Int_t eventBin = GetEventMixBin();
955   
956   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
957   if(eventBin < 0) return;
958   
959   fhEventMBBin->Fill(eventBin);
960   
961   TObjArray * mixEventTracks = new TObjArray;
962   
963   if(fUseMixStoredInReader) 
964   {
965     fListMixTrackEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
966   }
967   
968   if(!fListMixTrackEvents[eventBin]) fListMixTrackEvents[eventBin] = new TList();
969   
970   //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
971   
972   TList * pool = fListMixTrackEvents[eventBin];
973   
974   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
975   {
976     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
977     
978     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
979     Float_t pt   = fTrackVector.Pt();
980     
981     //Select only hadrons in pt range
982     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
983     
984     AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),0);
985     mixedTrack->SetDetectorTag(kCTS);
986     mixedTrack->SetChargedBit(track->Charge()>0);
987     mixEventTracks->Add(mixedTrack);
988   }
989   
990   fhNtracksMB->Fill(mixEventTracks->GetEntriesFast(),eventBin);
991   
992   //Set the event number where the last event was added, to avoid double pool filling
993   GetReader()->SetLastTracksMixedEvent(GetEventNumber());
994   
995   //printf("Add event to pool with %d tracks \n ",mixEventTracks->GetEntries());
996   pool->AddFirst(mixEventTracks);
997   mixEventTracks = 0;
998   
999   //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
1000
1001   if(pool->GetSize() > GetNMaxEvMix())
1002   {//Remove last event
1003     TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
1004     pool->RemoveLast() ;
1005     delete tmp ;
1006   }
1007 }
1008
1009 //_____________________________________________________________
1010 void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
1011 {
1012   // Mixed event pool filling for neutral clusters
1013   // Right now only for EMCAL and in isolation case
1014   
1015   //printf("FillNeutralEventMixPool for %s\n",GetInputAODName().Data());
1016   
1017   if(fUseMixStoredInReader && GetReader()->GetLastCaloMixedEvent() == GetEventNumber())
1018   {
1019     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
1020     return ; // pool filled previously for another trigger
1021   }
1022   
1023   TObjArray * pl = GetEMCALClusters();
1024   //if (GetAODObjArrayName.Contains("PHOS") )pl    = GetPHOSClusters();
1025   //else                                     pl    = GetEMCALClusters();
1026   
1027   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
1028   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
1029   
1030   if(!inputHandler) return ;
1031     
1032   // Do mixing only with MB event (or the chosen mask), if not skip
1033   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
1034   
1035   Int_t eventBin = GetEventMixBin();
1036   
1037   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
1038   if(eventBin < 0) return;
1039   
1040   TObjArray * mixEventCalo = new TObjArray;
1041   
1042   if(fUseMixStoredInReader) 
1043   {
1044     fListMixCaloEvents[eventBin] = GetReader()->GetListWithMixedEventsForCalo(eventBin);
1045   }
1046   
1047   if(!fListMixCaloEvents[eventBin]) fListMixCaloEvents[eventBin] = new TList();
1048   
1049   TList * poolCalo = fListMixCaloEvents[eventBin];
1050   
1051   for(Int_t ipr = 0;ipr <  pl->GetEntriesFast() ; ipr ++ )
1052   {
1053     AliVCluster * calo = (AliVCluster *) (pl->At(ipr)) ;
1054   
1055     // remove matched clusters
1056     if( IsTrackMatched( calo, GetReader()->GetInputEvent() ) ) continue ;
1057     
1058     //Cluster momentum calculation
1059     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1060     {
1061       calo->GetMomentum(fMomentum,GetVertex(0)) ;
1062     }//Assume that come from vertex in straight line
1063     else
1064     {
1065       Double_t vertex[]={0,0,0};
1066       calo->GetMomentum(fMomentum,vertex) ;
1067     }
1068     
1069     Float_t pt = fMomentum.Pt();
1070     //Select only clusters in pt range
1071     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1072     
1073     AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(fMomentum);
1074     mixedCalo->SetDetectorTag(kEMCAL);
1075     mixEventCalo->Add(mixedCalo);
1076   }
1077   
1078   fhNclustersMB->Fill(mixEventCalo->GetEntriesFast(),eventBin);
1079   
1080   //Set the event number where the last event was added, to avoid double pool filling
1081   GetReader()->SetLastCaloMixedEvent(GetEventNumber());
1082   
1083   //printf("Add event to pool with %d clusters \n ",mixEventCalo->GetEntries());
1084   poolCalo->AddFirst(mixEventCalo);
1085   mixEventCalo = 0;
1086   
1087   //printf("Pool size %d, max %d\n",poolCalo->GetSize(), GetNMaxEvMix());
1088   
1089   if(poolCalo->GetSize() > GetNMaxEvMix())
1090   {//Remove last event
1091     TClonesArray * tmp = static_cast<TClonesArray*>(poolCalo->Last()) ;
1092     poolCalo->RemoveLast() ;
1093     delete tmp ;
1094   }  
1095 }
1096
1097 //_________________________________________________________________________________________________________________
1098 Bool_t AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow(AliAODPWG4ParticleCorrelation * particle)
1099 {
1100   // Select events where the leading charged particle in the opposite hemisphere
1101   // to the trigger particle is in a window centered at 180 from the trigger
1102   
1103   Float_t etaTrig    = particle->Eta();
1104   Float_t ptTrig     = particle->Pt();
1105   Float_t phiTrig    = particle->Phi();
1106   if(phiTrig < 0 ) phiTrig+= TMath::TwoPi();
1107
1108   Float_t ptLeadHad  = 0 ;
1109   Float_t dphiLeadHad= -100 ;
1110   Float_t phiLeadHad = -100 ;
1111   Float_t etaLeadHad = -100 ;
1112   Int_t   nTrack     = 0;
1113
1114   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
1115   {
1116     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1117     
1118     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1119     
1120     Float_t pt   = fTrackVector.Pt();
1121     Float_t phi  = fTrackVector.Phi() ;
1122     if(phi < 0 ) phi+= TMath::TwoPi();
1123     
1124     Float_t deltaPhi = phiTrig-phi;
1125     //
1126     // Calculate deltaPhi shift so that for the particles on the opposite side
1127     // it is defined between 90 and 270 degrees
1128     // Shift [-360,-90]  to [0, 270]
1129     // and [270,360] to [-90,0]
1130     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1131     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
1132
1133     if(pt > ptLeadHad && deltaPhi > TMath::PiOver2()) // in opposite hemisphere
1134     {
1135       ptLeadHad  = pt ;
1136       phiLeadHad = phi;
1137       dphiLeadHad= deltaPhi;
1138       etaLeadHad = fTrackVector.Eta();
1139       nTrack++;
1140     }
1141   }// track loop
1142   
1143   if(fFillLeadHadOppositeHisto)
1144   {
1145     if(nTrack == 0)
1146     {
1147       fhPtNoLeadingOppositeHadron    ->Fill(ptTrig);
1148       fhEtaPhiNoLeadingOppositeHadron->Fill(etaTrig,phiTrig);
1149     }
1150     else
1151     {
1152       fhPtLeadingOppositeHadron       ->Fill(ptTrig,  ptLeadHad);
1153       fhPtDiffPhiLeadingOppositeHadron->Fill(ptTrig,dphiLeadHad);
1154       fhPtDiffEtaLeadingOppositeHadron->Fill(ptTrig, etaLeadHad-etaTrig);
1155     }
1156   }
1157   
1158   if(GetDebug() > 1 )
1159   {
1160     printf("AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow() pT %2.2f, phi %2.2f, eta %2.2f, nTracks away %d, total tracks %d\n",
1161            ptLeadHad,phiLeadHad*TMath::RadToDeg(),etaLeadHad,nTrack, GetTrackMultiplicity());
1162     
1163     printf("\t  pT trig %2.2f, Dphi (trigger-hadron) %2.2f, Deta (trigger-hadron) %2.2f\n",
1164            ptTrig, dphiLeadHad*TMath::RadToDeg(), etaLeadHad-etaTrig);
1165     printf("\t cuts pT: min %2.2f, max %2.2f; DPhi: min %2.2f, max %2.2f\n",fMinLeadHadPt,fMaxLeadHadPt,fMinLeadHadPhi*TMath::RadToDeg(),fMaxLeadHadPhi*TMath::RadToDeg());
1166   }
1167   
1168   // reject the trigger if the leading hadron is not in the requested pt or phi window and
1169   
1170   if( nTrack == 0 ) return kFALSE; // No track found in opposite hemisphere
1171   
1172   if( ptLeadHad < fMinLeadHadPt || ptLeadHad > fMaxLeadHadPt ) return kFALSE;
1173   
1174   //printf("Accept leading hadron pT \n");
1175   
1176   if( dphiLeadHad < fMinLeadHadPhi || dphiLeadHad > fMaxLeadHadPhi ) return kFALSE;
1177   
1178   //printf("Accept leading hadron phi \n");
1179   
1180   
1181   return kTRUE ;
1182 }
1183
1184 //____________________________________________________________
1185 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
1186 {
1187   //Save parameters used for analysis
1188   TString parList ; //this will be list of parameters used for this analysis.
1189   const Int_t buffersize = 560;
1190   char onePar[buffersize] ;
1191   
1192   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---:") ;
1193   parList+=onePar ;     
1194   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f; ", fMinAssocPt,   fMaxAssocPt) ;
1195   parList+=onePar ;
1196   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f; ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
1197   parList+=onePar ;
1198   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f; ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
1199   parList+=onePar ;
1200   snprintf(onePar,buffersize,"Isolated Trigger?  %d;",    fSelectIsolated) ;
1201   parList+=onePar ;
1202   snprintf(onePar,buffersize,"Several UE?  %d;",          fMakeSeveralUE) ;
1203   parList+=onePar ;
1204   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s;", fPi0AODBranchName.Data());
1205   parList+=onePar ;
1206   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d;", fPi0Trigger, fDecayTrigger) ;
1207   parList+=onePar ;
1208   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d;",
1209            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
1210   parList+=onePar ;
1211   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
1212   parList+=onePar ;
1213   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
1214     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
1215   }
1216   parList+=onePar ;
1217   
1218   //Get parameters set in base class.
1219   parList += GetBaseParametersList() ;
1220   
1221   //Get parameters set in FiducialCut class (not available yet)
1222   //parlist += GetFidCut()->GetFidCutParametersList() 
1223   
1224   return new TObjString(parList) ;  
1225   
1226
1227
1228 //________________________________________________________________
1229 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
1230 {
1231   // Create histograms to be saved in output file and
1232   // store them in fOutputContainer
1233   
1234   TList * outputContainer = new TList() ;
1235   outputContainer->SetName("CorrelationHistos") ;
1236   
1237   Int_t   nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t  nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
1238   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
1239   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
1240
1241   Int_t  ndeltaphibins = GetHistogramRanges()->GetHistoDeltaPhiBins(); Int_t   ndeltaetabins = GetHistogramRanges()->GetHistoDeltaEtaBins();
1242   Float_t deltaphimax  = GetHistogramRanges()->GetHistoDeltaPhiMax();  Float_t deltaetamax   = GetHistogramRanges()->GetHistoDeltaEtaMax();
1243   Float_t deltaphimin  = GetHistogramRanges()->GetHistoDeltaPhiMin();  Float_t deltaetamin   = GetHistogramRanges()->GetHistoDeltaEtaMin();
1244
1245   Int_t ntrbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins();
1246   Int_t trmax   = GetHistogramRanges()->GetHistoTrackMultiplicityMax();  Int_t clmax   = GetHistogramRanges()->GetHistoNClustersMax();
1247   Int_t trmin   = GetHistogramRanges()->GetHistoTrackMultiplicityMin();  Int_t clmin   = GetHistogramRanges()->GetHistoNClustersMin();
1248
1249   Int_t  nxeztbins = GetHistogramRanges()->GetHistoRatioBins();  Int_t  nhbpbins = GetHistogramRanges()->GetHistoHBPBins();
1250   Float_t xeztmax  = GetHistogramRanges()->GetHistoRatioMax();   Float_t hbpmax  = GetHistogramRanges()->GetHistoHBPMax();
1251   Float_t xeztmin  = GetHistogramRanges()->GetHistoRatioMin();   Float_t hbpmin  = GetHistogramRanges()->GetHistoHBPMin();
1252   
1253   Int_t nMixBins = GetNCentrBin()*GetNZvertBin()*GetNRPBin();
1254   
1255   TString nameMC[]     = {"Photon","Pi0","Pi0Decay","EtaDecay","OtherDecay","Electron","Hadron","Pi0DecayLostPair"};
1256   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1257   
1258   // For vz dependent histograms, if option ON
1259   Int_t   nz  = 1  ;
1260   if(fCorrelVzBin) nz = GetNZvertBin();
1261   TString sz  = "" ;
1262   TString tz  = "" ;
1263   
1264   // Fill histograms for neutral clusters in mixing?
1265   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
1266   Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
1267   
1268   fhPtTriggerInput  = new TH1F("hPtTriggerInput","Input trigger #it{p}_{T}", nptbins,ptmin,ptmax);
1269   fhPtTriggerInput->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1270   outputContainer->Add(fhPtTriggerInput);
1271   
1272   if( fM02MaxCut > 0 && fM02MinCut > 0 )
1273   {
1274     fhPtTriggerSSCut  = new TH1F("hPtTriggerSSCut","Trigger #it{p}_{T} after #lambda^{2}_{0} cut", nptbins,ptmin,ptmax);
1275     fhPtTriggerSSCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1276     outputContainer->Add(fhPtTriggerSSCut);
1277   }
1278   
1279   if( OnlyIsolated() )
1280   {
1281     fhPtTriggerIsoCut  = new TH1F("hPtTriggerIsoCut","Trigger #it{p}_{T} after isolation (and #lambda^{2}_{0} cut)", nptbins,ptmin,ptmax);
1282     fhPtTriggerIsoCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1283     outputContainer->Add(fhPtTriggerIsoCut);
1284   }
1285   
1286   fhPtTriggerFidCut  = new TH1F("hPtTriggerFidCut","Trigger #it{p}_{T} after fiducial (isolation and #lambda^{2}_{0}) cut", nptbins,ptmin,ptmax);
1287   fhPtTriggerFidCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1288   outputContainer->Add(fhPtTriggerFidCut);
1289   
1290   fhPtTrigger  = new TH1F("hPtTrigger","#it{p}_{T} distribution of trigger particles (after opposite hadron leading cut and rest)", nptbins,ptmin,ptmax);
1291   fhPtTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1292   outputContainer->Add(fhPtTrigger);
1293   
1294   if(IsDataMC())
1295   {
1296     for(Int_t i=0; i < fgkNmcTypes; i++)
1297     {
1298       fhPtTriggerMC[i]  = new TH1F(Form("hPtTrigger_MC%s",nameMC[i].Data()),
1299                                    Form("#it{p}_{T} distribution of trigger particles, trigger origin is %s",nameMC[i].Data()),
1300                                    nptbins,ptmin,ptmax);
1301       fhPtTriggerMC[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1302       outputContainer->Add(fhPtTriggerMC[i]);
1303     }
1304   }
1305   
1306   if(fDecayTrigger)
1307   {
1308     for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1309     {
1310       fhPtDecayTrigger[ibit]  = new TH1F(Form("hPtDecayTrigger_bit%d",fDecayBits[ibit]),
1311                                          Form("#it{p}_{T} distribution of trigger particles, decay Bit %d",fDecayBits[ibit]),
1312                                          nptbins,ptmin,ptmax);
1313       fhPtDecayTrigger[ibit]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1314       outputContainer->Add(fhPtDecayTrigger[ibit]);
1315       
1316       if(IsDataMC())
1317       {
1318         for(Int_t i=0; i < fgkNmcTypes; i++)
1319         {
1320           fhPtDecayTriggerMC[ibit][i]  = new TH1F(Form("hPtDecayTrigger_bit%d_MC%s",fDecayBits[ibit], nameMC[i].Data()),
1321                                             Form("#it{p}_{T} distribution of trigger particles, decay Bit %d, trigger origin is %s",fDecayBits[ibit], nameMC[i].Data()),
1322                                             nptbins,ptmin,ptmax);
1323           fhPtDecayTriggerMC[ibit][i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1324           outputContainer->Add(fhPtDecayTriggerMC[ibit][i]);
1325         }
1326       }
1327     }
1328   }
1329   
1330   if(fCorrelVzBin)
1331   {
1332     fhPtTriggerVzBin  = new TH2F("hPtTriggerVzBin","#it{p}_{T} distribution of trigger particles vs vz bin", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
1333     fhPtTriggerVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1334     fhPtTriggerVzBin->SetYTitle("#it{v}_{#it{z}} bin");
1335     outputContainer->Add(fhPtTriggerVzBin);
1336   }
1337   
1338   fhPtTriggerBin  = new TH2F ("hPtTriggerBin","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
1339   fhPtTriggerBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1340   fhPtTriggerBin->SetYTitle("Bin");
1341   outputContainer->Add(fhPtTriggerBin);
1342   
1343   fhPhiTrigger  = new TH2F ("hPhiTrigger","#phi distribution of trigger Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1344   fhPhiTrigger->SetYTitle("#phi (rad)");
1345   outputContainer->Add(fhPhiTrigger);
1346   
1347   fhEtaTrigger  = new TH2F ("hEtaTrigger","#eta distribution of trigger",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1348   fhEtaTrigger->SetYTitle("#eta ");
1349   outputContainer->Add(fhEtaTrigger);
1350   
1351   if(IsHighMultiplicityAnalysisOn())
1352   {
1353     fhPtTriggerCentrality   = new TH2F("hPtTriggerCentrality","Trigger particle #it{p}_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
1354     fhPtTriggerCentrality->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1355     fhPtTriggerCentrality->SetYTitle("Centrality (%)");
1356     outputContainer->Add(fhPtTriggerCentrality) ;
1357     
1358     fhPtTriggerEventPlane  = new TH2F("hPtTriggerEventPlane","Trigger particle #it{p}_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
1359     fhPtTriggerEventPlane->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1360     fhPtTriggerEventPlane->SetXTitle("EP angle (rad)");
1361     outputContainer->Add(fhPtTriggerEventPlane) ;
1362     
1363     fhTriggerEventPlaneCentrality  = new TH2F("hTriggerEventPlaneCentrality","Trigger particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
1364     fhTriggerEventPlaneCentrality->SetXTitle("Centrality (%)");
1365     fhTriggerEventPlaneCentrality->SetYTitle("EP angle (rad)");
1366     outputContainer->Add(fhTriggerEventPlaneCentrality) ;
1367   }
1368
1369   // Leading hadron in oposite side
1370   if(fFillLeadHadOppositeHisto)
1371   {
1372     fhPtLeadingOppositeHadron  = new TH2F("hPtTriggerPtLeadingOppositeHadron","Leading hadron opposite to trigger vs trigger #it{p}_{T}",
1373                                           nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1374     fhPtLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1375     fhPtLeadingOppositeHadron->SetYTitle("#it{p}_{T}^{lead hadron} (GeV/#it{c})");
1376     outputContainer->Add(fhPtLeadingOppositeHadron);
1377
1378     fhPtNoLeadingOppositeHadron  = new TH1F("hPtTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #it{p}_{T}",
1379                                           nptbins,ptmin,ptmax);
1380     fhPtNoLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1381     outputContainer->Add(fhPtNoLeadingOppositeHadron);
1382
1383     fhEtaPhiNoLeadingOppositeHadron  = new TH2F("hEtaPhiTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #eta:#phi",
1384                                             netabins,etamin,etamax,nphibins,phimin,phimax);
1385     fhEtaPhiNoLeadingOppositeHadron->SetXTitle("#eta");
1386     fhEtaPhiNoLeadingOppositeHadron->SetYTitle("#phi");
1387     outputContainer->Add(fhEtaPhiNoLeadingOppositeHadron);
1388
1389     
1390     fhPtDiffPhiLeadingOppositeHadron  = new TH2F("hPtTriggerDiffPhiTriggerLeadingOppositeHadron","#phi_{trigger}-#phi_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1391                                                  nptbins,ptmin,ptmax,ndeltaphibins,deltaphimin,deltaphimax);
1392     fhPtDiffPhiLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1393     fhPtDiffPhiLeadingOppositeHadron->SetYTitle("#phi_{trigger}-#phi_{leading opposite hadron} (rad)");
1394     outputContainer->Add(fhPtDiffPhiLeadingOppositeHadron);
1395     
1396     fhPtDiffEtaLeadingOppositeHadron  = new TH2F("hPtTriggerDiffEtaTriggerPhiLeadingOppositeHadron","#eta_{trigger}-#eta_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1397                                                  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1398     fhPtDiffEtaLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1399     fhPtDiffEtaLeadingOppositeHadron->SetYTitle("#eta_{trigger}-#eta_{leading opposite hadron}");
1400     outputContainer->Add(fhPtDiffEtaLeadingOppositeHadron);
1401   }
1402   
1403   //Correlation with charged hadrons
1404   
1405   fhDeltaPhiDeltaEtaCharged  = new TH2F
1406   ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
1407    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1408   fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
1409   fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1410   
1411   fhDeltaPhiDeltaEtaChargedPtA3GeV  = new TH2F
1412   ("hDeltaPhiDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}, #it{p}_{TA}>3 GeV/#it{c}}",
1413    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1414   fhDeltaPhiDeltaEtaChargedPtA3GeV->SetXTitle("#Delta #phi (rad)");
1415   fhDeltaPhiDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1416   
1417   fhPhiCharged  = new TH2F
1418   ("hPhiCharged","#phi_{h^{#pm}}  vs #it{p}_{T #pm}",
1419    nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
1420   fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
1421   fhPhiCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1422   
1423   fhEtaCharged  = new TH2F
1424   ("hEtaCharged","#eta_{h^{#pm}}  vs #it{p}_{T #pm}",
1425    nptbins,ptmin,ptmax,100,-1.,1.);
1426   fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1427   fhEtaCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1428   
1429   fhDeltaPhiCharged  = new TH2F
1430   ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
1431    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1432   fhDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
1433   fhDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1434   
1435   fhDeltaPhiChargedPtA3GeV  = new TH2F
1436   ("hDeltaPhiChargedPtA3GeV","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1437    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1438   fhDeltaPhiChargedPtA3GeV->SetYTitle("#Delta #phi (rad)");
1439   fhDeltaPhiChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1440   
1441   
1442   fhDeltaPhiChargedPt  = new TH2F
1443   ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",
1444    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1445   fhDeltaPhiChargedPt->SetYTitle("#Delta #phi (rad)");
1446   fhDeltaPhiChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1447   
1448   fhDeltaEtaCharged  = new TH2F
1449   ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}",
1450    nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1451   fhDeltaEtaCharged->SetYTitle("#Delta #eta");
1452   fhDeltaEtaCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1453   
1454   fhDeltaEtaChargedPtA3GeV  = new TH2F
1455   ("hDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1456    nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1457   fhDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1458   fhDeltaEtaChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1459   
1460   fhXECharged  =
1461   new TH2F("hXECharged","#it{x}_{#it{E}} for charged tracks",
1462            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1463   fhXECharged->SetYTitle("#it{x}_{#it{E}}");
1464   fhXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1465   
1466   fhXECharged_Cone2  =
1467   new TH2F("hXECharged_Cone2","#it{x}_{#it{E}} for charged tracks in cone 2 (5#pi/6-7#pi/6)",
1468            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1469   fhXECharged_Cone2->SetYTitle("#it{x}_{#it{E}}");
1470   fhXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1471   
1472   fhXEPosCharged  =
1473   new TH2F("hXEPositiveCharged","#it{x}_{#it{E}} for positive charged tracks",
1474            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1475   fhXEPosCharged->SetYTitle("#it{x}_{#it{E}}");
1476   fhXEPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1477   
1478   fhXENegCharged  =
1479   new TH2F("hXENegativeCharged","#it{x}_{#it{E}} for negative charged tracks",
1480            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1481   fhXENegCharged->SetYTitle("#it{x}_{#it{E}}");
1482   fhXENegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1483   
1484   fhPtHbpXECharged  =
1485   new TH2F("hHbpXECharged","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
1486            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1487   fhPtHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1488   fhPtHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1489   
1490   fhPtHbpXECharged_Cone2  =
1491   new TH2F("hHbpXECharged_Cone2","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons in cone 2 (5#pi/6-7#pi/6)",
1492            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1493   fhPtHbpXECharged_Cone2->SetYTitle("ln(1/#it{x}_{#it{E}})");
1494   fhPtHbpXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1495   
1496   fhZTCharged  =
1497   new TH2F("hZTCharged","#it{z}_{T} for charged tracks",
1498            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1499   fhZTCharged->SetYTitle("#it{z}_{T}");
1500   fhZTCharged->SetXTitle("#it{p}_{T trigger}");
1501   
1502   fhZTPosCharged  =
1503   new TH2F("hZTPositiveCharged","#it{z}_{T} for positive charged tracks",
1504            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1505   fhZTPosCharged->SetYTitle("#it{z}_{T}");
1506   fhZTPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1507   
1508   fhZTNegCharged  =
1509   new TH2F("hZTNegativeCharged","#it{z}_{T} for negative charged tracks",
1510            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1511   fhZTNegCharged->SetYTitle("#it{z}_{T}");
1512   fhZTNegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1513   
1514   fhPtHbpZTCharged  =
1515   new TH2F("hHbpZTCharged","#xi = ln(1/#it{z}_{T}) with charged hadrons",
1516            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1517   fhPtHbpZTCharged->SetYTitle("ln(1/#it{z}_{T})");
1518   fhPtHbpZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1519   
1520   fhPtTrigPout  =
1521   new TH2F("hPtTrigPout","Pout with triggers",
1522            nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
1523   fhPtTrigPout->SetYTitle("#it{p}_{out} (GeV/#it{c})");
1524   fhPtTrigPout->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1525   
1526   fhPtTrigCharged  =
1527   new TH2F("hPtTrigCharged","trigger and charged tracks pt distribution",
1528            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1529   fhPtTrigCharged->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1530   fhPtTrigCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1531   
1532   outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
1533   outputContainer->Add(fhDeltaPhiDeltaEtaChargedPtA3GeV);
1534   outputContainer->Add(fhPhiCharged) ;
1535   outputContainer->Add(fhEtaCharged) ;
1536   outputContainer->Add(fhDeltaPhiCharged) ;
1537   outputContainer->Add(fhDeltaPhiChargedPtA3GeV) ;
1538   outputContainer->Add(fhDeltaEtaCharged) ;
1539   outputContainer->Add(fhDeltaEtaChargedPtA3GeV) ;
1540   outputContainer->Add(fhDeltaPhiChargedPt) ;
1541   
1542   outputContainer->Add(fhXECharged) ;
1543   outputContainer->Add(fhXECharged_Cone2) ;
1544   
1545   if(IsDataMC())
1546   {
1547     for(Int_t i=0; i < fgkNmcTypes; i++)
1548     {
1549       
1550       fhDeltaPhiChargedMC[i]  = new TH2F(Form("hDeltaPhiCharged_MC%s",nameMC[i].Data()),
1551                                          Form("#Delta #phi for charged tracks, trigger origin is %s",nameMC[i].Data()),
1552                                          nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1553       fhDeltaPhiChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1554       fhDeltaPhiChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1555       outputContainer->Add(fhDeltaPhiChargedMC[i]) ;
1556       
1557       fhXEChargedMC[i]  = new TH2F(Form("hXECharged_MC%s",nameMC[i].Data()),
1558                                    Form("#it{x}_{#it{E}} for charged tracks, trigger origin is %s",nameMC[i].Data()),
1559                                    nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1560       fhXEChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1561       fhXEChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1562       outputContainer->Add(fhXEChargedMC[i]) ;
1563     }
1564   }
1565   
1566   outputContainer->Add(fhXEPosCharged) ;
1567   outputContainer->Add(fhXENegCharged) ;
1568   outputContainer->Add(fhPtHbpXECharged) ;
1569   outputContainer->Add(fhPtHbpXECharged_Cone2) ;
1570   
1571   outputContainer->Add(fhZTCharged) ;
1572   outputContainer->Add(fhZTPosCharged) ;
1573   outputContainer->Add(fhZTNegCharged) ;
1574   outputContainer->Add(fhPtHbpZTCharged) ;
1575   
1576   outputContainer->Add(fhPtTrigPout) ;
1577   outputContainer->Add(fhPtTrigCharged) ;
1578   
1579   TString right = "";
1580   if(fMakeSeveralUE) right = "Right";
1581   
1582   fhUePart  =  new TH1F("hUePart","UE particles distribution vs pt trig",
1583                         nptbins,ptmin,ptmax);
1584   fhUePart->SetYTitle("dNch");
1585   fhUePart->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1586   
1587   fhDeltaPhiUeChargedPt  = new TH2F
1588   (Form("hDeltaPhiUe%sChargedPt",right.Data()),"#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}}",
1589    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1590   fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi (rad)");
1591   fhDeltaPhiUeChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1592   
1593   fhXEUeCharged  =
1594   new TH2F(Form("hXEUeCharged%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event",
1595            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1596   fhXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
1597   fhXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1598   
1599   fhPtHbpXEUeCharged  =
1600   new TH2F(Form("hHbpXEUeCharged%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for Underlying Event",
1601            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1602   fhPtHbpXEUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1603   fhPtHbpXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1604   
1605   fhZTUeCharged  =
1606   new TH2F(Form("hZTUeCharged%s",right.Data()),"#it{z}_{T} for Underlying Event",
1607            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1608   fhZTUeCharged->SetYTitle("#it{z}_{T}");
1609   fhZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1610   
1611   fhPtHbpZTUeCharged  =
1612   new TH2F(Form("hHbpZTUeCharged%s",right.Data()),"#xi = ln(1/#it{z}_{T}) for Underlying Event",
1613            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1614   fhPtHbpZTUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1615   fhPtHbpZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1616   
1617   outputContainer->Add(fhUePart);
1618   outputContainer->Add(fhDeltaPhiUeChargedPt) ;
1619   outputContainer->Add(fhXEUeCharged) ;
1620   outputContainer->Add(fhPtHbpXEUeCharged) ;
1621   outputContainer->Add(fhZTUeCharged) ;
1622   outputContainer->Add(fhPtHbpZTUeCharged) ;
1623   
1624   if(fMakeSeveralUE)
1625   {
1626     fhDeltaPhiUeLeftCharged  = new TH2F
1627     ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left side range of trigger particles",
1628      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1629     fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi (rad)");
1630     fhDeltaPhiUeLeftCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1631     outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
1632     
1633     fhDeltaPhiUeLeftUpCharged  = new TH2F
1634     ("hDeltaPhiUeLeftUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Up side range of trigger particles",
1635      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1636     fhDeltaPhiUeLeftUpCharged->SetYTitle("#Delta #phi (rad)");
1637     fhDeltaPhiUeLeftUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1638     outputContainer->Add(fhDeltaPhiUeLeftUpCharged) ;
1639     
1640     fhDeltaPhiUeRightUpCharged  = new TH2F
1641     ("hDeltaPhiUeRightUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Up side range of trigger particles",
1642      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1643     fhDeltaPhiUeRightUpCharged->SetYTitle("#Delta #phi (rad)");
1644     fhDeltaPhiUeRightUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1645     outputContainer->Add(fhDeltaPhiUeRightUpCharged) ;
1646     
1647     fhDeltaPhiUeLeftDownCharged  = new TH2F
1648     ("hDeltaPhiUeLeftDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Down side range of trigger particles",
1649      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1650     fhDeltaPhiUeLeftDownCharged->SetYTitle("#Delta #phi (rad)");
1651     fhDeltaPhiUeLeftDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1652     outputContainer->Add(fhDeltaPhiUeLeftDownCharged) ;
1653     
1654     fhDeltaPhiUeRightDownCharged  = new TH2F
1655     ("hDeltaPhiUeRightDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Down side range of trigger particles",
1656      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1657     fhDeltaPhiUeRightDownCharged->SetYTitle("#Delta #phi (rad)");
1658     fhDeltaPhiUeRightDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1659     outputContainer->Add(fhDeltaPhiUeRightDownCharged) ;
1660     
1661     fhXEUeLeftCharged  =
1662     new TH2F("hXEUeChargedLeft","#it{x}_{#it{E}} with UE left side of trigger",
1663              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1664     fhXEUeLeftCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1665     fhXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1666     outputContainer->Add(fhXEUeLeftCharged) ;
1667     
1668     fhXEUeLeftUpCharged  =
1669     new TH2F("hXEUeChargedLeftUp","#it{x}_{#it{E}} with UE left Up side of trigger",
1670              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1671     fhXEUeLeftUpCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1672     fhXEUeLeftUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1673     outputContainer->Add(fhXEUeLeftUpCharged) ;
1674     
1675     fhXEUeRightUpCharged  =
1676     new TH2F("hXEUeChargedRightUp","#it{x}_{#it{E} h^{#pm}} with UE right Up side of trigger",
1677              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1678     fhXEUeRightUpCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1679     fhXEUeRightUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1680     outputContainer->Add(fhXEUeRightUpCharged) ;
1681     
1682     fhXEUeLeftDownCharged  =
1683     new TH2F("hXEUeChargedLeftDown","#it{x}_{#it{E}} with UE left Down side of trigger",
1684              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1685     fhXEUeLeftDownCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1686     fhXEUeLeftDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1687     outputContainer->Add(fhXEUeLeftDownCharged) ;
1688     
1689     fhXEUeRightDownCharged  =
1690     new TH2F("hXEUeChargedRightDown","#it{x}_{#it{E} h^{#pm}} with UE right Down side of trigger",
1691              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1692     fhXEUeRightDownCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1693     fhXEUeRightDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1694     outputContainer->Add(fhXEUeRightDownCharged) ;
1695     
1696     fhPtHbpXEUeLeftCharged  =
1697     new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/#it{x}_{#it{E}}) with charged UE left side of trigger",
1698              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1699     fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1700     fhPtHbpXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1701     outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
1702     
1703     fhZTUeLeftCharged  =
1704     new TH2F("hZTUeChargedLeft","#it{z}_{trigger h^{#pm}} = #it{p}_{T Ueh^{#pm}} / #it{p}_{T trigger} with UE left side of trigger",
1705              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1706     fhZTUeLeftCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1707     fhZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1708     outputContainer->Add(fhZTUeLeftCharged) ;
1709     
1710     fhPtHbpZTUeLeftCharged  =
1711     new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/#it{z}_{T}) with charged UE left side of trigger",
1712              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1713     fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/#it{z}_{T})");
1714     fhPtHbpZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1715     outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
1716   }
1717   
1718   if(IsPileUpAnalysisOn())
1719   {
1720     fhDeltaPhiChargedOtherBC  = new TH2F
1721     ("hDeltaPhiChargedOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC!=0",
1722      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1723     fhDeltaPhiChargedOtherBC->SetYTitle("#Delta #phi (rad)");
1724     fhDeltaPhiChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1725     
1726     fhDeltaPhiChargedPtA3GeVOtherBC  = new TH2F
1727     ("hDeltaPhiChargedPtA3GeVOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC!=0",
1728      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1729     fhDeltaPhiChargedPtA3GeVOtherBC->SetYTitle("#Delta #phi (rad)");
1730     fhDeltaPhiChargedPtA3GeVOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1731     
1732     fhPtTrigChargedOtherBC  =
1733     new TH2F("hPtTrigChargedOtherBC","trigger and charged tracks pt distribution, track BC!=0",
1734              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1735     fhPtTrigChargedOtherBC->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1736     fhPtTrigChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1737     
1738     fhXEChargedOtherBC  =
1739     new TH2F("hXEChargedOtherBC","#it{x}_{#it{E}} for charged tracks, track BC!=0",
1740              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1741     fhXEChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
1742     fhXEChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1743     
1744     fhXEUeChargedOtherBC  =
1745     new TH2F("hXEUeChargedOtherBC","#it{x}_{#it{E}} for Underlying Event, track BC!=0",
1746              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1747     fhXEUeChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
1748     fhXEUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1749     
1750     fhZTChargedOtherBC  =
1751     new TH2F("hZTChargedOtherBC","#it{z}_{T} for charged tracks, track BC!=0",
1752              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1753     fhZTChargedOtherBC->SetYTitle("#it{z}_{T}");
1754     fhZTChargedOtherBC->SetXTitle("#it{p}_{T trigger}");
1755     
1756     fhZTUeChargedOtherBC  =
1757     new TH2F("hZTUeChargedOtherBC","#it{z}_{T} for Underlying Event, track BC!=0",
1758              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1759     fhZTUeChargedOtherBC->SetYTitle("#it{z}_{T}");
1760     fhZTUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1761     
1762     outputContainer->Add(fhDeltaPhiChargedOtherBC) ;
1763     outputContainer->Add(fhDeltaPhiChargedPtA3GeVOtherBC) ;
1764     outputContainer->Add(fhXEChargedOtherBC) ;
1765     outputContainer->Add(fhXEUeChargedOtherBC) ;
1766     outputContainer->Add(fhZTChargedOtherBC) ;
1767     outputContainer->Add(fhZTUeChargedOtherBC) ;
1768     outputContainer->Add(fhPtTrigChargedOtherBC) ;
1769     
1770     fhDeltaPhiChargedBC0  = new TH2F
1771     ("hDeltaPhiChargedBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
1772      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1773     fhDeltaPhiChargedBC0->SetYTitle("#Delta #phi (rad)");
1774     fhDeltaPhiChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1775     
1776     fhDeltaPhiChargedPtA3GeVBC0  = new TH2F
1777     ("hDeltaPhiChargedPtA3GeVBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
1778      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1779     fhDeltaPhiChargedPtA3GeVBC0->SetYTitle("#Delta #phi (rad)");
1780     fhDeltaPhiChargedPtA3GeVBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1781     
1782     fhPtTrigChargedBC0  =
1783     new TH2F("hPtTrigChargedBC0","trigger and charged tracks pt distribution, track BC==0",
1784              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1785     fhPtTrigChargedBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1786     fhPtTrigChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1787     
1788     fhXEChargedBC0  =
1789     new TH2F("hXEChargedBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
1790              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1791     fhXEChargedBC0->SetYTitle("#it{x}_{#it{E}}");
1792     fhXEChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1793     
1794     fhXEUeChargedBC0  =
1795     new TH2F("hXEUeChargedBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
1796              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1797     fhXEUeChargedBC0->SetYTitle("#it{x}_{#it{E}}");
1798     fhXEUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1799     
1800     fhZTChargedBC0  =
1801     new TH2F("hZTChargedBC0","#it{z}_{T} for charged tracks, track BC==0",
1802              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1803     fhZTChargedBC0->SetYTitle("#it{z}_{T}");
1804     fhZTChargedBC0->SetXTitle("#it{p}_{T trigger}");
1805     
1806     fhZTUeChargedBC0  =
1807     new TH2F("hZTUeChargedBC0","#it{z}_{T} for Underlying Event, track BC==0",
1808              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1809     fhZTUeChargedBC0->SetYTitle("#it{z}_{T}");
1810     fhZTUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1811     
1812     outputContainer->Add(fhDeltaPhiChargedBC0) ;
1813     outputContainer->Add(fhDeltaPhiChargedPtA3GeVBC0) ;
1814     outputContainer->Add(fhXEChargedBC0) ;
1815     outputContainer->Add(fhXEUeChargedBC0) ;
1816     outputContainer->Add(fhZTChargedBC0) ;
1817     outputContainer->Add(fhZTUeChargedBC0) ;
1818     outputContainer->Add(fhPtTrigChargedBC0) ;
1819     
1820     fhPtTriggerVtxBC0  = new TH1F("hPtTriggerVtxBC0","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax);
1821     fhPtTriggerVtxBC0->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1822     
1823     fhDeltaPhiChargedVtxBC0  = new TH2F
1824     ("hDeltaPhiChargedVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
1825      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1826     fhDeltaPhiChargedVtxBC0->SetYTitle("#Delta #phi (rad)");
1827     fhDeltaPhiChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1828     
1829     fhDeltaPhiChargedPtA3GeVVtxBC0  = new TH2F
1830     ("hDeltaPhiChargedPtA3GeVVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
1831      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1832     fhDeltaPhiChargedPtA3GeVVtxBC0->SetYTitle("#Delta #phi (rad)");
1833     fhDeltaPhiChargedPtA3GeVVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1834     
1835     fhPtTrigChargedVtxBC0  =
1836     new TH2F("hPtTrigChargedVtxBC0","trigger and charged tracks pt distribution, track BC==0",
1837              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1838     fhPtTrigChargedVtxBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1839     fhPtTrigChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1840     
1841     fhXEChargedVtxBC0  =
1842     new TH2F("hXEChargedVtxBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
1843              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1844     fhXEChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
1845     fhXEChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1846     
1847     fhXEUeChargedVtxBC0  =
1848     new TH2F("hXEUeChargedVtxBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
1849              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1850     fhXEUeChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
1851     fhXEUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1852     
1853     fhZTChargedVtxBC0  =
1854     new TH2F("hZTChargedVtxBC0","#it{z}_{T} for charged tracks, track BC==0",
1855              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1856     fhZTChargedVtxBC0->SetYTitle("#it{z}_{T}");
1857     fhZTChargedVtxBC0->SetXTitle("#it{p}_{T trigger}");
1858     
1859     fhZTUeChargedVtxBC0  =
1860     new TH2F("hZTUeChargedVtxBC0","#it{z}_{T} for Underlying Event, track BC==0",
1861              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1862     fhZTUeChargedVtxBC0->SetYTitle("#it{z}_{T}");
1863     fhZTUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1864     
1865     outputContainer->Add(fhPtTriggerVtxBC0);
1866     outputContainer->Add(fhDeltaPhiChargedVtxBC0) ;
1867     outputContainer->Add(fhDeltaPhiChargedPtA3GeVVtxBC0) ;
1868     outputContainer->Add(fhXEChargedVtxBC0) ;
1869     outputContainer->Add(fhXEUeChargedVtxBC0) ;
1870     outputContainer->Add(fhZTChargedVtxBC0) ;
1871     outputContainer->Add(fhZTUeChargedVtxBC0) ;
1872     outputContainer->Add(fhPtTrigChargedVtxBC0) ;
1873     
1874     for(Int_t i = 0 ; i < 7 ; i++)
1875     {
1876       fhPtTriggerPileUp[i]  = new TH1F(Form("hPtTriggerPileUp%s",pileUpName[i].Data()),
1877                                        Form("#it{p}_{T} distribution of trigger particles, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1878       fhPtTriggerPileUp[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1879       outputContainer->Add(fhPtTriggerPileUp[i]);
1880       
1881       fhDeltaPhiChargedPileUp[i]  = new TH2F(Form("hDeltaPhiChargedPileUp%s",pileUpName[i].Data()),
1882                                              Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
1883                                              nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1884       fhDeltaPhiChargedPileUp[i]->SetYTitle("#Delta #phi (rad)");
1885       fhDeltaPhiChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1886       outputContainer->Add(fhDeltaPhiChargedPileUp[i]) ;
1887       
1888       fhDeltaPhiChargedPtA3GeVPileUp[i]  = new TH2F(Form("hDeltaPhiChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
1889                                                     Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
1890                                                     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1891       fhDeltaPhiChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #phi (rad)");
1892       fhDeltaPhiChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1893       outputContainer->Add(fhDeltaPhiChargedPtA3GeVPileUp[i]) ;
1894       
1895       fhDeltaEtaChargedPileUp[i]  = new TH2F(Form("hDeltaEtaChargedPileUp%s",pileUpName[i].Data()),
1896                                              Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
1897                                              nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1898       fhDeltaEtaChargedPileUp[i]->SetYTitle("#Delta #eta");
1899       fhDeltaEtaChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1900       outputContainer->Add(fhDeltaEtaChargedPileUp[i]) ;
1901       
1902       fhDeltaEtaChargedPtA3GeVPileUp[i]  = new TH2F(Form("hDeltaEtaChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
1903                                                     Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
1904                                                     nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1905       fhDeltaEtaChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #eta");
1906       fhDeltaEtaChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1907       outputContainer->Add(fhDeltaEtaChargedPtA3GeVPileUp[i]) ;
1908       
1909       fhXEChargedPileUp[i]  = new TH2F(Form("hXEChargedPileUp%s",pileUpName[i].Data()),
1910                                        Form("#it{x}_{#it{E}} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
1911                                        nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1912       fhXEChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
1913       fhXEChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1914       outputContainer->Add(fhXEChargedPileUp[i]) ;
1915       
1916       fhXEUeChargedPileUp[i]  = new TH2F(Form("hXEUeChargedPileUp%s",pileUpName[i].Data()),
1917                                          Form("#it{x}_{#it{E}} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
1918                                          nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1919       fhXEUeChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
1920       fhXEUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1921       outputContainer->Add(fhXEUeChargedPileUp[i]) ;
1922       
1923       fhZTChargedPileUp[i]  = new TH2F(Form("hZTChargedPileUp%s",pileUpName[i].Data()),
1924                                        Form("#it{z}_{T} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
1925                                        nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1926       fhZTChargedPileUp[i]->SetYTitle("#it{z}_{T}");
1927       fhZTChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1928       outputContainer->Add(fhZTChargedPileUp[i]) ;
1929       
1930       fhZTUeChargedPileUp[i]  = new TH2F(Form("hZTUeChargedPileUp%s",pileUpName[i].Data()),
1931                                          Form("#it{z}_{T} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
1932                                          nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1933       fhZTUeChargedPileUp[i]->SetYTitle("#it{z}_{T}");
1934       fhZTUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1935       outputContainer->Add(fhZTUeChargedPileUp[i]) ;
1936       
1937       fhPtTrigChargedPileUp[i]  = new TH2F(Form("hPtTrigChargedPileUp%s",pileUpName[i].Data()),
1938                                            Form("trigger and charged tracks pt distribution, %s Pile-Up event",pileUpName[i].Data()),
1939                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1940       fhPtTrigChargedPileUp[i]->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1941       fhPtTrigChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1942       outputContainer->Add(fhPtTrigChargedPileUp[i]) ;
1943       
1944     }
1945   }
1946
1947   if(IsHighMultiplicityAnalysisOn())
1948   {
1949     Int_t nMultiBins = GetNCentrBin();
1950     fhDeltaPhiChargedMult = new TH2F*[nMultiBins] ;
1951     fhDeltaEtaChargedMult = new TH2F*[nMultiBins] ;
1952     fhXEMult              = new TH2F*[nMultiBins] ;
1953     fhXEUeMult            = new TH2F*[nMultiBins] ;
1954     fhZTMult              = new TH2F*[nMultiBins] ;
1955     fhZTUeMult            = new TH2F*[nMultiBins] ;
1956     
1957     for(Int_t im=0; im<nMultiBins; im++)
1958     {
1959       fhDeltaPhiChargedMult[im]  = new TH2F
1960       (Form("hDeltaPhiCharged_Mult%d",im),Form("#Delta #phi charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1961       fhDeltaPhiChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1962       fhDeltaPhiChargedMult[im]->SetYTitle("#Delta #phi (rad)");
1963       
1964       fhDeltaEtaChargedMult[im]  = new TH2F
1965       (Form("hDeltaEtaCharged_Mult%d",im),Form("#Delta #eta charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
1966       fhDeltaEtaChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1967       fhDeltaEtaChargedMult[im]->SetYTitle("#Delta #eta");
1968       
1969       fhXEMult[im]  = new TH2F
1970       (Form("hXECharged_Mult%d",im),Form("#it{x}_{E} charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1971       fhXEMult[im]->SetYTitle("#it{x}_{E}");
1972       fhXEMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1973       
1974       fhXEUeMult[im]  = new TH2F
1975       (Form("hXEUeCharged_Mult%d",im),Form("#it{x}_{E} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1976       fhXEUeMult[im]->SetYTitle("#it{x}_{E}");
1977       fhXEUeMult[im]->SetXTitle("#it{p}_{T trigger}(GeV/#it{c})");
1978       
1979       fhZTMult[im]  = new TH2F
1980       (Form("hZTCharged_Mult%d",im),Form("#it{z}_{T} charged  Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1981       fhZTMult[im]->SetYTitle("#it{z}_{T}");
1982       fhZTMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1983       
1984       fhZTUeMult[im]  = new TH2F
1985       (Form("hZTUeCharged_Mult%d",im),Form("#it{z}_{T} UE charged  Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1986       fhZTUeMult[im]->SetYTitle("#it{z}_{T}");
1987       fhZTUeMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1988       
1989       outputContainer->Add(fhDeltaPhiChargedMult[im]) ;
1990       outputContainer->Add(fhDeltaEtaChargedMult[im]) ;
1991       outputContainer->Add(fhXEMult  [im]);
1992       outputContainer->Add(fhXEUeMult[im]);
1993       outputContainer->Add(fhZTMult  [im]);
1994       outputContainer->Add(fhZTUeMult[im]);
1995     }
1996   }
1997   
1998   if(fFillBradHisto)
1999   {
2000     fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger #it{p}_{T} vs associated hadron #it{p}_{T} from background",
2001                                    nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
2002     fhAssocPtBkg->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2003     fhAssocPtBkg->SetYTitle("#it{p}_{T associated} (GeV/#it{c})");
2004     outputContainer->Add(fhAssocPtBkg) ;
2005     
2006     fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} ",
2007                               nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2008     fhDeltaPhiBrad->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2009     fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2010     outputContainer->Add(fhDeltaPhiBrad) ;
2011   }
2012   
2013   fhDeltaPhiDeltaEtaAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2014   fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
2015   fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
2016   fhDeltaPhiAssocPtBinDEta0  = new TH2F*[fNAssocPtBins*nz];
2017   if(fFillMomImbalancePtAssocBinsHisto)
2018   {
2019     fhXEAssocPtBin           = new TH2F*[fNAssocPtBins*nz];
2020     fhZTAssocPtBin           = new TH2F*[fNAssocPtBins*nz];
2021   }
2022
2023   if(fCorrelVzBin)
2024   {
2025     fhXEVZ = new TH2F*[nz];
2026     fhZTVZ = new TH2F*[nz];
2027   }
2028   
2029   if(fFillBradHisto)
2030     fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2031   
2032
2033   fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
2034   if(fFillEtaGapsHisto)fhDeltaPhiAssocPtBinDEta08       = new TH2F*[fNAssocPtBins*nz];
2035   if(fDecayTrigger)    fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2036   
2037   
2038   if(fHMPIDCorrelation)
2039   {
2040     fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins*nz];
2041     fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins*nz];
2042   }
2043
2044   for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2045   {
2046     for(Int_t z = 0 ; z < nz ; z++)
2047     {
2048       Int_t bin = i*nz+z;
2049       
2050       if(fCorrelVzBin)
2051       {
2052         sz = Form("_vz%d",z);
2053         tz = Form(", #it{v}_{#it{z}} bin %d",z);
2054       }
2055       
2056       //printf("iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2057       
2058       fhDeltaPhiDeltaEtaAssocPtBin[bin]  = new TH2F(Form("hDeltaPhiDeltaEtaPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2059                                                     Form("#Delta #phi vs #Delta #eta for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2060                                                     ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
2061       fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2062       fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetYTitle("#Delta #eta");
2063       
2064       fhDeltaPhiAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2065                                            Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2066                                            nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2067       fhDeltaPhiAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2068       fhDeltaPhiAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2069       
2070       outputContainer->Add(fhDeltaPhiDeltaEtaAssocPtBin[bin]) ;
2071       outputContainer->Add(fhDeltaPhiAssocPtBin[bin]) ;
2072       
2073       if(fFillEtaGapsHisto)
2074       {
2075         fhDeltaPhiAssocPtBinDEta08[bin] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2076                                                    Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2077                                                    nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2078         fhDeltaPhiAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2079         fhDeltaPhiAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2080         
2081         fhDeltaPhiAssocPtBinDEta0[bin] = new TH2F(Form("hDeltaPhiDeltaEta0PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2082                                                   Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0.", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2083                                                   nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2084         fhDeltaPhiAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2085         fhDeltaPhiAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2086         
2087         outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[bin]) ;
2088         outputContainer->Add(fhDeltaPhiAssocPtBinDEta0[bin]) ;
2089       }
2090       
2091       if(fDecayTrigger)
2092       {
2093         fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s_bit%d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data(),fDecayBits[0]),
2094                                                          Form("#Delta #phi vs #it{p}_{T trigger} tagged as decay for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, Bit %d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data(),fDecayBits[0]),
2095                                                          nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2096         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2097         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2098         
2099         outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[bin]) ;
2100       }
2101       
2102       if(fFillBradHisto)
2103       {
2104         fhDeltaPhiBradAssocPtBin[bin] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2105                                                  Form("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2106                                                  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2107         fhDeltaPhiBradAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2108         fhDeltaPhiBradAssocPtBin[bin]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2109         outputContainer->Add(fhDeltaPhiBradAssocPtBin[bin]) ;
2110       }
2111       
2112       if(fHMPIDCorrelation)
2113       {
2114         fhDeltaPhiAssocPtBinHMPID[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2115                                                   Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track having HMPID signal", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2116                                                   nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2117         fhDeltaPhiAssocPtBinHMPID[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})" );
2118         fhDeltaPhiAssocPtBinHMPID[bin]->SetYTitle("#Delta #phi (rad)");
2119         
2120         fhDeltaPhiAssocPtBinHMPIDAcc[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2121                                                      Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track within 5<phi<20 deg", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2122                                                      nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2123         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2124         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetYTitle("#Delta #phi (rad)");
2125         
2126         outputContainer->Add(fhDeltaPhiAssocPtBinHMPID   [bin]) ;
2127         outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[bin]) ;
2128       }
2129     }
2130   }
2131   
2132   if(fFillMomImbalancePtAssocBinsHisto)
2133   {
2134     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2135     {
2136       fhXEAssocPtBin[i]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2137                                          Form("#it{x}_{#it{E}} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2138                                          nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2139       fhXEAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2140       fhXEAssocPtBin[i]->SetYTitle("#it{x}_{#it{E}}");
2141       
2142       fhZTAssocPtBin[i]       = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2143                                          Form("#it{z}_{T} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2144                                          nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2145       fhZTAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2146       fhZTAssocPtBin[i]->SetYTitle("#it{z}_{T}");
2147       
2148       
2149       outputContainer->Add(fhXEAssocPtBin[i]);
2150       outputContainer->Add(fhZTAssocPtBin[i]);
2151     }
2152   }
2153
2154   if(fCorrelVzBin)
2155   {
2156     for(Int_t z = 0 ; z < nz ; z++)
2157     {
2158       sz = Form("_vz%d",z);
2159       tz = Form(", #it{v}_{#it{z}} bin %d",z);
2160       
2161       fhXEVZ[z]       = new TH2F(Form("hXE%s", sz.Data()),
2162                                  Form("#it{x}_{#it{E}} vs #it{p}_{T trigger}%s", tz.Data()),
2163                                  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2164       fhXEVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2165       fhXEVZ[z]->SetYTitle("#it{x}_{#it{E}}");
2166       
2167       fhZTVZ[z]       = new TH2F(Form("hZT%s",sz.Data()),
2168                                  Form("#it{z}_{T} vs #it{p}_{T trigger}%s", tz.Data()),
2169                                  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2170       fhZTVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2171       fhZTVZ[z]->SetYTitle("#it{z}_{T}");
2172       
2173       outputContainer->Add(fhXEVZ[z]);
2174       outputContainer->Add(fhZTVZ[z]);
2175     }
2176   }
2177
2178  
2179   if(fPi0Trigger)
2180   {
2181     fhPtPi0DecayRatio  = new TH2F
2182     ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
2183      nptbins,ptmin,ptmax, 100,0.,2.);
2184     fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
2185     fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
2186     outputContainer->Add(fhPtPi0DecayRatio) ;
2187     
2188     fhDeltaPhiPi0DecayCharged  = new TH2F
2189     ("hDeltaPhiPi0DecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
2190      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2191     fhDeltaPhiPi0DecayCharged->SetYTitle("#Delta #phi (rad)");
2192     fhDeltaPhiPi0DecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2193     
2194     fhXEPi0DecayCharged  =
2195     new TH2F("hXEPi0DecayCharged","#it{x}_{#it{E}}  Decay",
2196              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2197     fhXEPi0DecayCharged->SetYTitle("#it{x}_{#it{E}}");
2198     fhXEPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2199     
2200     fhZTPi0DecayCharged  =
2201     new TH2F("hZTPi0DecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
2202              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2203     fhZTPi0DecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
2204     fhZTPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2205     
2206     outputContainer->Add(fhDeltaPhiPi0DecayCharged) ;
2207     outputContainer->Add(fhXEPi0DecayCharged) ;
2208     outputContainer->Add(fhZTPi0DecayCharged) ;
2209   }
2210   
2211   if(fDecayTrigger)
2212   {
2213     for(Int_t ibit = 0; ibit< fNDecayBits; ibit++)
2214     {
2215       fhDeltaPhiDecayCharged[ibit]  = new TH2F
2216       (Form("hDeltaPhiDecayCharged_bit%d",fDecayBits[ibit]),
2217        Form("#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2218        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2219       fhDeltaPhiDecayCharged[ibit]->SetYTitle("#Delta #phi (rad)");
2220       fhDeltaPhiDecayCharged[ibit]->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2221       
2222       fhXEDecayCharged[ibit]  =
2223       new TH2F(Form("hXEDecayCharged_bit%d",fDecayBits[ibit]),
2224                Form("#it{x}_{#it{E}}  Decay, Bit %d",fDecayBits[ibit]),
2225                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2226       fhXEDecayCharged[ibit]->SetYTitle("#it{x}_{#it{E}}");
2227       fhXEDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2228       
2229       fhZTDecayCharged[ibit]  =
2230       new TH2F(Form("hZTDecayCharged_bit%d",fDecayBits[ibit]),
2231                Form("#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2232                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2233       fhZTDecayCharged[ibit]->SetYTitle("#it{z}_{decay h^{#pm}}");
2234       fhZTDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2235       
2236       outputContainer->Add(fhDeltaPhiDecayCharged[ibit]) ;
2237       outputContainer->Add(fhXEDecayCharged[ibit]) ;
2238       outputContainer->Add(fhZTDecayCharged[ibit]) ;
2239     }
2240   }
2241   
2242   //Correlation with neutral hadrons
2243   if(fNeutralCorr)
2244   {
2245     fhDeltaPhiDeltaEtaNeutral  = new TH2F
2246     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
2247      ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
2248     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi (rad)");
2249     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
2250           
2251     fhPhiNeutral  = new TH2F
2252     ("hPhiNeutral","#phi_{#pi^{0}}  vs #it{p}_{T #pi^{0}}",
2253      nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2254     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
2255     fhPhiNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2256     
2257     fhEtaNeutral  = new TH2F
2258     ("hEtaNeutral","#eta_{#pi^{0}}  vs #it{p}_{T #pi^{0}}",
2259      nptbins,ptmin,ptmax,200,-1.,1.);
2260     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
2261     fhEtaNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2262     
2263     fhDeltaPhiNeutral  = new TH2F
2264     ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T trigger}",
2265      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2266     fhDeltaPhiNeutral->SetYTitle("#Delta #phi (rad)");
2267     fhDeltaPhiNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2268     
2269     fhDeltaPhiNeutralPt  = new TH2F
2270     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2271      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2272     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi (rad)");
2273     fhDeltaPhiNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2274     
2275     fhDeltaEtaNeutral  = new TH2F
2276     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs #it{p}_{T trigger}",
2277      nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2278     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
2279     fhDeltaEtaNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2280     
2281     fhXENeutral  =
2282     new TH2F("hXENeutral","#it{x}_{#it{E}} for #pi^{0} associated",
2283              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2284     fhXENeutral->SetYTitle("#it{x}_{#it{E}}");
2285     fhXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2286     
2287     fhPtHbpXENeutral  =
2288     new TH2F("hHbpXENeutral","#xi = ln(1/#it{x}_{#it{E}})for #pi^{0} associated",
2289              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2290     fhPtHbpXENeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2291     fhPtHbpXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2292     
2293     fhZTNeutral  =
2294     new TH2F("hZTNeutral","#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2295              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2296     fhZTNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2297     fhZTNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2298     
2299     fhPtHbpZTNeutral  =
2300     new TH2F("hHbpZTNeutral","#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2301              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2302     fhPtHbpZTNeutral->SetYTitle("ln(1/#it{z}_{T})");
2303     fhPtHbpZTNeutral->SetXTitle("#it{p}_{T trigger}");
2304     
2305     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
2306     outputContainer->Add(fhPhiNeutral) ;
2307     outputContainer->Add(fhEtaNeutral) ;
2308     outputContainer->Add(fhDeltaPhiNeutral) ;
2309     outputContainer->Add(fhDeltaPhiNeutralPt) ;
2310     outputContainer->Add(fhDeltaEtaNeutral) ;
2311     outputContainer->Add(fhXENeutral) ;
2312     outputContainer->Add(fhPtHbpXENeutral) ;
2313     outputContainer->Add(fhZTNeutral) ;
2314     outputContainer->Add(fhPtHbpZTNeutral) ;
2315     
2316     fhDeltaPhiUeNeutralPt  = new TH2F
2317     (Form("hDeltaPhiUe%sNeutralPt",right.Data()),"#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2318      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2319     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi (rad)");
2320     fhDeltaPhiUeNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2321     
2322     fhXEUeNeutral  =
2323     new TH2F(Form("hXEUeNeutral%s",right.Data()),"#it{x}_{#it{E}} for #pi^{0} associated",
2324              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2325     fhXEUeNeutral->SetYTitle("#it{x}_{#it{E}}");
2326     fhXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2327     
2328     fhPtHbpXEUeNeutral  =
2329     new TH2F(Form("hHbpXEUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2330              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2331     fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2332     fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2333     
2334     fhZTUeNeutral  =
2335     new TH2F(Form("hZTUeNeutral%s",right.Data()),"#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2336              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2337     fhZTUeNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2338     fhZTUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2339     
2340     fhPtHbpZTUeNeutral  =
2341     new TH2F(Form("hHbpZTUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2342              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2343     fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{z}_{T})");
2344     fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2345
2346     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
2347     outputContainer->Add(fhXEUeNeutral) ;
2348     outputContainer->Add(fhPtHbpXEUeNeutral) ;
2349     outputContainer->Add(fhZTUeNeutral) ;
2350     outputContainer->Add(fhPtHbpZTUeNeutral) ;
2351
2352     if(fMakeSeveralUE)
2353     {
2354       fhDeltaPhiUeLeftNeutral  = new TH2F
2355       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs #it{p}_{T h^{0}} with neutral UE left side range of trigger particles",
2356        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2357       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi (rad)");
2358       fhDeltaPhiUeLeftNeutral->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2359       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
2360       
2361       fhXEUeLeftNeutral  =
2362       new TH2F("hXEUeNeutralLeft","#it{x}_{#it{E}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2363                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2364       fhXEUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2365       fhXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2366       outputContainer->Add(fhXEUeLeftNeutral) ;
2367       
2368       fhPtHbpXEUeLeftNeutral  =
2369       new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/#it{x}_{#it{E}}) with neutral UE left side of trigger",
2370                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2371       fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2372       fhPtHbpXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2373       outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
2374       
2375       fhZTUeLeftNeutral  =
2376       new TH2F("hZTUeNeutralLeft","#it{z}_{trigger h^{0}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2377                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2378       fhZTUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2379       fhZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2380       outputContainer->Add(fhZTUeLeftNeutral) ;
2381       
2382       fhPtHbpZTUeLeftNeutral  =
2383       new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/#it{z}_{T}) with neutral UE left side of trigger",
2384                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2385       fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/#it{z}_{T})");
2386       fhPtHbpZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger}");
2387       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
2388     }
2389     
2390     if(fPi0Trigger)
2391     {
2392       fhDeltaPhiPi0DecayNeutral  = new TH2F
2393       ("hDeltaPhiPi0DecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
2394        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2395       fhDeltaPhiPi0DecayNeutral->SetYTitle("#Delta #phi (rad)");
2396       fhDeltaPhiPi0DecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2397       
2398       fhXEPi0DecayNeutral  =
2399       new TH2F("hXEPi0DecayNeutral","#it{x}_{#it{E}} for decay trigger",
2400                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2401       fhXEPi0DecayNeutral->SetYTitle("#it{x}_{#it{E}}");
2402       fhXEPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2403       
2404       fhZTPi0DecayNeutral  =
2405       new TH2F("hZTPi0DecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
2406                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2407       fhZTPi0DecayNeutral->SetYTitle("#it{z}_{h^{0}}");
2408       fhZTPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2409       
2410       outputContainer->Add(fhDeltaPhiPi0DecayNeutral) ;
2411       outputContainer->Add(fhXEPi0DecayNeutral) ;
2412       outputContainer->Add(fhZTPi0DecayNeutral) ;
2413     }
2414   }//Correlation with neutral hadrons
2415   
2416   // If data is MC, fill more histograms, depending on origin
2417   if(IsDataMC())
2418   {
2419     for(Int_t i= fMCGenTypeMin; i <= fMCGenTypeMax; i++)
2420     {
2421       fhMCPtTrigger[i]  = new TH1F (Form("hMCPtTrigger_%s",nameMC[i].Data()),
2422                                  Form("MC %s: generated trigger #it{p}_{T}",nameMC[i].Data()),
2423                                  nptbins,ptmin,ptmax);
2424       fhMCPtTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2425       
2426       fhMCPhiTrigger[i]  = new TH2F (Form("hMCPhiTrigger_%s",nameMC[i].Data()),
2427                                      Form("MC %s: generated trigger #phi",nameMC[i].Data()),
2428                                      nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2429       fhMCPhiTrigger[i]->SetYTitle("#phi (rad)");
2430       fhMCPhiTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2431       
2432       fhMCEtaTrigger[i]  = new TH2F (Form("hMCEtaTrigger_%s",nameMC[i].Data()),
2433                                      Form("MC %s: generated trigger #eta",nameMC[i].Data()),
2434                                      nptbins,ptmin,ptmax, netabins,etamin,etamax);
2435       fhMCEtaTrigger[i]->SetYTitle("#eta");
2436       fhMCEtaTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2437       
2438       if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2439       {
2440         fhMCPtTriggerNotLeading[i]  = new TH1F (Form("hMCPtTriggerNotLeading_%s",nameMC[i].Data()),
2441                                                 Form("MC %s: generated trigger #it{p}_{T}, when not leading of primaries",nameMC[i].Data()),
2442                                                 nptbins,ptmin,ptmax);
2443         fhMCPtTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2444         
2445         fhMCPhiTriggerNotLeading[i]  = new TH2F (Form("hMCPhiTriggerNotLeading_%s",nameMC[i].Data()),
2446                                                  Form("MC %s: generated trigger #phi, when not leading of primaries",nameMC[i].Data()),
2447                                                  nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2448         fhMCPhiTriggerNotLeading[i]->SetYTitle("#phi (rad)");
2449         fhMCPhiTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2450         
2451         
2452         fhMCEtaTriggerNotLeading[i]  = new TH2F (Form("hMCEtaTriggerNotLeading_%s",nameMC[i].Data()),
2453                                                  Form("MC %s: generated triogger #eta, when not leading of primaries",nameMC[i].Data()),
2454                                                  nptbins,ptmin,ptmax, netabins,etamin,etamax);
2455         fhMCEtaTriggerNotLeading[i]->SetYTitle("#eta ");
2456         fhMCEtaTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2457       }
2458       
2459       fhMCEtaCharged[i]  = new TH2F (Form("hMCEtaCharged_%s",nameMC[i].Data()),
2460                                      Form("MC %s: #eta_{h^{#pm}}  vs #it{p}_{T #pm}",nameMC[i].Data()),
2461                                      nptbins,ptmin,ptmax,100,-1.,1.);
2462       fhMCEtaCharged[i]->SetYTitle("#eta_{h^{#pm}} (rad)");
2463       fhMCEtaCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2464       
2465       fhMCPhiCharged[i]  = new TH2F(Form("hMCPhiCharged_%s",nameMC[i].Data()),
2466                                     Form("MC %s: phi_{h^{#pm}}  vs #it{p}_{T #pm}",nameMC[i].Data()),
2467                                     nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2468       fhMCPhiCharged[i]->SetYTitle("MC #phi_{h^{#pm}} (rad)");
2469       fhMCPhiCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2470       
2471       fhMCDeltaPhiDeltaEtaCharged[i]  = new TH2F (Form("hMCDeltaPhiDeltaEtaCharged_%s",nameMC[i].Data()),
2472                                                   Form("MC %s: phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",nameMC[i].Data()),
2473                                                   ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2474       fhMCDeltaPhiDeltaEtaCharged[i]->SetXTitle("#Delta #phi (rad)");
2475       fhMCDeltaPhiDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2476       
2477       fhMCDeltaEtaCharged[i]  = new TH2F (Form("hMCDeltaEtaCharged_%s",nameMC[i].Data()),
2478                                           Form("MC %s: #eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger} and #it{p}_{T assoc}",nameMC[i].Data()),
2479                                           nptbins,ptmin,ptmax,ndeltaetabins ,deltaetamin,deltaetamax);
2480       fhMCDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2481       fhMCDeltaEtaCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2482       
2483       fhMCDeltaPhiCharged[i]  = new TH2F  (Form("hMCDeltaPhiCharged_%s",nameMC[i].Data()),
2484                                            Form("MC %s: #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",nameMC[i].Data()),
2485                                            nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2486       fhMCDeltaPhiCharged[i]->SetYTitle("#Delta #phi (rad)");
2487       fhMCDeltaPhiCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2488
2489       fhMCDeltaPhiChargedPt[i]  = new TH2F (Form("hMCDeltaPhiChargedPt_%s",nameMC[i].Data()),
2490                                             Form("MC %s: #phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",nameMC[i].Data()),
2491                                             nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2492       fhMCDeltaPhiChargedPt[i]->SetYTitle("#Delta #phi (rad)");
2493       fhMCDeltaPhiChargedPt[i]->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2494       
2495       fhMCPtXECharged[i]  = new TH2F (Form("hMCPtXECharged_%s",nameMC[i].Data()),
2496                                       Form("MC %s: #it{x}_{#it{E}} with charged hadrons",nameMC[i].Data()),
2497                                       nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2498       fhMCPtXECharged[i]->SetYTitle("#it{x}_{#it{E}}");
2499       fhMCPtXECharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2500       
2501       fhMCPtHbpXECharged[i]  = new TH2F(Form("hMCHbpXECharged_%s",nameMC[i].Data()),
2502                                         Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",nameMC[i].Data()),
2503                                         nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2504       fhMCPtHbpXECharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2505       fhMCPtHbpXECharged[i]->SetXTitle("#it{p}_{T trigger}");
2506       
2507       fhMCPtZTCharged[i]  = new TH2F(Form("hMCPtZTCharged_%s",nameMC[i].Data()),
2508                                      Form("MC %s: #it{z}_{T} with charged hadrons",nameMC[i].Data()),
2509                                      nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2510       fhMCPtZTCharged[i]->SetYTitle("#it{z}_{T}");
2511       fhMCPtZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2512       
2513       fhMCPtHbpZTCharged[i]  = new TH2F(Form("hMCHbpZTCharged_%s",nameMC[i].Data()),
2514                                         Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons",nameMC[i].Data()),
2515                                              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2516       fhMCPtHbpZTCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2517       fhMCPtHbpZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2518       
2519       fhMCPtTrigPout[i]  = new TH2F(Form("hMCPtTrigPout_%s",nameMC[i].Data()),
2520                                     Form("MC %s: #it{p}_{out} with triggers",nameMC[i].Data()),
2521                                     nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
2522       fhMCPtTrigPout[i]->SetYTitle("#it{p}_{out} (GeV/#it{c})");
2523       fhMCPtTrigPout[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2524       
2525       fhMCPtAssocDeltaPhi[i]  = new TH2F(Form("hMCPtAssocDeltaPhi_%s",nameMC[i].Data()),
2526                                          Form("MC %s: #Delta #phi with associated charged hadrons",nameMC[i].Data()),
2527                                          nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2528       fhMCPtAssocDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
2529       fhMCPtAssocDeltaPhi[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2530       
2531       outputContainer->Add(fhMCPtTrigger[i]);
2532       outputContainer->Add(fhMCPhiTrigger[i]);
2533       outputContainer->Add(fhMCEtaTrigger[i]);
2534       
2535       if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2536       {
2537         outputContainer->Add(fhMCPtTriggerNotLeading[i]);
2538         outputContainer->Add(fhMCPhiTriggerNotLeading[i]);
2539         outputContainer->Add(fhMCEtaTriggerNotLeading[i]);
2540       }
2541       
2542       outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged[i]);
2543       outputContainer->Add(fhMCPhiCharged[i]) ;
2544       outputContainer->Add(fhMCEtaCharged[i]) ;
2545       outputContainer->Add(fhMCDeltaEtaCharged[i]) ;
2546       outputContainer->Add(fhMCDeltaPhiCharged[i]) ;
2547       
2548       outputContainer->Add(fhMCDeltaPhiChargedPt[i]) ;
2549       outputContainer->Add(fhMCPtXECharged[i]) ;
2550       outputContainer->Add(fhMCPtZTCharged[i]) ;
2551       outputContainer->Add(fhMCPtHbpXECharged[i]) ;
2552       outputContainer->Add(fhMCPtHbpZTCharged[i]) ;
2553       outputContainer->Add(fhMCPtTrigPout[i]) ;
2554       outputContainer->Add(fhMCPtAssocDeltaPhi[i]) ;
2555
2556       // Underlying event
2557       
2558       fhMCUePart[i]  =
2559       new TH1F(Form("hMCUePart_%s",nameMC[i].Data()),
2560                Form("MC %s: UE particles distribution vs #it{p}_{T trigger}",nameMC[i].Data()),
2561                nptbins,ptmin,ptmax);
2562       fhMCUePart[i]->SetYTitle("#it{dN}^{ch}");
2563       fhMCUePart[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2564       
2565       fhMCPtXEUeCharged[i]  =
2566       new TH2F(Form("hMCPtXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2567                Form("MC %s: #it{x}_{#it{E}} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2568                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2569       fhMCPtXEUeCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2570       fhMCPtXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2571       
2572       fhMCPtHbpXEUeCharged[i] =
2573       new TH2F(Form("hMCPtHbpXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2574                Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2575                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2576       fhMCPtHbpXEUeCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2577       fhMCPtHbpXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2578       
2579       fhMCPtZTUeCharged[i] =
2580       new TH2F(Form("hMCPtZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2581                Form("MC %s: #it{z}_{T} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2582                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2583       fhMCPtZTUeCharged[i]->SetYTitle("#it{z}_{T}");
2584       fhMCPtZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2585       
2586       fhMCPtHbpZTUeCharged[i] =
2587       new TH2F(Form("hMCPtHbpZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2588                Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2589                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2590       fhMCPtHbpZTUeCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2591       fhMCPtHbpZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2592       
2593       outputContainer->Add(fhMCUePart[i]);
2594       outputContainer->Add(fhMCPtXEUeCharged[i]) ;
2595       outputContainer->Add(fhMCPtZTUeCharged[i]) ;
2596       outputContainer->Add(fhMCPtHbpZTUeCharged[i]);
2597       outputContainer->Add(fhMCPtHbpXEUeCharged[i]);
2598
2599       if(fMakeSeveralUE)
2600       {
2601         fhMCPtXEUeLeftCharged[i]  = new TH2F(Form("hMCPtXEUeChargedLeft_%s",nameMC[i].Data()),
2602                                              Form("MC %s: #it{x}_{#it{E}} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2603                                              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2604         fhMCPtXEUeLeftCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2605         fhMCPtXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2606         
2607         fhMCPtHbpXEUeLeftCharged[i] = new TH2F(Form("hMCPtHbpXEUeChargedLeft_%s",nameMC[i].Data()),
2608                                                Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2609                                                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2610         fhMCPtHbpXEUeLeftCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2611         fhMCPtHbpXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2612         
2613         fhMCPtZTUeLeftCharged[i]  = new TH2F(Form("hMCPtZTUeChargedLeft_%s",nameMC[i].Data()),
2614                                              Form("MC %s: #it{z}_{T} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2615                                              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2616         fhMCPtZTUeLeftCharged[i]->SetYTitle("#it{z}_{T}");
2617         fhMCPtZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2618         
2619         fhMCPtHbpZTUeLeftCharged[i] = new TH2F(Form("hMCPtHbpZTUeChargedLeft_%s",nameMC[i].Data()),
2620                                                Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2621                                                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2622         fhMCPtHbpZTUeLeftCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2623         fhMCPtHbpZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2624         
2625         outputContainer->Add(fhMCPtXEUeLeftCharged[i]) ;
2626         outputContainer->Add(fhMCPtZTUeLeftCharged[i]) ;
2627         outputContainer->Add(fhMCPtHbpXEUeLeftCharged[i]);
2628         outputContainer->Add(fhMCPtHbpZTUeLeftCharged[i]) ;
2629         
2630       }
2631     }
2632   } //for MC histogram
2633   
2634   if(DoOwnMix())
2635   {
2636     //create event containers
2637     
2638     if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
2639     {
2640       Int_t nvz = GetNZvertBin();
2641       Int_t nrp = GetNRPBin();
2642       Int_t nce = GetNCentrBin();
2643       
2644       fListMixTrackEvents= new TList*[nvz*nrp*nce] ;
2645       
2646       for( Int_t ice = 0 ; ice < nce ; ice++ )
2647       {
2648         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
2649         {
2650           for( Int_t irp = 0 ; irp < nrp ; irp++ )
2651           {
2652             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
2653             
2654             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
2655             //       ic,iz, irp, bin);
2656             
2657             fListMixTrackEvents[bin] = new TList() ;
2658             fListMixTrackEvents[bin]->SetOwner(kFALSE);
2659           }
2660         }
2661       }
2662     }
2663     
2664     fhPtTriggerMixed  = new TH1F ("hPtTriggerMixed","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax);
2665     fhPtTriggerMixed->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2666     
2667     if(fCorrelVzBin)
2668     {
2669       fhPtTriggerMixedVzBin  = new TH2F ("hPtTriggerMixedVzBin","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
2670       fhPtTriggerMixedVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2671       fhPtTriggerMixedVzBin->SetYTitle("#it{v}_{#it{z}} bin");
2672       outputContainer->Add(fhPtTriggerMixedVzBin);
2673     }
2674     
2675     fhPtTriggerMixedBin  = new TH2F ("hPtTriggerMixedBin","#it{p}_{T} distribution of trigger particles vs mixing bin", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
2676     fhPtTriggerMixedBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2677     fhPtTriggerMixedBin->SetYTitle("Bin");
2678     
2679     fhPhiTriggerMixed  = new TH2F ("hPhiTriggerMixed","#phi distribution of trigger Particles, used for mixing",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2680     fhPhiTriggerMixed->SetYTitle("#phi (rad)");
2681     
2682     fhEtaTriggerMixed  = new TH2F ("hEtaTriggerMixed","#eta distribution of trigger, used for mixing",nptbins,ptmin,ptmax, netabins,etamin,etamax);
2683     fhEtaTriggerMixed->SetYTitle("#eta ");
2684     
2685     outputContainer->Add(fhPtTriggerMixed);
2686     outputContainer->Add(fhPtTriggerMixedBin);
2687     outputContainer->Add(fhPhiTriggerMixed);
2688     outputContainer->Add(fhEtaTriggerMixed);
2689     
2690     // Fill the cluster pool only in isolation analysis or if requested
2691     if( neutralMix && (!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForCaloExists())))
2692     {
2693       Int_t nvz = GetNZvertBin();
2694       Int_t nrp = GetNRPBin();
2695       Int_t nce = GetNCentrBin();
2696       
2697       fListMixCaloEvents= new TList*[nvz*nrp*nce] ;
2698       
2699       for( Int_t ice = 0 ; ice < nce ; ice++ )
2700       {
2701         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
2702         {
2703           for( Int_t irp = 0 ; irp < nrp ; irp++ )
2704           {
2705             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
2706             
2707             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
2708             //       ic,iz, irp, bin);
2709             
2710             fListMixCaloEvents[bin] = new TList() ;
2711             fListMixCaloEvents[bin]->SetOwner(kFALSE);
2712           }
2713         }
2714       }
2715     }
2716     
2717     //Init the list in the reader if not done previously
2718     if(fUseMixStoredInReader)
2719     {
2720       if( !GetReader()->ListWithMixedEventsForTracksExists() )
2721         GetReader()->SetListWithMixedEventsForTracks(fListMixTrackEvents);
2722       
2723       if( !GetReader()->ListWithMixedEventsForCaloExists()   )
2724         GetReader()->SetListWithMixedEventsForCalo  (fListMixCaloEvents );
2725     }
2726     
2727     fhEventBin=new TH1I("hEventBin","Number of triggers per bin(cen,vz,rp)",
2728                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2729                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2730     fhEventBin->SetXTitle("event bin");
2731     outputContainer->Add(fhEventBin) ;
2732     
2733     fhEventMixBin=new TH1I("hEventMixBin","Number of triggers mixed per event bin(cen,vz,rp)",
2734                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2735                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2736     fhEventMixBin->SetXTitle("event bin");
2737     outputContainer->Add(fhEventMixBin) ;
2738
2739     fhEventMBBin=new TH1I("hEventMBBin","Number of min bias events per bin(cen,vz,rp)",
2740                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2741                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2742     fhEventMBBin->SetXTitle("event bin");
2743     outputContainer->Add(fhEventMBBin) ;
2744
2745     fhNtracksMB=new TH2F("hNtracksMBEvent","Number of filtered tracks in MB event per event bin",ntrbins,trmin,trmax,
2746                          GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2747                          GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2748     fhNtracksMB->SetYTitle("event bin");
2749     fhNtracksMB->SetXTitle("#it{N}_{track}");
2750     outputContainer->Add(fhNtracksMB);
2751
2752     if( neutralMix )
2753     {
2754       fhNclustersMB=new TH2F("hNclustersMBEvent","Number of filtered clusters in MB events per event bin",nclbins,clmin,clmax,
2755                              GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2756                              GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2757       fhNclustersMB->SetYTitle("event bin");
2758       fhNclustersMB->SetXTitle("#it{N}_{cluster}");
2759       outputContainer->Add(fhNclustersMB);
2760     }
2761     
2762     fhMixDeltaPhiCharged  = new TH2F
2763     ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
2764      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2765     fhMixDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
2766     fhMixDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2767     outputContainer->Add(fhMixDeltaPhiCharged);
2768     
2769     fhMixDeltaPhiDeltaEtaCharged  = new TH2F
2770     ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
2771      ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2772     fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
2773     fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
2774     outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
2775     
2776     fhMixXECharged  =
2777     new TH2F("hMixXECharged","Mixed event : #it{x}_{#it{E}} for charged tracks",
2778              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2779     fhMixXECharged->SetYTitle("#it{x}_{#it{E}}");
2780     fhMixXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2781     outputContainer->Add(fhMixXECharged);
2782     
2783     fhMixXEUeCharged  =
2784     new TH2F("hMixXEUeCharged","Mixed event : #it{x}_{#it{E}} for charged tracks in Ue region",
2785              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2786     fhMixXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
2787     fhMixXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2788     outputContainer->Add(fhMixXEUeCharged);
2789     
2790     fhMixHbpXECharged  =
2791     new TH2F("hMixHbpXECharged","mixed event : #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
2792              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2793     fhMixHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2794     fhMixHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2795     outputContainer->Add(fhMixHbpXECharged);
2796     
2797     fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins*nz];
2798     fhMixDeltaPhiChargedAssocPtBinDEta08   = new TH2F*[fNAssocPtBins*nz];
2799     fhMixDeltaPhiChargedAssocPtBinDEta0    = new TH2F*[fNAssocPtBins*nz];
2800     fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2801     
2802     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2803     {
2804       for(Int_t z = 0 ; z < nz ; z++)
2805       {
2806         Int_t bin = i*nz+z;
2807         
2808         if(fCorrelVzBin)
2809         {
2810           sz = Form("_vz%d",z);
2811           tz = Form(", #it{v}_{#it{z}} bin %d",z);
2812         }
2813         
2814         //printf("MIX : iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2815         
2816         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2817                                                                Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2818                                                                ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2819         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2820         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetYTitle("#Delta #eta");
2821         
2822         outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]);
2823         
2824         fhMixDeltaPhiChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2825                                                        Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2826                                                        nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2827         fhMixDeltaPhiChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2828         fhMixDeltaPhiChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2829         
2830         outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[bin]);
2831         
2832         if(fFillEtaGapsHisto)
2833         {
2834           fhMixDeltaPhiChargedAssocPtBinDEta08[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2835                                                                Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2836                                                                nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2837           fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2838           fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2839           
2840           fhMixDeltaPhiChargedAssocPtBinDEta0[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2841                                                               Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2842                                                               nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2843           fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2844           fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2845           
2846           outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[bin]);
2847           outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta0[bin]);
2848         }
2849       }
2850     }
2851   }
2852   
2853   return outputContainer;
2854   
2855 }
2856
2857 //_____________________________________________________________________________________________________________________
2858 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(Int_t indexPhoton1, Int_t indexPhoton2, Int_t idetector)
2859 {
2860   // Get the momentum of the pi0/eta assigned decay photons
2861   // In case of pi0/eta trigger, we may want to check their decay correlation,
2862   // get their decay children
2863   
2864   if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
2865   
2866   if(GetDebug() > 1)
2867     printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
2868   
2869   TObjArray * clusters  = 0x0 ;
2870   if(idetector==kEMCAL) clusters = GetEMCALClusters() ;
2871   else                  clusters = GetPHOSClusters()  ;
2872   
2873   for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
2874   {
2875     AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
2876     
2877     if(photon->GetID()==indexPhoton1) photon->GetMomentum(fDecayMom1,GetVertex(0)) ;
2878     if(photon->GetID()==indexPhoton2) photon->GetMomentum(fDecayMom2,GetVertex(0)) ;
2879     
2880     if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", fDecayMom1.Pt(), fDecayMom2.Pt());
2881     
2882   } //cluster loop        
2883   
2884   return kTRUE;
2885   
2886
2887
2888 //_____________________________________________________________
2889 Int_t AliAnaParticleHadronCorrelation::GetMCTagHistogramIndex(Int_t mcTag)
2890 {
2891   // Index of MC histograms depending on MC origin
2892   
2893   if     ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||        
2894            GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) return 0;
2895   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))           return 1;    
2896   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      return 2;
2897   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  ||
2898            GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      return 3;
2899   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    return 4;
2900   else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))      return 5;
2901   else                                                                                    return 6;
2902   
2903 }
2904
2905 //_________________________________________
2906 void AliAnaParticleHadronCorrelation::Init()
2907 {
2908   //Init
2909   //Do some checks
2910   
2911   if(!GetReader()->IsCTSSwitchedOn())
2912     AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
2913 }
2914
2915 //____________________________________________________
2916 void AliAnaParticleHadronCorrelation::InitParameters()
2917 {
2918   
2919   //Initialize the parameters of the analysis.
2920   SetInputAODName("Particle");
2921   SetAODObjArrayName("Hadrons");  
2922   AddToHistogramsName("AnaHadronCorr_");
2923   
2924   SetPtCutRange(0.,300);
2925   fDeltaPhiMinCut       = TMath::DegToRad()*120.;
2926   fDeltaPhiMaxCut       = TMath::DegToRad()*240. ;
2927   fSelectIsolated       = kFALSE;
2928   fMakeSeveralUE        = kFALSE;
2929   fUeDeltaPhiMinCut     = TMath::DegToRad()*60.;
2930   fUeDeltaPhiMaxCut     = TMath::DegToRad()*120 ;
2931   
2932   fNeutralCorr          = kFALSE ;
2933   fPi0Trigger           = kFALSE ;
2934   fDecayTrigger         = kFALSE ;
2935   fHMPIDCorrelation     = kFALSE ;
2936   
2937   fMakeAbsoluteLeading  = kTRUE;
2938   fMakeNearSideLeading  = kFALSE;
2939
2940   fNAssocPtBins         = 9   ;
2941   fAssocPtBinLimit[0]   = 0.2 ; 
2942   fAssocPtBinLimit[1]   = 0.5 ; 
2943   fAssocPtBinLimit[2]   = 1.0 ; 
2944   fAssocPtBinLimit[3]   = 2.0 ; 
2945   fAssocPtBinLimit[4]   = 3.0 ; 
2946   fAssocPtBinLimit[5]   = 4.0 ; 
2947   fAssocPtBinLimit[6]   = 5.0 ;
2948   fAssocPtBinLimit[7]   = 6.0 ;
2949   fAssocPtBinLimit[8]   = 7.0 ;
2950   fAssocPtBinLimit[9]   = 8.0 ;
2951   fAssocPtBinLimit[10]  = 9.0 ; 
2952   fAssocPtBinLimit[11]  = 10.0 ; 
2953   fAssocPtBinLimit[12]  = 12.0 ; 
2954   fAssocPtBinLimit[13]  = 14.0 ; 
2955   fAssocPtBinLimit[14]  = 16.0 ; 
2956   fAssocPtBinLimit[15]  = 20.0 ; 
2957   fAssocPtBinLimit[16]  = 30.0 ;
2958   fAssocPtBinLimit[17]  = 40.0 ;
2959   fAssocPtBinLimit[18]  = 50.0 ;
2960   fAssocPtBinLimit[19]  = 200.0 ;
2961   
2962   fUseMixStoredInReader = kTRUE;
2963   
2964   fM02MinCut   = -1 ;
2965   fM02MaxCut   = -1 ;
2966   
2967   fSelectLeadingHadronAngle = kFALSE;
2968   fFillLeadHadOppositeHisto = kFALSE;
2969   fMinLeadHadPhi = 150*TMath::DegToRad();
2970   fMaxLeadHadPhi = 210*TMath::DegToRad();
2971
2972   fMinLeadHadPt  = 1;
2973   fMaxLeadHadPt  = 100;
2974   
2975   fMCGenTypeMin = 0;
2976   fMCGenTypeMax = 6;
2977   
2978   fNDecayBits = 1;
2979   fDecayBits[0] = AliNeutralMesonSelection::kPi0;
2980   fDecayBits[1] = AliNeutralMesonSelection::kEta;
2981   fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
2982   fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
2983 }
2984
2985 //_________________________________________________________________________
2986 Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
2987 {
2988   // Check if the what of the selected triggers is leading particle comparing
2989   // with all the triggers, all the tracks or all the clusters (if requested for the clusters).
2990   
2991   Double_t ptTrig      = GetMinPt();
2992   Double_t phiTrig     = 0 ;
2993   fLeadingTriggerIndex =-1 ;
2994   Int_t index          =-1 ;
2995   AliAODPWG4ParticleCorrelation* pLeading = 0;
2996
2997   // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
2998   
2999   for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3000   {
3001     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3002     particle->SetLeadingParticle(kFALSE); // set it later
3003
3004     // Vertex cut in case of mixing
3005     Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3006     if(check ==  0) continue;
3007     if(check == -1) return kFALSE; // not sure if it is correct.
3008     
3009     // find the leading particles with highest momentum
3010     if (particle->Pt() > ptTrig)
3011     {
3012       ptTrig   = particle->Pt() ;
3013       phiTrig  = particle->Phi();
3014       index    = iaod     ;
3015       pLeading = particle ;
3016     }
3017   }// finish search of leading trigger particle on the AOD branch.
3018   
3019   if(index < 0) return kFALSE;
3020   
3021   //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3022   
3023   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3024   
3025   // Compare if it is the leading of all tracks
3026   
3027   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3028   {
3029     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3030     
3031     if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
3032        track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3)   ) continue ;
3033     
3034     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3035     Float_t pt   = fTrackVector.Pt();
3036     Float_t phi  = fTrackVector.Phi() ;
3037     if(phi < 0) phi+=TMath::TwoPi();
3038     
3039     //jump out this event if near side associated particle pt larger than trigger
3040     if (fMakeNearSideLeading)
3041     {
3042       Float_t deltaPhi = phiTrig-phi;
3043       if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3044       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3045       
3046       if(pt > ptTrig && deltaPhi < TMath::PiOver2())  return kFALSE;
3047     }
3048     //jump out this event if there is any other particle with pt larger than trigger
3049     else
3050     {
3051       if(pt > ptTrig)  return kFALSE ;
3052     }
3053    }// track loop
3054
3055   // Compare if it is leading of all calorimeter clusters
3056   
3057   if(fCheckLeadingWithNeutralClusters)
3058   {
3059     // Select the calorimeter cluster list
3060     TObjArray * nePl = 0x0;
3061     if      (pLeading->GetDetectorTag() == kPHOS )
3062       nePl = GetPHOSClusters();
3063     else
3064       nePl = GetEMCALClusters();
3065     
3066     if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3067     
3068     for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3069     {
3070       AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3071       
3072       if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3073       
3074       cluster->GetMomentum(fMomentum,GetVertex(0));
3075       
3076       Float_t pt   = fMomentum.Pt();
3077       Float_t phi  = fMomentum.Phi() ;
3078       if(phi < 0) phi+=TMath::TwoPi();
3079       
3080       if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3081       
3082       //jump out this event if near side associated particle pt larger than trigger
3083       // not really needed for calorimeter, unless DCal is included
3084       if (fMakeNearSideLeading)
3085       {
3086         Float_t deltaPhi = phiTrig-phi;
3087         if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3088         if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3089         
3090         if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3091       }
3092       //jump out this event if there is any other particle with pt larger than trigger
3093       else
3094       {
3095         if(pt > ptTrig)  return kFALSE ;
3096       }
3097     }// cluster loop
3098   } // check neutral clusters
3099   
3100   fLeadingTriggerIndex = index ;
3101   pLeading->SetLeadingParticle(kTRUE);
3102
3103   if( GetDebug() > 1 ) printf("\t particle AOD with index %d is leading with pT %2.2f\n", fLeadingTriggerIndex, pLeading->Pt());
3104   
3105   return kTRUE;
3106   
3107 }
3108
3109 //_________________________________________________________________
3110 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
3111 {  
3112   //Particle-Hadron Correlation Analysis, fill histograms
3113   
3114   if(!GetInputAODBranch())
3115   {
3116     AliFatal(Form("No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data()));
3117     return ; // coverity
3118   }
3119   
3120   Int_t naod = GetInputAODBranch()->GetEntriesFast();
3121   if( naod == 0 )
3122   {
3123     if(GetDebug() > 1)
3124       printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No particle AOD found! \n");
3125     
3126     return ; // no trigger particles found.
3127   }
3128
3129   if(GetDebug() > 1)
3130   {
3131     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
3132     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", naod);
3133     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In CTS aod entries %d\n",   GetCTSTracks()->GetEntriesFast());
3134   }
3135   
3136   //------------------------------------------------------
3137   // Find leading trigger if analysis request only leading,
3138   // if there is no leading trigger, then skip the event
3139   
3140   Int_t iaod = 0 ;
3141   if( fMakeAbsoluteLeading || fMakeNearSideLeading )
3142   {
3143     Bool_t leading = IsTriggerTheEventLeadingParticle();
3144     
3145     if(GetDebug() > 1)
3146       printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - AOD Leading trigger? %d, with index %d\n",leading,fLeadingTriggerIndex);
3147     
3148     if(!leading)
3149     {
3150       if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Leading was requested and not found\n");
3151       return ;
3152     }
3153     else
3154     {
3155       // Select only the leading in the trigger AOD loop
3156       naod = fLeadingTriggerIndex+1 ;
3157       iaod = fLeadingTriggerIndex   ;
3158     }
3159   }
3160
3161   //------------------------------------------------------
3162   // Get event multiplicity and bins
3163   
3164   Float_t cen         = GetEventCentrality();
3165   Float_t ep          = GetEventPlaneAngle();
3166   if(IsHighMultiplicityAnalysisOn()) fhTriggerEventPlaneCentrality->Fill(cen,ep);
3167
3168   Int_t   mixEventBin = GetEventMixBin();
3169   Int_t   vzbin       = GetEventVzBin();
3170
3171   //------------------------------------------------------
3172   // Loop on trigger AOD
3173   
3174   for( iaod = 0; iaod < naod; iaod++ )
3175   {
3176     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3177     
3178     //
3179     // Trigger particle selection criteria:
3180     //
3181     Float_t pt = particle->Pt();
3182     
3183     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3184
3185     fhPtTriggerInput->Fill(pt);
3186     
3187     //
3188     // check if it was a calorimeter cluster
3189     // and if the shower shape cut was requested apply it.
3190     // Not needed if already done at the particle identification level,
3191     // but for isolation studies, it is preferred not to remove so we do it here
3192     //
3193     Int_t clID1  = particle->GetCaloLabel(0) ;
3194     Int_t clID2  = particle->GetCaloLabel(1) ; // for photon clusters should not be set.
3195     if( GetDebug() > 1 ) printf("%s Trigger : id1 %d, id2 %d, min %f, max %f, det %d\n",
3196            GetInputAODName().Data(),clID1,clID2,fM02MinCut,fM02MaxCut,particle->GetDetectorTag());
3197     
3198     if(clID1 > 0 && clID2 < 0 && fM02MaxCut > 0 && fM02MinCut > 0)
3199     {
3200 //      Int_t iclus = -1;
3201 //      TObjArray* clusters = 0x0;
3202 //      if     (particle->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
3203 //      else if(particle->GetDetectorTag() == kPHOS ) clusters = GetPHOSClusters();
3204 //      
3205 //      if(clusters)
3206 //      {
3207 //        AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
3208 //        Float_t m02 = cluster->GetM02();
3209 //        if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3210 //      }
3211       
3212       Float_t m02 = particle->GetM02();
3213       if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3214       
3215       fhPtTriggerSSCut->Fill(pt);
3216     }
3217     
3218     //
3219     // Check if the particle is isolated or if we want to take the isolation into account
3220     // This bool is set in AliAnaParticleIsolation
3221     //
3222     if(OnlyIsolated())
3223     {
3224       if( !particle->IsIsolated() ) continue;
3225       fhPtTriggerIsoCut->Fill(pt);
3226     }
3227     
3228     //
3229     // Check if trigger is in fiducial region
3230     //
3231     if(IsFiducialCutOn())
3232     {
3233       Bool_t in = GetFiducialCut()->IsInFiducialCut(particle->Eta(),particle->Phi(),particle->GetDetectorTag()) ;
3234       if(! in ) continue ;
3235     }
3236     
3237     fhPtTriggerFidCut->Fill(pt);
3238     
3239     //---------------------------------------
3240     // Make correlation
3241     
3242     // Find the leading hadron in the opposite hemisphere to the triggeer
3243     // and accept the trigger if leading is in defined window.
3244     Bool_t okLeadHad = kTRUE;
3245     if(fSelectLeadingHadronAngle || fFillLeadHadOppositeHisto)
3246     {
3247       okLeadHad = FindLeadingOppositeHadronInWindow(particle);
3248       if(!okLeadHad && fSelectLeadingHadronAngle) continue;
3249     }
3250     
3251     //
3252     // Charged particles correlation
3253     //
3254     MakeChargedCorrelation(particle);
3255     
3256     // MC
3257     Int_t mcIndex = -1;
3258     Int_t mcTag   = particle->GetTag();
3259     Bool_t lostDecayPair = kFALSE;
3260     if(IsDataMC())
3261     {
3262       mcIndex = GetMCTagHistogramIndex(mcTag);
3263       lostDecayPair = GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost);
3264       MakeMCChargedCorrelation(particle->GetLabel(), mcIndex,lostDecayPair);
3265     }
3266     
3267     // Do own mixed event with charged,
3268     // add event and remove previous or fill the mixed histograms
3269     if(DoOwnMix())
3270       MakeChargedMixCorrelation(particle);
3271
3272     //
3273     // Neutral particles correlation
3274     //
3275     if(fNeutralCorr)
3276       MakeNeutralCorrelation(particle);
3277     
3278     //----------------------------------------------------------------
3279     // Fill trigger pT related histograms if not absolute leading
3280     
3281     //
3282     // pT of the trigger, vs trigger origin if MC
3283     //
3284     fhPtTrigger->Fill(pt);
3285     if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3286     {
3287       fhPtTriggerMC[mcIndex]->Fill(pt);
3288       if( lostDecayPair && mcIndex==2 )
3289         fhPtTriggerMC[7]->Fill(pt);
3290     }
3291     
3292     if(fDecayTrigger)
3293     {
3294       Int_t decayTag = particle->DecayTag();
3295       if(decayTag < 0) decayTag = 0;
3296       
3297       for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
3298       {
3299         if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3300         {
3301           fhPtDecayTrigger[ibit]->Fill(pt);
3302           
3303           if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3304           {
3305             fhPtDecayTriggerMC[ibit][mcIndex]->Fill(pt);
3306             if(lostDecayPair && mcInd