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