AliPhysics  095eea3 (095eea3)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaParticleHadronCorrelation.cxx
Go to the documentation of this file.
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 // --- ROOT system ---
17 #include <TClass.h>
18 #include <TMath.h>
19 #include <TH2F.h>
20 #include <TDatabasePDG.h>
21 #include "TParticle.h"
22 #include <TClonesArray.h>
23 #include <TList.h>
24 #include <TObjString.h>
25 
26 //---- ANALYSIS system ----
27 #include "AliNeutralMesonSelection.h"
28 #include "AliAnaParticleHadronCorrelation.h"
29 #include "AliCaloTrackReader.h"
30 #include "AliAODPWG4ParticleCorrelation.h"
31 #include "AliFiducialCut.h"
32 #include "AliVTrack.h"
33 #include "AliVCluster.h"
34 #include "AliMCAnalysisUtils.h"
35 #include "TParticle.h"
36 #include "AliStack.h"
37 #include "AliAODMCParticle.h"
38 #include "AliMixedEvent.h"
39 #include "AliAnalysisManager.h"
40 #include "AliInputEventHandler.h"
41 #include "AliEventplane.h"
42 #include "AliIsolationCut.h"
43 #include "AliESDEvent.h"
44 #include "AliAODEvent.h"
45 
46 // --- Detectors ---
47 #include "AliEMCALGeometry.h"
48 #include "AliPHOSGeoUtils.h"
49 
51 ClassImp(AliAnaParticleHadronCorrelation) ;
53 
54 
55 //___________________________________________________________________
57 //___________________________________________________________________
58 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
59 AliAnaCaloTrackCorrBaseClass(),
60 fFillAODWithReferences(0), fCheckLeadingWithNeutralClusters(0),
61 fMaxAssocPt(1000.), fMinAssocPt(0.),
62 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
63 fSelectIsolated(0), fMakeSeveralUE(0),
64 fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.),
65 fPi0AODBranchName(""), fAODNamepTInConeHisto(""), fNeutralCorr(0),
66 fPi0Trigger(0), fDecayTrigger(0),
67 fNDecayBits(0), fDecayBits(),
68 fNBkgBin(0),
69 fMakeAbsoluteLeading(0), fMakeNearSideLeading(0),
70 fLeadingTriggerIndex(-1), fHMPIDCorrelation(0), fFillBradHisto(0),
71 fNAssocPtBins(0), fAssocPtBinLimit(),
72 fCorrelVzBin(0),
73 fListMixTrackEvents(), fListMixCaloEvents(),
74 fUseMixStoredInReader(0), fFillNeutralEventMixPool(0),
75 fM02MaxCut(0), fM02MinCut(0),
76 fSelectLeadingHadronAngle(0), fFillLeadHadOppositeHisto(0),
77 fMinLeadHadPhi(0), fMaxLeadHadPhi(0),
78 fMinLeadHadPt(0), fMaxLeadHadPt(0),
79 fFillEtaGapsHisto(1), fFillMomImbalancePtAssocBinsHisto(0),
80 fFillInvMassHisto(0), fFillBkgBinsHisto(0),
81 fFillTaggedDecayHistograms(0), fDecayTagsM02Cut(0),
82 fMCGenTypeMin(0), fMCGenTypeMax(0),
83 fTrackVector(), fMomentum(), fMomentumIM(),
84 fDecayMom1(), fDecayMom2(),
85 //Histograms
86 fhPtTriggerInput(0), fhPtTriggerSSCut(0),
87 fhPtTriggerIsoCut(0), fhPtTriggerFidCut(0),
88 fhPtTrigger(0), fhPtTriggerVtxBC0(0),
89 fhPtTriggerVzBin(0), fhPtTriggerBin(0),
90 fhPhiTrigger(0), fhEtaTrigger(0),
91 fhPtTriggerMC(),
92 fhPtDecayTrigger(), fhPtDecayTriggerMC(),
93 fhPtTriggerCentrality(0), fhPtTriggerEventPlane(0),
94 fhTriggerEventPlaneCentrality(0),
95 fhPtTriggerMixed(0), fhPtTriggerMixedVzBin(0), fhPtTriggerMixedBin(0),
96 fhPhiTriggerMixed(0), fhEtaTriggerMixed(0),
97 fhPtLeadingOppositeHadron(0), fhPtDiffPhiLeadingOppositeHadron(0), fhPtDiffEtaLeadingOppositeHadron(0),
98 fhPtNoLeadingOppositeHadron(0), fhEtaPhiNoLeadingOppositeHadron(0),
99 fhDeltaPhiDeltaEtaCharged(0),
100 fhPhiCharged(0), fhEtaCharged(0),
101 fhDeltaPhiCharged(0), fhDeltaEtaCharged(0),
102 fhDeltaPhiChargedPt(0), fhDeltaPhiUeChargedPt(0),
103 fhUePart(0),
104 fhXECharged(0), fhXECharged_Cone2(0), fhXEUeCharged(0),
105 fhXEUeChargedSmallCone(0), fhXEUeChargedMediumCone(0), fhXEUeChargedLargeCone(0),
106 fhXEPosCharged(0), fhXENegCharged(0),
107 fhPtHbpXECharged(0), fhPtHbpXECharged_Cone2(0), fhPtHbpXEUeCharged(0),
108 fhZTCharged(0), fhZTUeCharged(0),
109 fhZTPosCharged(0), fhZTNegCharged(0),
110 fhPtHbpZTCharged(0), fhPtHbpZTUeCharged(0),
111 fhXEChargedMC(), fhDeltaPhiChargedMC(),
112 fhXEUeChargedRightMC(), fhXEUeChargedLeftMC(),
113 fhDeltaPhiDeltaEtaChargedPtA3GeV(0),
114 fhDeltaPhiChargedPtA3GeV(0), fhDeltaEtaChargedPtA3GeV(0),
115 // Pile-Up
116 fhDeltaPhiChargedPileUp(), fhDeltaEtaChargedPileUp(),
117 fhDeltaPhiChargedPtA3GeVPileUp(), fhDeltaEtaChargedPtA3GeVPileUp(),
118 fhXEChargedPileUp(), fhXEUeChargedPileUp(),
119 fhZTChargedPileUp(), fhZTUeChargedPileUp(),
120 fhPtTrigChargedPileUp(),
121 fhDeltaPhiChargedOtherBC(), fhDeltaPhiChargedPtA3GeVOtherBC(),
122 fhXEChargedOtherBC(), fhXEUeChargedOtherBC(),
123 fhZTChargedOtherBC(), fhZTUeChargedOtherBC(),
124 fhPtTrigChargedOtherBC(),
125 fhDeltaPhiChargedBC0(), fhDeltaPhiChargedPtA3GeVBC0(),
126 fhXEChargedBC0(), fhXEUeChargedBC0(),
127 fhZTChargedBC0(), fhZTUeChargedBC0(),
128 fhPtTrigChargedBC0(),
129 fhDeltaPhiChargedVtxBC0(), fhDeltaPhiChargedPtA3GeVVtxBC0(),
130 fhXEChargedVtxBC0(), fhXEUeChargedVtxBC0(),
131 fhZTChargedVtxBC0(), fhZTUeChargedVtxBC0(),
132 fhPtTrigChargedVtxBC0(),
133 fhDeltaPhiUeLeftCharged(0),
134 fhDeltaPhiUeLeftUpCharged(0), fhDeltaPhiUeRightUpCharged(0),
135 fhDeltaPhiUeLeftDownCharged(0), fhDeltaPhiUeRightDownCharged(0),
136 fhXEUeLeftCharged(0),
137 fhXEUeLeftUpCharged(0), fhXEUeRightUpCharged(0),
138 fhXEUeLeftDownCharged(0), fhXEUeRightDownCharged(0),
139 fhPtHbpXEUeLeftCharged(0), fhZTUeLeftCharged(0),
140 fhPtHbpZTUeLeftCharged(0),
141 fhPtTrigPout(0), fhPtTrigCharged(0),
142 fhDeltaPhiChargedMult(0x0), fhDeltaEtaChargedMult(0x0),
143 fhXEMult(0x0), fhXEUeMult(0x0),
144 fhZTMult(0x0), fhZTUeMult(0x0),
145 fhAssocPtBkg(0), fhDeltaPhiDeltaEtaAssocPtBin(0),
146 fhDeltaPhiAssocPtBin(0),
147 fhDeltaPhiAssocPtBinDEta08(0), fhDeltaPhiAssocPtBinDEta0(0),
148 fhDeltaPhiAssocPtBinHMPID(0), fhDeltaPhiAssocPtBinHMPIDAcc(0),
149 fhDeltaPhiBradAssocPtBin(0), fhDeltaPhiBrad(0),
150 fhXEAssocPtBin(0), fhZTAssocPtBin(0),
151 fhXEVZ(0), fhZTVZ(0),
152 fhDeltaPhiDeltaEtaNeutral(0),
153 fhPhiNeutral(0), fhEtaNeutral(0),
154 fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),
155 fhDeltaPhiNeutralPt(0), fhDeltaPhiUeNeutralPt(0),
156 fhXENeutral(0), fhXEUeNeutral(0),
157 fhPtHbpXENeutral(0), fhPtHbpXEUeNeutral(0),
158 fhZTNeutral(0), fhZTUeNeutral(0),
159 fhPtHbpZTNeutral(0), fhPtHbpZTUeNeutral(0),
160 fhDeltaPhiUeLeftNeutral(0), fhXEUeLeftNeutral(0),
161 fhPtHbpXEUeLeftNeutral(0), fhZTUeLeftNeutral(0),
162 fhPtHbpZTUeLeftNeutral(0), fhPtPi0DecayRatio(0),
163 fhDeltaPhiPi0DecayCharged(0), fhXEPi0DecayCharged(0), fhZTPi0DecayCharged(0),
164 fhDeltaPhiPi0DecayNeutral(0), fhXEPi0DecayNeutral(0), fhZTPi0DecayNeutral(0),
165 fhDeltaPhiDecayCharged(), fhXEDecayCharged(), fhZTDecayCharged(),
166 fhDeltaPhiDecayChargedAssocPtBin(),
167 fhMCPtTrigger(), fhMCPhiTrigger(), fhMCEtaTrigger(),
168 fhMCPtTriggerNotLeading(), fhMCPhiTriggerNotLeading(), fhMCEtaTriggerNotLeading(),
169 fhMCEtaCharged(), fhMCPhiCharged(),
170 fhMCDeltaEtaCharged(), fhMCDeltaPhiCharged(),
171 fhMCDeltaPhiDeltaEtaCharged(), fhMCDeltaPhiChargedPt(),
172 fhMCPtXECharged(), fhMCPtXEUeCharged(),
173 fhMCPtXEUeLeftCharged(),
174 fhMCPtHbpXECharged(), fhMCPtHbpXEUeCharged(),
175 fhMCPtHbpXEUeLeftCharged(),
176 fhMCUePart(),
177 fhMCPtZTCharged(), fhMCPtZTUeCharged(),
178 fhMCPtZTUeLeftCharged(),
179 fhMCPtHbpZTCharged(), fhMCPtHbpZTUeCharged(),
180 fhMCPtHbpZTUeLeftCharged(),
181 fhMCPtTrigPout(), fhMCPtAssocDeltaPhi(),
182 // Mixing
183 fhNEventsTrigger(0), fhNtracksMB(0), fhNclustersMB(0),
184 fhMixDeltaPhiCharged(0), fhMixDeltaPhiDeltaEtaCharged(0),
185 fhMixXECharged(0), fhMixXEUeCharged(0), fhMixHbpXECharged(0),
186 fhMixDeltaPhiChargedAssocPtBin(),
187 fhMixDeltaPhiChargedAssocPtBinDEta08(),
188 fhMixDeltaPhiChargedAssocPtBinDEta0(),
189 fhMixDeltaPhiDeltaEtaChargedAssocPtBin(),
190 fhEventBin(0), fhEventMixBin(0), fhEventMBBin(0),
191 fhMassPtTrigger(0), fhMCMassPtTrigger(),
192 fhPtLeadInConeBin(), fhPtSumInConeBin(),
193 fhPtLeadConeBinDecay(), fhSumPtConeBinDecay(),
194 fhPtLeadConeBinMC(), fhSumPtConeBinMC(),
195 fhTrackResolution(0), fhTrackResolutionUE(0)
196 {
197  InitParameters();
198 
199  for(Int_t i = 0; i < fgkNmcTypes; i++)
200  {
201  fhPtTriggerMC[i] = 0;
202  fhXEChargedMC[i] = 0;
203  fhDeltaPhiChargedMC[i] = 0;
204  fhMCMassPtTrigger [i] = 0 ;
205  fhXEUeChargedRightMC[i] = 0;
206  fhXEUeChargedLeftMC[i] = 0;
207 
208  for(Int_t ib = 0; ib < AliNeutralMesonSelection::fgkMaxNDecayBits; ib++) fhPtDecayTriggerMC[ib][i] = 0;
209  }
210 
211  for(Int_t ib = 0; ib < AliNeutralMesonSelection::fgkMaxNDecayBits; ib++) fhPtDecayTrigger[ib] = 0;
212 
213  for(Int_t i = 0; i < 7; i++)
214  {
215  fhPtTriggerPileUp [i] = 0 ;
216  fhDeltaPhiChargedPileUp [i] = 0 ; fhDeltaEtaChargedPileUp [i] = 0 ;
217  fhXEChargedPileUp [i] = 0 ; fhXEUeChargedPileUp [i] = 0 ;
218  fhZTChargedPileUp [i] = 0 ; fhZTUeChargedPileUp [i] = 0 ;
219  fhPtTrigChargedPileUp [i] = 0 ;
220  fhDeltaPhiChargedPtA3GeVPileUp[i] = 0 ; fhDeltaEtaChargedPtA3GeVPileUp[i] = 0 ;
221  }
222 }
223 
224 //_________________________________________________________________
226 //_________________________________________________________________
227 AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation()
228 {
229  if(DoOwnMix())
230  {
231  if(fListMixTrackEvents)
232  {
233  for(Int_t iz=0; iz < GetNZvertBin(); iz++)
234  {
235  for(Int_t ic=0; ic < GetNCentrBin(); ic++)
236  {
237  for(Int_t irp=0; irp<GetNRPBin(); irp++)
238  {
239  Int_t bin = GetEventMixBin(ic, iz, irp);
240  fListMixTrackEvents[bin]->Delete() ;
241  delete fListMixTrackEvents[bin] ;
242  }
243  }
244  }
245  }
246 
247  delete[] fListMixTrackEvents;
248 
249  if(fListMixCaloEvents)
250  {
251  for(Int_t iz=0; iz < GetNZvertBin(); iz++)
252  {
253  for(Int_t ic=0; ic < GetNCentrBin(); ic++)
254  {
255  for(Int_t irp=0; irp<GetNRPBin(); irp++)
256  {
257  Int_t bin = GetEventMixBin(ic, iz, irp);
258  fListMixCaloEvents[bin]->Delete() ;
259  delete fListMixCaloEvents[bin] ;
260  }
261  }
262  }
263  }
264 
265  delete[] fListMixCaloEvents;
266  }
267 }
268 
269 //____________________________________________________________________________________________________________________________________
271 //____________________________________________________________________________________________________________________________________
272 void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Float_t ptAssoc, Float_t ptTrig, Int_t bin,
273  Float_t phiAssoc, Float_t phiTrig, Float_t deltaPhi,
274  Float_t etaAssoc, Float_t etaTrig,
275  Int_t decayTag, Float_t hmpidSignal, Int_t outTOF,
276  Int_t cen, Int_t mcTag)
277 {
278  Float_t deltaEta = etaTrig-etaAssoc;
279  Float_t deltaPhiOrg = phiTrig-phiAssoc;
280 
281  fhEtaCharged ->Fill(ptAssoc, etaAssoc, GetEventWeight());
282  fhPhiCharged ->Fill(ptAssoc, phiAssoc, GetEventWeight());
283  fhDeltaEtaCharged ->Fill(ptTrig , deltaEta, GetEventWeight());
284  fhDeltaPhiCharged ->Fill(ptTrig , deltaPhi, GetEventWeight());
285  fhDeltaPhiChargedPt->Fill(ptAssoc, deltaPhi, GetEventWeight());
286  fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta, GetEventWeight());
287 
288  if(ptAssoc > 3 )
289  {
290  fhDeltaEtaChargedPtA3GeV ->Fill(ptTrig , deltaEta, GetEventWeight());
291  fhDeltaPhiChargedPtA3GeV ->Fill(ptTrig , deltaPhi, GetEventWeight());
292  fhDeltaPhiDeltaEtaChargedPtA3GeV->Fill(deltaPhi, deltaEta, GetEventWeight());
293  }
294 
295  // Pile up studies
296 
297  if(IsPileUpAnalysisOn())
298  {
299  if (outTOF==1)
300  {
301  fhDeltaPhiChargedOtherBC->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
302  if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVOtherBC->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
303  }
304  else if(outTOF==0)
305  {
306  fhDeltaPhiChargedBC0->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
307  if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVBC0->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
308  }
309 
310  Int_t vtxBC = GetReader()->GetVertexBC();
311  if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
312  {
313  fhDeltaPhiChargedVtxBC0->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
314  if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVVtxBC0->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
315  }
316 
317  if(GetReader()->IsPileUpFromSPD())
318  {
319  fhDeltaEtaChargedPileUp[0]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
320  fhDeltaPhiChargedPileUp[0]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
321  }
322  if(GetReader()->IsPileUpFromEMCal())
323  {
324  fhDeltaEtaChargedPileUp[1]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
325  fhDeltaPhiChargedPileUp[1]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
326  }
327  if(GetReader()->IsPileUpFromSPDOrEMCal())
328  {
329  fhDeltaEtaChargedPileUp[2]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
330  fhDeltaPhiChargedPileUp[2]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
331  }
332  if(GetReader()->IsPileUpFromSPDAndEMCal())
333  {
334  fhDeltaEtaChargedPileUp[3]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
335  fhDeltaPhiChargedPileUp[3]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
336  }
337  if(GetReader()->IsPileUpFromSPDAndNotEMCal())
338  {
339  fhDeltaEtaChargedPileUp[4]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
340  fhDeltaPhiChargedPileUp[4]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
341  }
342  if(GetReader()->IsPileUpFromEMCalAndNotSPD())
343  {
344  fhDeltaEtaChargedPileUp[5]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
345  fhDeltaPhiChargedPileUp[5]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
346  }
347  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
348  {
349  fhDeltaEtaChargedPileUp[6]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
350  fhDeltaPhiChargedPileUp[6]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
351  }
352 
353  if(ptAssoc > 3 )
354  {
355  if(GetReader()->IsPileUpFromSPD())
356  {
357  fhDeltaEtaChargedPtA3GeVPileUp[0]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
358  fhDeltaPhiChargedPtA3GeVPileUp[0]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
359  }
360  if(GetReader()->IsPileUpFromEMCal())
361  {
362  fhDeltaEtaChargedPtA3GeVPileUp[1]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
363  fhDeltaPhiChargedPtA3GeVPileUp[1]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
364  }
365  if(GetReader()->IsPileUpFromSPDOrEMCal())
366  {
367  fhDeltaEtaChargedPtA3GeVPileUp[2]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
368  fhDeltaPhiChargedPtA3GeVPileUp[2]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
369  }
370  if(GetReader()->IsPileUpFromSPDAndEMCal())
371  {
372  fhDeltaEtaChargedPtA3GeVPileUp[3]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
373  fhDeltaPhiChargedPtA3GeVPileUp[3]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
374  }
375  if(GetReader()->IsPileUpFromSPDAndNotEMCal())
376  {
377  fhDeltaEtaChargedPtA3GeVPileUp[4]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
378  fhDeltaPhiChargedPtA3GeVPileUp[4]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
379  }
380  if(GetReader()->IsPileUpFromEMCalAndNotSPD())
381  {
382  fhDeltaEtaChargedPtA3GeVPileUp[5]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
383  fhDeltaPhiChargedPtA3GeVPileUp[5]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
384  }
385  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
386  {
387  fhDeltaEtaChargedPtA3GeVPileUp[6]->Fill(ptTrig, deltaEta, GetEventWeight()) ;
388  fhDeltaPhiChargedPtA3GeVPileUp[6]->Fill(ptTrig, deltaPhi, GetEventWeight()) ;
389  }
390  }
391  }
392 
393  if(IsDataMC())
394  {
395  Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
396  fhDeltaPhiChargedMC[mcIndex]->Fill(ptTrig , deltaPhi, GetEventWeight());
397 
398  if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
399  {
400  // check index in GetMCTagIndexHistogram
401  if ( mcIndex == 2 ) fhDeltaPhiChargedMC[8]->Fill(ptTrig, deltaPhi, GetEventWeight()); // pi0 decay
402  else if( mcIndex == 4 ) fhDeltaPhiChargedMC[9]->Fill(ptTrig, deltaPhi, GetEventWeight()); // eta decay
403  }
404  }
405 
406  if(fDecayTrigger && decayTag > 0)
407  {
408  for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
409  {
410  if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
411  fhDeltaPhiDecayCharged[ibit]->Fill(ptTrig, deltaPhi, GetEventWeight());
412  }
413  }
414 
415  Double_t dphiBrad = -100;
416  if(fFillBradHisto)
417  {
418  dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
419  if( TMath::Abs(dphiBrad) > 0.325 && TMath::Abs(dphiBrad) < 0.475 ) // Hardcoded values, BAD, FIXME
420  {
421  fhAssocPtBkg->Fill(ptTrig, ptAssoc, GetEventWeight());
422  }
423 
424  if( dphiBrad < -1./3 ) dphiBrad += 2;
425  fhDeltaPhiBrad->Fill(ptTrig, dphiBrad, GetEventWeight());
426  }
427 
428  // Fill histograms in bins of associated particle pT
429  if ( bin >= 0 )
430  {
431  fhDeltaPhiDeltaEtaAssocPtBin[bin]->Fill(deltaPhi,deltaEta, GetEventWeight());
432 
433  fhDeltaPhiAssocPtBin [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
434 
435  if(fFillEtaGapsHisto)
436  {
437  if(TMath::Abs(deltaEta)> 0.8)
438  fhDeltaPhiAssocPtBinDEta08[bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
439 
440  if(TMath::Abs(deltaEta)< 0.01)
441  fhDeltaPhiAssocPtBinDEta0 [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
442  }
443 
444  if (fFillBradHisto)
445  fhDeltaPhiBradAssocPtBin [bin]->Fill(ptTrig, dphiBrad, GetEventWeight());
446 
447  if(fDecayTrigger && decayTag > 0 && fNDecayBits > 0 &&
448  GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
449  {
450  fhDeltaPhiDecayChargedAssocPtBin[bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
451  }
452 
453  if(fHMPIDCorrelation)
454  {
455  if( hmpidSignal > 0 )
456  {
457  //printf("Track pt %f with HMPID signal %f \n",pt,hmpidSignal);
458  fhDeltaPhiAssocPtBinHMPID[bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
459  }
460 
461  if(phiAssoc > 5*TMath::DegToRad() && phiAssoc < 20*TMath::DegToRad())
462  {
463  //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
464  fhDeltaPhiAssocPtBinHMPIDAcc[bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
465  }
466  }
467  }
468 
469  // Fill different multiplicity/centrality histogram
470  if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
471  {
472  fhDeltaPhiChargedMult[cen]->Fill(ptTrig, deltaPhi, GetEventWeight());
473  fhDeltaEtaChargedMult[cen]->Fill(ptTrig, deltaEta, GetEventWeight());
474  }
475 }
476 
477 //___________________________________________________________________________________________________________________________________
479 //___________________________________________________________________________________________________________________________________
480 Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float_t mcAssocPt, Float_t mcAssocPhi, Float_t mcAssocEta,
481  Float_t mcTrigPt, Float_t mcTrigPhi, Float_t mcTrigEta,
482  Int_t histoIndex, Bool_t lostDecayPair)
483 {
484  Bool_t lead = kTRUE;
485 
486  // In case we requested the trigger to be a leading particle,
487  // check if this is true at the MC level.
488  // Not sure if it is correct to skip or not skip this.
489  // Absolute leading?
490  if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt ) lead = kFALSE; // skip event
491 
492  // Skip this event if near side associated particle pt larger than trigger
493  if( mcAssocPhi < 0 ) mcAssocPhi+=TMath::TwoPi();
494 
495  Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi;
496  if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
497  if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
498 
499  if( fMakeNearSideLeading)
500  {
501  if( mcAssocPt > mcTrigPt && mcdeltaPhi < TMath::PiOver2() ) lead = kFALSE; // skip event
502  }
503 
504  //
505  // Select only hadrons in pt range
506  if ( mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt ) return lead ; // exclude but continue
507  if ( mcAssocPt < GetReader()->GetCTSPtMin()) return lead ;
508 
509 
510  //
511  // Remove trigger itself for correlation when use charged triggers
512  if(TMath::Abs(mcAssocPt -mcTrigPt ) < 1e-6 &&
513  mcdeltaPhi < 1e-6 &&
514  TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6) return lead ; // exclude but continue
515 
516  Float_t mcxE =-mcAssocPt/mcTrigPt*TMath::Cos(mcdeltaPhi);// -(mcAssocPx*pxprim+mcAssocPy*pyprim)/(mcTrigPt*mcTrigPt);
517  Float_t mchbpXE =-100 ;
518  if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
519 
520  Float_t mczT = mcAssocPt/mcTrigPt ;
521  Float_t mchbpZT =-100 ;
522  if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
523 
524  Double_t mcpout = mcAssocPt*TMath::Sin(mcdeltaPhi) ;
525 
526  AliDebug(1,Form("Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f",
527  mcAssocPt,mcAssocPhi, mcTrigPhi,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut));
528 
529  // Fill Histograms
530  fhMCEtaCharged [histoIndex]->Fill(mcAssocPt, mcAssocEta, GetEventWeight());
531  fhMCPhiCharged [histoIndex]->Fill(mcAssocPt, mcAssocPhi, GetEventWeight());
532  fhMCDeltaPhiCharged[histoIndex]->Fill(mcTrigPt , mcdeltaPhi, GetEventWeight());
533  fhMCPtAssocDeltaPhi[histoIndex]->Fill(mcAssocPt, mcdeltaPhi, GetEventWeight());
534 
535  fhMCDeltaEtaCharged [histoIndex]->Fill(mcTrigPt , mcTrigEta-mcAssocEta, GetEventWeight());
536  fhMCDeltaPhiDeltaEtaCharged[histoIndex]->Fill(mcdeltaPhi, mcTrigEta-mcAssocEta, GetEventWeight());
537 
538  // Delta phi cut for correlation
539  if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
540  {
541  fhMCDeltaPhiChargedPt[histoIndex]->Fill(mcAssocPt, mcdeltaPhi, GetEventWeight());
542  fhMCPtXECharged [histoIndex]->Fill(mcTrigPt , mcxE , GetEventWeight());
543  fhMCPtHbpXECharged [histoIndex]->Fill(mcTrigPt , mchbpXE , GetEventWeight());
544  fhMCPtZTCharged [histoIndex]->Fill(mcTrigPt , mczT , GetEventWeight());
545  fhMCPtHbpZTCharged [histoIndex]->Fill(mcTrigPt , mchbpZT , GetEventWeight());
546  fhMCPtTrigPout [histoIndex]->Fill(mcTrigPt , mcpout , GetEventWeight());
547  }
548 
549  if(lostDecayPair)
550  {
551  // check index in GetMCTagIndexHistogram
552  if(histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
553  {
554  // pi0 decay
555  fhMCEtaCharged [8]->Fill(mcAssocPt, mcAssocEta, GetEventWeight());
556  fhMCPhiCharged [8]->Fill(mcAssocPt, mcAssocPhi, GetEventWeight());
557  fhMCDeltaPhiCharged[8]->Fill(mcTrigPt , mcdeltaPhi, GetEventWeight());
558  fhMCPtAssocDeltaPhi[8]->Fill(mcAssocPt, mcdeltaPhi, GetEventWeight());
559 
560  fhMCDeltaEtaCharged [8]->Fill(mcTrigPt , mcTrigEta-mcAssocEta, GetEventWeight());
561  fhMCDeltaPhiDeltaEtaCharged[8]->Fill(mcdeltaPhi, mcTrigEta-mcAssocEta, GetEventWeight());
562 
563  //delta phi cut for correlation
564  if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
565  {
566  fhMCDeltaPhiChargedPt[8]->Fill(mcAssocPt, mcdeltaPhi, GetEventWeight());
567  fhMCPtXECharged [8]->Fill(mcTrigPt , mcxE , GetEventWeight());
568  fhMCPtHbpXECharged [8]->Fill(mcTrigPt , mchbpXE , GetEventWeight());
569  fhMCPtZTCharged [8]->Fill(mcTrigPt , mczT , GetEventWeight());
570  fhMCPtHbpZTCharged [8]->Fill(mcTrigPt , mchbpZT , GetEventWeight());
571  fhMCPtTrigPout [8]->Fill(mcTrigPt , mcpout , GetEventWeight());
572  }
573  } // pi0 decay lost pair
574 
575  if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
576  {
577  // eta decay
578  fhMCEtaCharged [9]->Fill(mcAssocPt, mcAssocEta, GetEventWeight());
579  fhMCPhiCharged [9]->Fill(mcAssocPt, mcAssocPhi, GetEventWeight());
580  fhMCDeltaPhiCharged[9]->Fill(mcTrigPt , mcdeltaPhi, GetEventWeight());
581  fhMCPtAssocDeltaPhi[9]->Fill(mcAssocPt, mcdeltaPhi, GetEventWeight());
582 
583  fhMCDeltaEtaCharged [9]->Fill(mcTrigPt , mcTrigEta-mcAssocEta, GetEventWeight());
584  fhMCDeltaPhiDeltaEtaCharged[9]->Fill(mcdeltaPhi, mcTrigEta-mcAssocEta, GetEventWeight());
585 
586  // Delta phi cut for correlation
587  if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
588  {
589  fhMCDeltaPhiChargedPt[9]->Fill(mcAssocPt, mcdeltaPhi, GetEventWeight());
590  fhMCPtXECharged [9]->Fill(mcTrigPt , mcxE , GetEventWeight());
591  fhMCPtHbpXECharged [9]->Fill(mcTrigPt , mchbpXE , GetEventWeight());
592  fhMCPtZTCharged [9]->Fill(mcTrigPt , mczT , GetEventWeight());
593  fhMCPtHbpZTCharged [9]->Fill(mcTrigPt , mchbpZT , GetEventWeight());
594  fhMCPtTrigPout [9]->Fill(mcTrigPt , mcpout , GetEventWeight()) ;
595  }
596  } // eta decay lost pair
597  }
598  // Underlying event
599 
600  // Right
601  if ( (mcdeltaPhi > fUeDeltaPhiMinCut) && (mcdeltaPhi < fUeDeltaPhiMaxCut) )
602  {
603  //Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
604  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
605  Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
606  Double_t mcUezT = mcAssocPt/mcTrigPt;
607 
608  if(mcUexE < 0.)
609  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
610  mcUexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
611 
612  fhMCPtXEUeCharged[histoIndex]->Fill(mcTrigPt, mcUexE, GetEventWeight());
613  if(mcUexE > 0) fhMCPtHbpXEUeCharged[histoIndex]->Fill(mcTrigPt, TMath::Log(1/mcUexE), GetEventWeight());
614 
615  fhMCPtZTUeCharged[histoIndex]->Fill(mcTrigPt, mcUezT, GetEventWeight());
616  if(mcUezT > 0) fhMCPtHbpZTUeCharged[histoIndex]->Fill(mcTrigPt, TMath::Log(1/mcUezT), GetEventWeight());
617 
618  fhMCUePart[histoIndex]->Fill(mcTrigPt, GetEventWeight());
619 
620  if(lostDecayPair)
621  {
622  // check index in GetMCTagIndexHistogram
623  if(histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
624  {
625  // pi0 decay
626  fhMCPtXEUeCharged[8]->Fill(mcTrigPt, mcUexE, GetEventWeight());
627  if(mcUexE > 0) fhMCPtHbpXEUeCharged[8]->Fill(mcTrigPt, TMath::Log(1/mcUexE), GetEventWeight());
628 
629  fhMCPtZTUeCharged[8]->Fill(mcTrigPt, mcUezT, GetEventWeight());
630  if(mcUezT > 0) fhMCPtHbpZTUeCharged[8]->Fill(mcTrigPt, TMath::Log(1/mcUezT), GetEventWeight());
631 
632  fhMCUePart[8]->Fill(mcTrigPt, GetEventWeight());
633  }
634  else if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
635  {
636  // eta decay
637  fhMCPtXEUeCharged[9]->Fill(mcTrigPt, mcUexE, GetEventWeight());
638  if(mcUexE > 0) fhMCPtHbpXEUeCharged[9]->Fill(mcTrigPt, TMath::Log(1/mcUexE), GetEventWeight());
639 
640  fhMCPtZTUeCharged[9]->Fill(mcTrigPt, mcUezT, GetEventWeight());
641  if(mcUezT > 0) fhMCPtHbpZTUeCharged[9]->Fill(mcTrigPt, TMath::Log(1/mcUezT), GetEventWeight());
642 
643  fhMCUePart[9]->Fill(mcTrigPt, GetEventWeight());
644  }
645  }
646  }
647 
648  if(fMakeSeveralUE)
649  {
650  // Left
651  if((mcdeltaPhi<-fUeDeltaPhiMinCut) || (mcdeltaPhi >2*fUeDeltaPhiMaxCut))
652  {
653  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
654  Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
655  Double_t mcUezT = mcAssocPt/mcTrigPt;
656 
657  if(mcUexE < 0.)
658  AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
659  mcUexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
660 
661  fhMCPtXEUeLeftCharged[histoIndex]->Fill(mcTrigPt, mcUexE, GetEventWeight());
662  if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[histoIndex]->Fill(mcTrigPt, TMath::Log(1/mcUexE), GetEventWeight());
663 
664  fhMCPtZTUeLeftCharged[histoIndex]->Fill(mcTrigPt, mcUezT, GetEventWeight());
665  if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[histoIndex]->Fill(mcTrigPt, TMath::Log(1/mcUezT), GetEventWeight());
666 
667  if(lostDecayPair)
668  {
669  // Check index in GetMCTagIndexHistogram
670  if(histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
671  {
672  // pi0 decay
673  fhMCPtXEUeLeftCharged[8]->Fill(mcTrigPt,mcUexE, GetEventWeight());
674  if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[8]->Fill(mcTrigPt, TMath::Log(1/mcUexE), GetEventWeight());
675 
676  fhMCPtZTUeLeftCharged[8]->Fill(mcTrigPt,mcUezT, GetEventWeight());
677  if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[8]->Fill(mcTrigPt, TMath::Log(1/mcUezT), GetEventWeight());
678  }
679  else if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
680  {
681  // eta decay
682  fhMCPtXEUeLeftCharged[9]->Fill(mcTrigPt, mcUexE, GetEventWeight());
683  if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[9]->Fill(mcTrigPt, TMath::Log(1/mcUexE), GetEventWeight());
684 
685  fhMCPtZTUeLeftCharged[9]->Fill(mcTrigPt, mcUezT, GetEventWeight());
686  if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[9]->Fill(mcTrigPt, TMath::Log(1/mcUezT), GetEventWeight());
687  }
688  }
689  }
690  }
691 
692  return lead;
693 }
694 
695 //______________________________________________________________________________________________________________
697 //______________________________________________________________________________________________________________
698 void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Float_t ptTrig, Float_t ptAssoc,
699  Float_t deltaPhi,
700  Int_t cen, Int_t charge,
701  Int_t bin, Int_t decayTag,
702  Int_t outTOF, Int_t mcTag)
703 
704 {
705  Float_t zT = ptAssoc/ptTrig ;
706  Float_t xE =-ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
707  Float_t pout = ptAssoc*TMath::Sin(deltaPhi) ;
708 
709  if(xE < 0.)
710  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
711  xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
712 
713  Float_t hbpXE = -100;
714  Float_t hbpZT = -100;
715 
716  if(xE > 0 ) hbpXE = TMath::Log(1./xE);
717  if(zT > 0 ) hbpZT = TMath::Log(1./zT);
718 
719  fhXECharged ->Fill(ptTrig, xE , GetEventWeight());
720  fhPtHbpXECharged ->Fill(ptTrig, hbpXE , GetEventWeight());
721  fhZTCharged ->Fill(ptTrig, zT , GetEventWeight());
722  fhPtHbpZTCharged ->Fill(ptTrig, hbpZT , GetEventWeight());
723  fhPtTrigPout ->Fill(ptTrig, pout , GetEventWeight());
724  fhPtTrigCharged ->Fill(ptTrig, ptAssoc, GetEventWeight());
725 
726  if((deltaPhi > 5*TMath::Pi()/6.) && (deltaPhi < 7*TMath::Pi()/6.))
727  {
728  fhXECharged_Cone2 ->Fill(ptTrig, xE , GetEventWeight());
729  fhPtHbpXECharged_Cone2 ->Fill(ptTrig, hbpXE, GetEventWeight());
730  }
731 
732  // MC
733  if(IsDataMC())
734  {
735  Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
736  fhXEChargedMC[mcIndex]->Fill(ptTrig, xE, GetEventWeight());
737 
738  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
739  {
740  // check index in GetMCTagIndexHistogram
741  if ( mcIndex == 2 ) fhXEChargedMC[8]->Fill(ptTrig, xE, GetEventWeight()); // pi0 decay
742  else if ( mcIndex == 4 ) fhXEChargedMC[9]->Fill(ptTrig, xE, GetEventWeight()); // eta decay
743  }
744  }
745 
746  // Pile-up studies
747  if(IsPileUpAnalysisOn())
748  {
749  if (outTOF==1)
750  {
751  fhXEChargedOtherBC ->Fill(ptTrig, xE , GetEventWeight());
752  fhZTChargedOtherBC ->Fill(ptTrig, zT , GetEventWeight());
753  fhPtTrigChargedOtherBC->Fill(ptTrig, ptAssoc, GetEventWeight());
754  }
755  else if(outTOF==0)
756  {
757  fhXEChargedBC0 ->Fill(ptTrig, xE , GetEventWeight());
758  fhZTChargedBC0 ->Fill(ptTrig, zT , GetEventWeight());
759  fhPtTrigChargedBC0->Fill(ptTrig, ptAssoc, GetEventWeight());
760  }
761 
762  Int_t vtxBC = GetReader()->GetVertexBC();
763  if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
764  {
765  fhXEChargedVtxBC0 ->Fill(ptTrig, xE , GetEventWeight());
766  fhZTChargedVtxBC0 ->Fill(ptTrig, zT , GetEventWeight());
767  fhPtTrigChargedVtxBC0->Fill(ptTrig, ptAssoc, GetEventWeight());
768  }
769 
770  if(GetReader()->IsPileUpFromSPD())
771  {
772  fhXEChargedPileUp [0]->Fill(ptTrig, xE , GetEventWeight());
773  fhZTChargedPileUp [0]->Fill(ptTrig, zT , GetEventWeight());
774  fhPtTrigChargedPileUp[0]->Fill(ptTrig, ptAssoc, GetEventWeight());
775  }
776  if(GetReader()->IsPileUpFromEMCal())
777  {
778  fhXEChargedPileUp [1]->Fill(ptTrig, xE , GetEventWeight());
779  fhZTChargedPileUp [1]->Fill(ptTrig, zT , GetEventWeight());
780  fhPtTrigChargedPileUp[1]->Fill(ptTrig, ptAssoc, GetEventWeight());
781  }
782  if(GetReader()->IsPileUpFromSPDOrEMCal())
783  {
784  fhXEChargedPileUp [2]->Fill(ptTrig, xE , GetEventWeight());
785  fhZTChargedPileUp [2]->Fill(ptTrig, zT , GetEventWeight());
786  fhPtTrigChargedPileUp[2]->Fill(ptTrig, ptAssoc, GetEventWeight());
787  }
788  if(GetReader()->IsPileUpFromSPDAndEMCal())
789  {
790  fhXEChargedPileUp [3]->Fill(ptTrig, xE , GetEventWeight());
791  fhZTChargedPileUp [3]->Fill(ptTrig, zT , GetEventWeight());
792  fhPtTrigChargedPileUp[3]->Fill(ptTrig, ptAssoc, GetEventWeight());
793  }
794  if(GetReader()->IsPileUpFromSPDAndNotEMCal())
795  {
796  fhXEChargedPileUp [4]->Fill(ptTrig, xE , GetEventWeight());
797  fhZTChargedPileUp [4]->Fill(ptTrig, zT , GetEventWeight());
798  fhPtTrigChargedPileUp[4]->Fill(ptTrig, ptAssoc, GetEventWeight());
799  }
800  if(GetReader()->IsPileUpFromEMCalAndNotSPD())
801  {
802  fhXEChargedPileUp [5]->Fill(ptTrig, xE , GetEventWeight());
803  fhZTChargedPileUp [5]->Fill(ptTrig, zT , GetEventWeight());
804  fhPtTrigChargedPileUp[5]->Fill(ptTrig, ptAssoc, GetEventWeight());
805  }
806  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
807  {
808  fhXEChargedPileUp [6]->Fill(ptTrig, xE , GetEventWeight());
809  fhZTChargedPileUp [6]->Fill(ptTrig, zT , GetEventWeight());
810  fhPtTrigChargedPileUp[6]->Fill(ptTrig, ptAssoc, GetEventWeight());
811  }
812  }
813 
814  if(fDecayTrigger && decayTag > 0)
815  {
816  for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
817  {
818  if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
819  {
820  fhXEDecayCharged[ibit]->Fill(ptTrig, xE, GetEventWeight());
821  fhZTDecayCharged[ibit]->Fill(ptTrig, zT, GetEventWeight());
822  }
823  }
824  } // photon decay pi0/eta trigger
825 
826  if(bin >= 0 && fFillMomImbalancePtAssocBinsHisto)//away side
827  {
828  fhXEAssocPtBin[bin]->Fill(ptTrig, xE, GetEventWeight()) ;
829  fhZTAssocPtBin[bin]->Fill(ptTrig, zT, GetEventWeight()) ;
830  }
831 
832  if(fCorrelVzBin)
833  {
834  Int_t vz = GetEventVzBin();
835  fhXEVZ[vz]->Fill(ptTrig, xE, GetEventWeight()) ;
836  fhZTVZ[vz]->Fill(ptTrig, zT, GetEventWeight()) ;
837  }
838 
839  if(charge > 0)
840  {
841  fhXEPosCharged->Fill(ptTrig, xE, GetEventWeight()) ;
842  fhZTPosCharged->Fill(ptTrig, zT, GetEventWeight()) ;
843  }
844  else
845  {
846  fhXENegCharged->Fill(ptTrig, xE, GetEventWeight()) ;
847  fhZTNegCharged->Fill(ptTrig, zT, GetEventWeight()) ;
848  }
849 
850  // Fill different multiplicity/centrality histogram
851  if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
852  {
853  fhXEMult[cen]->Fill(ptTrig, xE, GetEventWeight());
854  fhZTMult[cen]->Fill(ptTrig, zT, GetEventWeight());
855  } // multiplicity/centrality events selection
856 }
857 
858 //___________________________________________________________________________________________________________________
860 //___________________________________________________________________________________________________________________
861 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(Float_t ptTrig, Float_t ptAssoc,
862  Float_t deltaPhi, Int_t cen, Int_t outTOF, Int_t mcTag)
863 {
864  fhUePart->Fill(ptTrig, GetEventWeight());
865 
866  fhDeltaPhiUeChargedPt->Fill(ptAssoc, deltaPhi, GetEventWeight());
867 
868  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
869  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
870  Double_t uezT = ptAssoc/ptTrig;
871 
872  if(uexE < 0.)
873  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
874  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
875 
876  fhXEUeCharged->Fill(ptTrig, uexE, GetEventWeight());
877  if(deltaPhi > TMath::DegToRad()*80 && deltaPhi < TMath::DegToRad()*100)fhXEUeChargedSmallCone->Fill(ptTrig, uexE, GetEventWeight());
878  if(deltaPhi > TMath::DegToRad()*70 && deltaPhi < TMath::DegToRad()*110)fhXEUeChargedMediumCone->Fill(ptTrig, uexE, GetEventWeight());
879  if(deltaPhi > TMath::DegToRad()*60 && deltaPhi < TMath::DegToRad()*120)fhXEUeChargedLargeCone->Fill(ptTrig, uexE, GetEventWeight());
880 
881  if(uexE > 0) fhPtHbpXEUeCharged->Fill(ptTrig, TMath::Log(1/uexE), GetEventWeight());
882 
883  fhZTUeCharged->Fill(ptTrig, uezT, GetEventWeight());
884  if(uezT > 0) fhPtHbpZTUeCharged->Fill(ptTrig, TMath::Log(1/uezT), GetEventWeight());
885 
886  if(IsDataMC())
887  {
888  Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
889  fhXEUeChargedRightMC[mcIndex]->Fill(ptTrig, uexE, GetEventWeight());
890 
891  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
892  {
893  // check index in GetMCTagIndexHistogram
894  if ( mcIndex == 2 ) fhXEUeChargedRightMC[8]->Fill(ptTrig, uexE, GetEventWeight()); // pi0 decay
895  else if ( mcIndex == 4 ) fhXEUeChargedRightMC[9]->Fill(ptTrig, uexE, GetEventWeight()); // eta decay
896  }
897  }
898 
899  // Pile-up studies
900 
901  if(IsPileUpAnalysisOn())
902  {
903  if (outTOF==1)
904  {
905  fhXEUeChargedOtherBC->Fill(ptTrig, uexE, GetEventWeight());
906  fhZTUeChargedOtherBC->Fill(ptTrig, uezT, GetEventWeight());
907  }
908  else if(outTOF==0)
909  {
910  fhXEUeChargedBC0->Fill(ptTrig, uexE, GetEventWeight());
911  fhZTUeChargedBC0->Fill(ptTrig, uezT, GetEventWeight());
912  }
913 
914  Int_t vtxBC = GetReader()->GetVertexBC();
915  if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
916  {
917  fhXEUeChargedVtxBC0->Fill(ptTrig, uexE, GetEventWeight());
918  fhZTUeChargedVtxBC0->Fill(ptTrig, uezT, GetEventWeight());
919  }
920 
921  if(GetReader()->IsPileUpFromSPD())
922  {
923  fhXEUeChargedPileUp[0]->Fill(ptTrig, uexE, GetEventWeight());
924  fhZTUeChargedPileUp[0]->Fill(ptTrig, uezT, GetEventWeight());
925  }
926  if(GetReader()->IsPileUpFromEMCal())
927  {
928  fhXEUeChargedPileUp[1]->Fill(ptTrig, uexE, GetEventWeight());
929  fhZTUeChargedPileUp[1]->Fill(ptTrig, uezT, GetEventWeight());
930  }
931  if(GetReader()->IsPileUpFromSPDOrEMCal())
932  {
933  fhXEUeChargedPileUp[2]->Fill(ptTrig, uexE, GetEventWeight());
934  fhZTUeChargedPileUp[2]->Fill(ptTrig, uezT, GetEventWeight());
935  }
936  if(GetReader()->IsPileUpFromSPDAndEMCal())
937  {
938  fhXEUeChargedPileUp[3]->Fill(ptTrig, uexE, GetEventWeight());
939  fhZTUeChargedPileUp[3]->Fill(ptTrig, uezT, GetEventWeight());
940  }
941  if(GetReader()->IsPileUpFromSPDAndNotEMCal())
942  {
943  fhXEUeChargedPileUp[4]->Fill(ptTrig, uexE, GetEventWeight());
944  fhZTUeChargedPileUp[4]->Fill(ptTrig, uezT, GetEventWeight());
945  }
946  if(GetReader()->IsPileUpFromEMCalAndNotSPD())
947  {
948  fhXEUeChargedPileUp[5]->Fill(ptTrig, uexE, GetEventWeight());
949  fhZTUeChargedPileUp[5]->Fill(ptTrig, uezT, GetEventWeight());
950  }
951  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
952  {
953  fhXEUeChargedPileUp[6]->Fill(ptTrig, uexE, GetEventWeight());
954  fhZTUeChargedPileUp[6]->Fill(ptTrig, uezT, GetEventWeight());
955  }
956  }
957 
958  // Fill different multiplicity/centrality histogram
959  if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
960  {
961  fhXEUeMult[cen]->Fill(ptTrig, uexE, GetEventWeight());
962  fhZTUeMult[cen]->Fill(ptTrig, uezT, GetEventWeight());
963  }
964 }
965 
966 //_______________________________________________________________________________________________
969 //_______________________________________________________________________________________________
970 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(Float_t ptTrig,
971  Float_t ptAssoc,
972  Float_t deltaPhi,
973  Int_t mcTag)
974 {
975  if((deltaPhi<-fUeDeltaPhiMinCut) || (deltaPhi >2*fUeDeltaPhiMaxCut))
976  {
977  fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi, GetEventWeight());
978 
979  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
980  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
981  Double_t uezT = ptAssoc/ptTrig;
982 
983  if(uexE < 0.)
984  AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
985  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
986 
987  fhXEUeLeftCharged->Fill(ptTrig, uexE, GetEventWeight());
988  if(uexE > 0) fhPtHbpXEUeLeftCharged->Fill(ptTrig, TMath::Log(1/uexE), GetEventWeight());
989 
990  fhZTUeLeftCharged->Fill(ptTrig, uezT, GetEventWeight());
991  if(uezT > 0) fhPtHbpZTUeLeftCharged->Fill(ptTrig, TMath::Log(1/uezT), GetEventWeight());
992 
993  fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi, GetEventWeight());
994 
995  if(IsDataMC())
996  {
997  Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
998  fhXEUeChargedLeftMC[mcIndex]->Fill(ptTrig, uexE, GetEventWeight());
999 
1000  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
1001  {
1002  // check index in GetMCTagIndexHistogram
1003  if ( mcIndex == 2 ) fhXEUeChargedLeftMC[8]->Fill(ptTrig, uexE, GetEventWeight()); // pi0 decay
1004  else if ( mcIndex == 4 ) fhXEUeChargedLeftMC[9]->Fill(ptTrig, uexE, GetEventWeight()); // eta decay
1005  }
1006  }
1007  }
1008 
1009  if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-TMath::Pi()/2))
1010  {
1011  fhDeltaPhiUeLeftDownCharged->Fill(ptAssoc, deltaPhi, GetEventWeight());
1012 
1013  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
1014  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
1015 
1016  if(uexE < 0.)
1017  AliWarning(Form("Careful!!, negative xE %2.2f for left-down UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
1018  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
1019 
1020  fhXEUeLeftDownCharged->Fill(ptTrig, uexE, GetEventWeight());
1021  }
1022 
1023  if((deltaPhi>2*fUeDeltaPhiMaxCut) && (deltaPhi <3*TMath::Pi()/2))
1024  {
1025  fhDeltaPhiUeLeftUpCharged->Fill(ptAssoc, deltaPhi, GetEventWeight());
1026 
1027  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
1028  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
1029 
1030  if(uexE < 0.)
1031  AliWarning(Form("Careful!!, negative xE %2.2f for left-up UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
1032  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
1033 
1034  fhXEUeLeftUpCharged->Fill(ptTrig, uexE, GetEventWeight());
1035  }
1036 
1037  if((deltaPhi > TMath::Pi()/2) && (deltaPhi < fUeDeltaPhiMaxCut))
1038  {
1039  fhDeltaPhiUeRightUpCharged->Fill(ptAssoc, deltaPhi, GetEventWeight());
1040 
1041  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
1042  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
1043 
1044  if(uexE < 0.)
1045  AliWarning(Form("Careful!!, negative xE %2.2f for right-up UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
1046  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
1047 
1048  fhXEUeRightUpCharged->Fill(ptTrig, uexE, GetEventWeight());
1049  }
1050 
1051  if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < TMath::Pi()/2))
1052  {
1053  fhDeltaPhiUeRightDownCharged->Fill(ptAssoc, deltaPhi, GetEventWeight());
1054 
1055  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
1056  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
1057 
1058  if(uexE < 0.)
1059  AliWarning(Form("Careful!!, negative xE %2.2f for right-down UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
1060  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
1061 
1062  fhXEUeRightDownCharged->Fill(ptTrig, uexE, GetEventWeight());
1063  }
1064 }
1065 
1066 //_____________________________________________________________________________________________________________________________________
1068 //_____________________________________________________________________________________________________________________________________
1069 void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float_t ptAssoc, Float_t phiAssoc, Bool_t bChargedOrNeutral)
1070 {
1071  // Calculate the correlation parameters
1072  Float_t ptDecay1 = fDecayMom1.Pt();
1073  Float_t ptDecay2 = fDecayMom2.Pt();
1074 
1075  Float_t zTDecay1 = -100, zTDecay2 = -100;
1076  if(ptDecay1 > 0) zTDecay1 = ptAssoc/ptDecay1 ;
1077  if(ptDecay2 > 0) zTDecay2 = ptAssoc/ptDecay2 ;
1078 
1079  Float_t deltaPhiDecay1 = fDecayMom1.Phi()-phiAssoc;
1080  if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
1081  if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
1082 
1083  Float_t deltaPhiDecay2 = fDecayMom2.Phi()-phiAssoc;
1084  if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
1085  if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
1086 
1087  Float_t xEDecay1 =-zTDecay1*TMath::Cos(deltaPhiDecay1); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1088  Float_t xEDecay2 =-zTDecay2*TMath::Cos(deltaPhiDecay2); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1089 
1090  if(bChargedOrNeutral) // correlate with charges
1091  {
1092  fhDeltaPhiPi0DecayCharged->Fill(ptDecay1, deltaPhiDecay1, GetEventWeight());
1093  fhDeltaPhiPi0DecayCharged->Fill(ptDecay2, deltaPhiDecay2, GetEventWeight());
1094 
1095  AliDebug(2,Form("deltaPhoton1 = %f, deltaPhoton2 = %f", deltaPhiDecay1, deltaPhiDecay2));
1096 
1097  if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
1098  {
1099  fhZTPi0DecayCharged->Fill(ptDecay1, zTDecay1, GetEventWeight());
1100  fhXEPi0DecayCharged->Fill(ptDecay1, xEDecay1, GetEventWeight());
1101  }
1102  if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
1103  {
1104  fhZTPi0DecayCharged->Fill(ptDecay2, zTDecay2, GetEventWeight());
1105  fhXEPi0DecayCharged->Fill(ptDecay2, xEDecay2, GetEventWeight());
1106  }
1107  }
1108  else // correlate with neutrals
1109  {
1110  fhDeltaPhiPi0DecayNeutral->Fill(ptDecay1, deltaPhiDecay1, GetEventWeight());
1111  fhDeltaPhiPi0DecayNeutral->Fill(ptDecay2, deltaPhiDecay2, GetEventWeight());
1112 
1113  AliDebug(2,Form("deltaPhoton1 = %f, deltaPhoton2 = %f", deltaPhiDecay1, deltaPhiDecay2));
1114 
1115  if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
1116  {
1117  fhZTPi0DecayNeutral->Fill(ptDecay1, zTDecay1, GetEventWeight());
1118  fhXEPi0DecayNeutral->Fill(ptDecay1, xEDecay1, GetEventWeight());
1119  }
1120  if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
1121  {
1122  fhZTPi0DecayNeutral->Fill(ptDecay2, zTDecay2, GetEventWeight());
1123  fhXEPi0DecayNeutral->Fill(ptDecay2, xEDecay2, GetEventWeight());
1124  }
1125  }
1126 }
1127 
1128 //_____________________________________________________________________________________________________________________________
1131 //_____________________________________________________________________________________________________________________________
1132 void AliAnaParticleHadronCorrelation::FillNeutralUnderlyingEventSidesHistograms(Float_t ptTrig, Float_t ptAssoc,
1133  Float_t zT, Float_t hbpZT,
1134  Float_t deltaPhi)
1135 {
1136  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
1137 
1138  Float_t xE =-ptAssoc/ptTrig*TMath::Cos(randomphi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1139  Float_t hbpXE = -100;
1140  if(xE > 0 ) hbpXE = TMath::Log(1./xE);
1141 
1142  if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
1143  {
1144  fhDeltaPhiUeLeftNeutral->Fill(ptAssoc, deltaPhi, GetEventWeight());
1145  fhXEUeLeftNeutral ->Fill(ptTrig , xE , GetEventWeight());
1146  fhPtHbpXEUeLeftNeutral ->Fill(ptTrig , hbpXE , GetEventWeight());
1147  fhZTUeLeftNeutral ->Fill(ptTrig , zT , GetEventWeight());
1148  fhPtHbpZTUeLeftNeutral ->Fill(ptTrig , hbpZT , GetEventWeight());
1149  }
1150 }
1151 
1152 //______________________________________________________
1154 //______________________________________________________
1155 void AliAnaParticleHadronCorrelation::FillEventMixPool()
1156 {
1157  if ( !DoOwnMix() ) return;
1158 
1159  FillChargedEventMixPool();
1160 
1161  // Do the cluster pool filling only if requested
1162  // or in case of isolation cut using clusters in the cone.
1163  Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
1164 
1165  if( !fFillNeutralEventMixPool && !isoCase) return;
1166 
1167  FillNeutralEventMixPool();
1168 }
1169 
1170 //_____________________________________________________________
1172 //_____________________________________________________________
1173 void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
1174 {
1175  if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
1176  {
1177  //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
1178  return ; // pool filled previously for another trigger
1179  }
1180 
1181  AliAnalysisManager * manager = AliAnalysisManager::GetAnalysisManager();
1182  AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
1183 
1184  if(!inputHandler) return ;
1185 
1186  // Do mixing only with MB event (or the chosen mask), if not skip
1187  if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
1188 
1189  Int_t eventBin = GetEventMixBin();
1190 
1191  //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
1192  if(eventBin < 0) return;
1193 
1194  fhEventMBBin->Fill(eventBin, GetEventWeight());
1195 
1196  TObjArray * mixEventTracks = new TObjArray;
1197 
1198  if(fUseMixStoredInReader)
1199  {
1200  fListMixTrackEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
1201  }
1202 
1203  if(!fListMixTrackEvents[eventBin]) fListMixTrackEvents[eventBin] = new TList();
1204 
1205  //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
1206 
1207  TList * pool = fListMixTrackEvents[eventBin];
1208 
1209  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
1210  {
1211  AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1212 
1213  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1214  Float_t pt = fTrackVector.Pt();
1215 
1216  // Select only hadrons in pt range
1217  if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1218 
1219  AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),0);
1220  mixedTrack->SetDetectorTag(kCTS);
1221  mixedTrack->SetChargedBit(track->Charge()>0);
1222  mixEventTracks->Add(mixedTrack);
1223  }
1224 
1225  fhNtracksMB->Fill(mixEventTracks->GetEntriesFast(), eventBin, GetEventWeight());
1226 
1227  // Set the event number where the last event was added, to avoid double pool filling
1228  GetReader()->SetLastTracksMixedEvent(GetEventNumber());
1229 
1230  //printf("Add event to pool with %d tracks \n ",mixEventTracks->GetEntries());
1231  pool->AddFirst(mixEventTracks);
1232  mixEventTracks = 0;
1233 
1234  //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
1235 
1236  if(pool->GetSize() > GetNMaxEvMix())
1237  {
1238  // Remove last event
1239  TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
1240  pool->RemoveLast() ;
1241  delete tmp ;
1242  }
1243 }
1244 
1245 //_____________________________________________________________
1248 //_____________________________________________________________
1249 void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
1250 {
1251  //printf("FillNeutralEventMixPool for %s\n",GetInputAODName().Data());
1252 
1253  if(fUseMixStoredInReader && GetReader()->GetLastCaloMixedEvent() == GetEventNumber())
1254  {
1255  //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
1256  return ; // pool filled previously for another trigger
1257  }
1258 
1259  TObjArray * pl = GetEMCALClusters();
1260  //if (GetAODObjArrayName.Contains("PHOS") )pl = GetPHOSClusters();
1261  //else pl = GetEMCALClusters();
1262 
1263  AliAnalysisManager * manager = AliAnalysisManager::GetAnalysisManager();
1264  AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
1265 
1266  if(!inputHandler) return ;
1267 
1268  // Do mixing only with MB event (or the chosen mask), if not skip
1269  if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
1270 
1271  Int_t eventBin = GetEventMixBin();
1272 
1273  // Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
1274  if(eventBin < 0) return;
1275 
1276  TObjArray * mixEventCalo = new TObjArray;
1277 
1278  if(fUseMixStoredInReader)
1279  {
1280  fListMixCaloEvents[eventBin] = GetReader()->GetListWithMixedEventsForCalo(eventBin);
1281  }
1282 
1283  if(!fListMixCaloEvents[eventBin]) fListMixCaloEvents[eventBin] = new TList();
1284 
1285  TList * poolCalo = fListMixCaloEvents[eventBin];
1286 
1287  for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
1288  {
1289  AliVCluster * calo = (AliVCluster *) (pl->At(ipr)) ;
1290 
1291  // Remove matched clusters
1292  if( IsTrackMatched( calo, GetReader()->GetInputEvent() ) ) continue ;
1293 
1294  // Cluster momentum calculation
1295  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1296  {
1297  calo->GetMomentum(fMomentum,GetVertex(0)) ;
1298  }// Assume that come from vertex in straight line
1299  else
1300  {
1301  Double_t vertex[]={0,0,0};
1302  calo->GetMomentum(fMomentum,vertex) ;
1303  }
1304 
1305  Float_t pt = fMomentum.Pt();
1306 
1307  // Select only clusters in pt range
1308  if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1309 
1310  AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(fMomentum);
1311  mixedCalo->SetDetectorTag(kEMCAL);
1312  mixEventCalo->Add(mixedCalo);
1313  }
1314 
1315  fhNclustersMB->Fill(mixEventCalo->GetEntriesFast(), eventBin, GetEventWeight());
1316 
1317  // Set the event number where the last event was added, to avoid double pool filling
1318  GetReader()->SetLastCaloMixedEvent(GetEventNumber());
1319 
1320  //printf("Add event to pool with %d clusters \n ",mixEventCalo->GetEntries());
1321  poolCalo->AddFirst(mixEventCalo);
1322  mixEventCalo = 0;
1323 
1324  //printf("Pool size %d, max %d\n",poolCalo->GetSize(), GetNMaxEvMix());
1325 
1326  if(poolCalo->GetSize() > GetNMaxEvMix())
1327  {
1328  // Remove last event
1329  TClonesArray * tmp = static_cast<TClonesArray*>(poolCalo->Last()) ;
1330  poolCalo->RemoveLast() ;
1331  delete tmp ;
1332  }
1333 }
1334 
1335 //_________________________________________________________________________________________________________________
1338 //_________________________________________________________________________________________________________________
1339 Bool_t AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow(AliAODPWG4ParticleCorrelation * particle)
1340 {
1341  Float_t etaTrig = particle->Eta();
1342  Float_t ptTrig = particle->Pt();
1343  Float_t phiTrig = particle->Phi();
1344  if(phiTrig < 0 ) phiTrig+= TMath::TwoPi();
1345 
1346  Float_t ptLeadHad = 0 ;
1347  Float_t dphiLeadHad= -100 ;
1348  Float_t phiLeadHad = -100 ;
1349  Float_t etaLeadHad = -100 ;
1350  Int_t nTrack = 0;
1351 
1352  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
1353  {
1354  AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1355 
1356  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1357 
1358  Float_t pt = fTrackVector.Pt();
1359  Float_t phi = fTrackVector.Phi() ;
1360  if(phi < 0 ) phi+= TMath::TwoPi();
1361 
1362  Float_t deltaPhi = phiTrig-phi;
1363  //
1364  // Calculate deltaPhi shift so that for the particles on the opposite side
1365  // it is defined between 90 and 270 degrees
1366  // Shift [-360,-90] to [0, 270]
1367  // and [270,360] to [-90,0]
1368  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1369  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
1370 
1371  if(pt > ptLeadHad && deltaPhi > TMath::PiOver2()) // in opposite hemisphere
1372  {
1373  ptLeadHad = pt ;
1374  phiLeadHad = phi;
1375  dphiLeadHad= deltaPhi;
1376  etaLeadHad = fTrackVector.Eta();
1377  nTrack++;
1378  }
1379  }// track loop
1380 
1381  if(fFillLeadHadOppositeHisto)
1382  {
1383  if(nTrack == 0)
1384  {
1385  fhPtNoLeadingOppositeHadron ->Fill(ptTrig , GetEventWeight());
1386  fhEtaPhiNoLeadingOppositeHadron->Fill(etaTrig, phiTrig, GetEventWeight());
1387  }
1388  else
1389  {
1390  fhPtLeadingOppositeHadron ->Fill(ptTrig, ptLeadHad, GetEventWeight());
1391  fhPtDiffPhiLeadingOppositeHadron->Fill(ptTrig, dphiLeadHad, GetEventWeight());
1392  fhPtDiffEtaLeadingOppositeHadron->Fill(ptTrig, etaLeadHad-etaTrig, GetEventWeight());
1393  }
1394  }
1395 
1396 
1397  AliDebug(1,Form("pT %2.2f, phi %2.2f, eta %2.2f, nTracks away %d, total tracks %d",
1398  ptLeadHad,phiLeadHad*TMath::RadToDeg(),etaLeadHad,nTrack, GetTrackMultiplicity()));
1399  AliDebug(1,Form("\t pT trig %2.2f, Dphi (trigger-hadron) %2.2f, Deta (trigger-hadron) %2.2f",
1400  ptTrig, dphiLeadHad*TMath::RadToDeg(), etaLeadHad-etaTrig));
1401  AliDebug(1,Form("\t cuts pT: min %2.2f, max %2.2f; DPhi: min %2.2f, max %2.2f",
1402  fMinLeadHadPt,fMaxLeadHadPt,fMinLeadHadPhi*TMath::RadToDeg(),fMaxLeadHadPhi*TMath::RadToDeg()));
1403 
1404 
1405  // reject the trigger if the leading hadron is not in the requested pt or phi window and
1406 
1407  if( nTrack == 0 ) return kFALSE; // No track found in opposite hemisphere
1408 
1409  if( ptLeadHad < fMinLeadHadPt || ptLeadHad > fMaxLeadHadPt ) return kFALSE;
1410 
1411  //printf("Accept leading hadron pT \n");
1412 
1413  if( dphiLeadHad < fMinLeadHadPhi || dphiLeadHad > fMaxLeadHadPhi ) return kFALSE;
1414 
1415  //printf("Accept leading hadron phi \n");
1416 
1417  return kTRUE ;
1418 }
1419 
1420 //____________________________________________________________
1422 //____________________________________________________________
1423 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
1424 {
1425  TString parList ; //this will be list of parameters used for this analysis.
1426  const Int_t buffersize = 560;
1427  char onePar[buffersize] ;
1428 
1429  snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---:") ;
1430  parList+=onePar ;
1431  snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f; ", fMinAssocPt, fMaxAssocPt) ;
1432  parList+=onePar ;
1433  snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f; ", fDeltaPhiMinCut, fDeltaPhiMaxCut) ;
1434  parList+=onePar ;
1435  snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron < %3.2f; ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
1436  parList+=onePar ;
1437  snprintf(onePar,buffersize,"Isolated Trigger? %d;", fSelectIsolated) ;
1438  parList+=onePar ;
1439  snprintf(onePar,buffersize,"Several UE? %d;", fMakeSeveralUE) ;
1440  parList+=onePar ;
1441  snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s;", fPi0AODBranchName.Data());
1442  parList+=onePar ;
1443  snprintf(onePar,buffersize,"Do Decay-hadron correlation ? pi0 %d, decay %d;", fPi0Trigger, fDecayTrigger) ;
1444  parList+=onePar ;
1445  snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d;",
1446  fMakeAbsoluteLeading, fMakeNearSideLeading) ;
1447  parList+=onePar ;
1448  snprintf(onePar,buffersize,"Associated particle pt bins %d: ", fNAssocPtBins) ;
1449  parList+=onePar ;
1450  for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
1451  snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
1452  }
1453  parList+=onePar ;
1454 
1455  //Get parameters set in base class.
1456  parList += GetBaseParametersList() ;
1457 
1458  //Get parameters set in FiducialCut class (not available yet)
1459  //parlist += GetFidCut()->GetFidCutParametersList()
1460 
1461  return new TObjString(parList) ;
1462 }
1463 
1464 //________________________________________________________________
1467 //________________________________________________________________
1468 TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
1469 {
1470  TList * outputContainer = new TList() ;
1471  outputContainer->SetName("CorrelationHistos") ;
1472 
1473  Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
1474  Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
1475  Float_t ptmin = GetHistogramRanges()->GetHistoPtMin(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1476 
1477  Int_t ndeltaphibins = GetHistogramRanges()->GetHistoDeltaPhiBins(); Int_t ndeltaetabins = GetHistogramRanges()->GetHistoDeltaEtaBins();
1478  Float_t deltaphimax = GetHistogramRanges()->GetHistoDeltaPhiMax(); Float_t deltaetamax = GetHistogramRanges()->GetHistoDeltaEtaMax();
1479  Float_t deltaphimin = GetHistogramRanges()->GetHistoDeltaPhiMin(); Float_t deltaetamin = GetHistogramRanges()->GetHistoDeltaEtaMin();
1480 
1481  Int_t ntrbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins();
1482  Int_t trmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t clmax = GetHistogramRanges()->GetHistoNClustersMax();
1483  Int_t trmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin(); Int_t clmin = GetHistogramRanges()->GetHistoNClustersMin();
1484 
1485  Int_t nxeztbins = GetHistogramRanges()->GetHistoRatioBins(); Int_t nhbpbins = GetHistogramRanges()->GetHistoHBPBins();
1486  Float_t xeztmax = GetHistogramRanges()->GetHistoRatioMax(); Float_t hbpmax = GetHistogramRanges()->GetHistoHBPMax();
1487  Float_t xeztmin = GetHistogramRanges()->GetHistoRatioMin(); Float_t hbpmin = GetHistogramRanges()->GetHistoHBPMin();
1488 
1489  Int_t nMixBins = GetNCentrBin()*GetNZvertBin()*GetNRPBin();
1490 
1491  Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
1492  Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
1493  Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
1494 
1495  TString nameMC[] = {"Photon","Pi0","Pi0Decay","Eta","EtaDecay","OtherDecay","Electron","Hadron","Pi0DecayLostPair","EtaDecayLostPair"};
1496  TString mcPartType[] = {"#gamma", "#pi^{0} (merged #gamma)", "#gamma_{#pi decay}","#eta (merged #gamma)" , "#gamma_{#eta decay}", "#gamma_{other decay}",
1497  "e^{#pm}" , "hadrons?" , "#gamma_{#pi decay} lost companion", "#gamma_{#eta decay} lost companion"} ;
1498  TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1499 
1500  TString parTitle = Form("#it{R} = %2.2f",GetIsolationCut()->GetConeSize());
1501 
1502  // For vz dependent histograms, if option ON
1503  Int_t nz = 1 ;
1504  if(fCorrelVzBin) nz = GetNZvertBin();
1505  TString sz = "" ;
1506  TString tz = "" ;
1507 
1508  // Fill histograms for neutral clusters in mixing?
1509  Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
1510  Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
1511 
1512  if( GetReader()->GetDataType() == AliCaloTrackReader::kESD )
1513  {
1514  fhTrackResolution = new TH2F ("hTrackResolution","Track resolution: #sigma_{#it{p}_{T}} vs #it{p}_{T}, away side, ESDs",
1515  nptbins,ptmin,ptmax,600,0,0.3);
1516  fhTrackResolution->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1517  fhTrackResolution->SetYTitle("#sigma_{#it{p}_{T}} / #it{p}_{T} (GeV/#it{c})");
1518  outputContainer->Add(fhTrackResolution);
1519 
1520  fhTrackResolutionUE = new TH2F ("hTrackResolutionUE","Track resolution: #sigma_{#it{p}_{T}} vs #it{p}_{T}, UE, ESDs",
1521  nptbins,ptmin,ptmax,600,0,0.3);
1522  fhTrackResolutionUE->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1523  fhTrackResolutionUE->SetYTitle("#sigma_{#it{p}_{T}} / #it{p}_{T} (GeV/#it{c})");
1524  outputContainer->Add(fhTrackResolutionUE);
1525  }
1526 
1527  fhPtTriggerInput = new TH1F("hPtTriggerInput","Input trigger #it{p}_{T}", nptbins,ptmin,ptmax);
1528  fhPtTriggerInput->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1529  outputContainer->Add(fhPtTriggerInput);
1530 
1531  if( fM02MaxCut > 0 && fM02MinCut > 0 )
1532  {
1533  fhPtTriggerSSCut = new TH1F("hPtTriggerSSCut","Trigger #it{p}_{T} after #lambda^{2}_{0} cut", nptbins,ptmin,ptmax);
1534  fhPtTriggerSSCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1535  outputContainer->Add(fhPtTriggerSSCut);
1536  }
1537 
1538  if( OnlyIsolated() )
1539  {
1540  fhPtTriggerIsoCut = new TH1F("hPtTriggerIsoCut","Trigger #it{p}_{T} after isolation (and #lambda^{2}_{0} cut)", nptbins,ptmin,ptmax);
1541  fhPtTriggerIsoCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1542  outputContainer->Add(fhPtTriggerIsoCut);
1543  }
1544 
1545  fhPtTriggerFidCut = new TH1F("hPtTriggerFidCut","Trigger #it{p}_{T} after fiducial (isolation and #lambda^{2}_{0}) cut", nptbins,ptmin,ptmax);
1546  fhPtTriggerFidCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1547  outputContainer->Add(fhPtTriggerFidCut);
1548 
1549  fhPtTrigger = new TH1F("hPtTrigger","#it{p}_{T} distribution of trigger particles (after opposite hadron leading cut and rest)", nptbins,ptmin,ptmax);
1550  fhPtTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1551  outputContainer->Add(fhPtTrigger);
1552 
1553  if(fFillInvMassHisto)
1554  {
1555  fhMassPtTrigger = new TH2F("hMassPtTrigger","2 photons invariant mass vs p_{T}^{trig}",
1556  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1557  fhMassPtTrigger->SetYTitle("M_{#gamma#gamma} (GeV/#it{c}^{2})");
1558  fhMassPtTrigger->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1559  outputContainer->Add(fhMassPtTrigger);
1560  }
1561 
1562  if(fFillBkgBinsHisto)
1563  {
1564 
1565  fhPtLeadInConeBin = new TH1F*[fNBkgBin] ;
1566 
1567  fhPtSumInConeBin = new TH1F*[fNBkgBin] ;
1568 
1569  fhPtLeadConeBinDecay = new TH1F*[fNBkgBin*fNDecayBits] ;
1570 
1571  fhSumPtConeBinDecay = new TH1F*[fNBkgBin*fNDecayBits] ;
1572 
1573  fhPtLeadConeBinMC = new TH1F*[fNBkgBin*fgkNmcTypes] ;
1574 
1575  fhSumPtConeBinMC = new TH1F*[fNBkgBin*fgkNmcTypes] ;
1576 
1577  for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1578  {
1579  fhPtLeadInConeBin[ibin] = new TH1F
1580  (Form("hPtLeadCone_Bin%d",ibin),
1581  Form("cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1582  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1583  fhPtLeadInConeBin[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
1584  fhPtLeadInConeBin[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1585  outputContainer->Add(fhPtLeadInConeBin[ibin]) ;
1586 
1587  fhPtSumInConeBin[ibin] = new TH1F
1588  (Form("hSumPtCone_Bin%d",ibin),
1589  Form("in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, %s",
1590  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1591  fhPtSumInConeBin[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
1592  fhPtSumInConeBin[ibin]->SetXTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
1593  outputContainer->Add(fhPtSumInConeBin[ibin]) ;
1594 
1595  if(fFillTaggedDecayHistograms)
1596  {
1597  for(Int_t idecay = 0; idecay < fNDecayBits; idecay++)
1598  {
1599  Int_t bindecay = ibin+idecay*fNBkgBin;
1600 
1601  fhPtLeadConeBinDecay[bindecay] = new TH1F
1602  (Form("hPtLeadCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
1603  Form("Decay bit %d, cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1604  fDecayBits[idecay],fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1605  fhPtLeadConeBinDecay[bindecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
1606  fhPtLeadConeBinDecay[bindecay]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1607  outputContainer->Add(fhPtLeadConeBinDecay[bindecay]) ;
1608 
1609  fhSumPtConeBinDecay[bindecay] = new TH1F
1610  (Form("hSumPtCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
1611  Form("Decay bit %d, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, %s",
1612  fDecayBits[idecay],fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1613  fhSumPtConeBinDecay[bindecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
1614  fhSumPtConeBinDecay[bindecay]->SetXTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
1615  outputContainer->Add(fhSumPtConeBinDecay[bindecay]) ;
1616  }
1617  }
1618 
1619  if(IsDataMC())
1620  {
1621  for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1622  {
1623  Int_t binmc = ibin+imc*fNBkgBin;
1624  fhPtLeadConeBinMC[binmc] = new TH1F
1625  (Form("hPtLeadCone_Bin%d_MC%s",ibin, nameMC[imc].Data()),
1626  Form("in cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, MC %s, %s",
1627  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
1628  fhPtLeadConeBinMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
1629  fhPtLeadConeBinMC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1630  outputContainer->Add(fhPtLeadConeBinMC[binmc]) ;
1631 
1632  fhSumPtConeBinMC[binmc] = new TH1F
1633  (Form("hSumPtCone_Bin%d_MC%s",ibin,nameMC[imc].Data()),
1634  Form("in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, MC %s, %s",
1635  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
1636  fhSumPtConeBinMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
1637  fhSumPtConeBinMC[binmc]->SetXTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
1638  outputContainer->Add(fhSumPtConeBinMC[binmc]) ;
1639  } // MC particle loop
1640  }
1641  }
1642  }
1643 
1644  if(IsDataMC())
1645  {
1646  for(Int_t i=0; i < fgkNmcTypes; i++)
1647  {
1648  fhPtTriggerMC[i] = new TH1F(Form("hPtTrigger_MC%s",nameMC[i].Data()),
1649  Form("#it{p}_{T} distribution of trigger particles, trigger origin is %s",nameMC[i].Data()),
1650  nptbins,ptmin,ptmax);
1651  fhPtTriggerMC[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1652  outputContainer->Add(fhPtTriggerMC[i]);
1653 
1654  if(fFillInvMassHisto)
1655  {
1656  fhMCMassPtTrigger[i] = new TH2F(Form("hMassPtTrigger_MC%s",nameMC[i].Data()),
1657  Form("2 photons invariant mass, trigger origin is %s",nameMC[i].Data()),
1658  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1659  fhMCMassPtTrigger[i]->SetYTitle("m_{#gamma#gamma}");
1660  fhMCMassPtTrigger[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1661  outputContainer->Add(fhMCMassPtTrigger[i]) ;
1662  }
1663  }
1664  }
1665 
1666  if(fDecayTrigger)
1667  {
1668  for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1669  {
1670  fhPtDecayTrigger[ibit] = new TH1F(Form("hPtDecayTrigger_bit%d",fDecayBits[ibit]),
1671  Form("#it{p}_{T} distribution of trigger particles, decay Bit %d",fDecayBits[ibit]),
1672  nptbins,ptmin,ptmax);
1673  fhPtDecayTrigger[ibit]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1674  outputContainer->Add(fhPtDecayTrigger[ibit]);
1675 
1676  if(IsDataMC())
1677  {
1678  for(Int_t i=0; i < fgkNmcTypes; i++)
1679  {
1680  fhPtDecayTriggerMC[ibit][i] = new TH1F(Form("hPtDecayTrigger_bit%d_MC%s",fDecayBits[ibit], nameMC[i].Data()),
1681  Form("#it{p}_{T} distribution of trigger particles, decay Bit %d, trigger origin is %s",fDecayBits[ibit], nameMC[i].Data()),
1682  nptbins,ptmin,ptmax);
1683  fhPtDecayTriggerMC[ibit][i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1684  outputContainer->Add(fhPtDecayTriggerMC[ibit][i]);
1685  }
1686  }
1687  }
1688  }
1689 
1690  if(fCorrelVzBin)
1691  {
1692  fhPtTriggerVzBin = new TH2F("hPtTriggerVzBin","#it{p}_{T} distribution of trigger particles vs vz bin", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
1693  fhPtTriggerVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1694  fhPtTriggerVzBin->SetYTitle("#it{v}_{#it{z}} bin");
1695  outputContainer->Add(fhPtTriggerVzBin);
1696  }
1697 
1698  fhPtTriggerBin = new TH2F ("hPtTriggerBin","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
1699  fhPtTriggerBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1700  fhPtTriggerBin->SetYTitle("Bin");
1701  outputContainer->Add(fhPtTriggerBin);
1702 
1703  fhPhiTrigger = new TH2F ("hPhiTrigger","#phi distribution of trigger Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1704  fhPhiTrigger->SetYTitle("#phi (rad)");
1705  outputContainer->Add(fhPhiTrigger);
1706 
1707  fhEtaTrigger = new TH2F ("hEtaTrigger","#eta distribution of trigger",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1708  fhEtaTrigger->SetYTitle("#eta ");
1709  outputContainer->Add(fhEtaTrigger);
1710 
1711  if(IsHighMultiplicityAnalysisOn())
1712  {
1713  fhPtTriggerCentrality = new TH2F("hPtTriggerCentrality","Trigger particle #it{p}_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
1714  fhPtTriggerCentrality->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1715  fhPtTriggerCentrality->SetYTitle("Centrality (%)");
1716  outputContainer->Add(fhPtTriggerCentrality) ;
1717 
1718  fhPtTriggerEventPlane = new TH2F("hPtTriggerEventPlane","Trigger particle #it{p}_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
1719  fhPtTriggerEventPlane->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1720  fhPtTriggerEventPlane->SetXTitle("EP angle (rad)");
1721  outputContainer->Add(fhPtTriggerEventPlane) ;
1722 
1723  fhTriggerEventPlaneCentrality = new TH2F("hTriggerEventPlaneCentrality","Trigger particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
1724  fhTriggerEventPlaneCentrality->SetXTitle("Centrality (%)");
1725  fhTriggerEventPlaneCentrality->SetYTitle("EP angle (rad)");
1726  outputContainer->Add(fhTriggerEventPlaneCentrality) ;
1727  }
1728 
1729  // Leading hadron in oposite side
1730  if(fFillLeadHadOppositeHisto)
1731  {
1732  fhPtLeadingOppositeHadron = new TH2F("hPtTriggerPtLeadingOppositeHadron","Leading hadron opposite to trigger vs trigger #it{p}_{T}",
1733  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1734  fhPtLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1735  fhPtLeadingOppositeHadron->SetYTitle("#it{p}_{T}^{lead hadron} (GeV/#it{c})");
1736  outputContainer->Add(fhPtLeadingOppositeHadron);
1737 
1738  fhPtNoLeadingOppositeHadron = new TH1F("hPtTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #it{p}_{T}",
1739  nptbins,ptmin,ptmax);
1740  fhPtNoLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1741  outputContainer->Add(fhPtNoLeadingOppositeHadron);
1742 
1743  fhEtaPhiNoLeadingOppositeHadron = new TH2F("hEtaPhiTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #eta:#phi",
1744  netabins,etamin,etamax,nphibins,phimin,phimax);
1745  fhEtaPhiNoLeadingOppositeHadron->SetXTitle("#eta");
1746  fhEtaPhiNoLeadingOppositeHadron->SetYTitle("#phi");
1747  outputContainer->Add(fhEtaPhiNoLeadingOppositeHadron);
1748 
1749 
1750  fhPtDiffPhiLeadingOppositeHadron = new TH2F("hPtTriggerDiffPhiTriggerLeadingOppositeHadron","#phi_{trigger}-#phi_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1751  nptbins,ptmin,ptmax,ndeltaphibins,deltaphimin,deltaphimax);
1752  fhPtDiffPhiLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1753  fhPtDiffPhiLeadingOppositeHadron->SetYTitle("#phi_{trigger}-#phi_{leading opposite hadron} (rad)");
1754  outputContainer->Add(fhPtDiffPhiLeadingOppositeHadron);
1755 
1756  fhPtDiffEtaLeadingOppositeHadron = new TH2F("hPtTriggerDiffEtaTriggerPhiLeadingOppositeHadron","#eta_{trigger}-#eta_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1757  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1758  fhPtDiffEtaLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1759  fhPtDiffEtaLeadingOppositeHadron->SetYTitle("#eta_{trigger}-#eta_{leading opposite hadron}");
1760  outputContainer->Add(fhPtDiffEtaLeadingOppositeHadron);
1761  }
1762 
1763  // Correlation with charged hadrons
1764 
1765  fhDeltaPhiDeltaEtaCharged = new TH2F
1766  ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
1767  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1768  fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
1769  fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1770 
1771  fhDeltaPhiDeltaEtaChargedPtA3GeV = new TH2F
1772  ("hDeltaPhiDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}, #it{p}_{TA}>3 GeV/#it{c}}",
1773  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1774  fhDeltaPhiDeltaEtaChargedPtA3GeV->SetXTitle("#Delta #phi (rad)");
1775  fhDeltaPhiDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1776 
1777  fhPhiCharged = new TH2F
1778  ("hPhiCharged","#phi_{h^{#pm}} vs #it{p}_{T #pm}",
1779  nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
1780  fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
1781  fhPhiCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1782 
1783  fhEtaCharged = new TH2F
1784  ("hEtaCharged","#eta_{h^{#pm}} vs #it{p}_{T #pm}",
1785  nptbins,ptmin,ptmax,100,-1.,1.);
1786  fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1787  fhEtaCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1788 
1789  fhDeltaPhiCharged = new TH2F
1790  ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
1791  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1792  fhDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
1793  fhDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1794 
1795  fhDeltaPhiChargedPtA3GeV = new TH2F
1796  ("hDeltaPhiChargedPtA3GeV","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1797  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1798  fhDeltaPhiChargedPtA3GeV->SetYTitle("#Delta #phi (rad)");
1799  fhDeltaPhiChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1800 
1801 
1802  fhDeltaPhiChargedPt = new TH2F
1803  ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",
1804  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1805  fhDeltaPhiChargedPt->SetYTitle("#Delta #phi (rad)");
1806  fhDeltaPhiChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1807 
1808  fhDeltaEtaCharged = new TH2F
1809  ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}",
1810  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1811  fhDeltaEtaCharged->SetYTitle("#Delta #eta");
1812  fhDeltaEtaCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1813 
1814  fhDeltaEtaChargedPtA3GeV = new TH2F
1815  ("hDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1816  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1817  fhDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1818  fhDeltaEtaChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1819 
1820  fhXECharged =
1821  new TH2F("hXECharged","#it{x}_{#it{E}} for charged tracks",
1822  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1823  fhXECharged->SetYTitle("#it{x}_{#it{E}}");
1824  fhXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1825 
1826  fhXECharged_Cone2 =
1827  new TH2F("hXECharged_Cone2","#it{x}_{#it{E}} for charged tracks in cone 2 (5#pi/6-7#pi/6)",
1828  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1829  fhXECharged_Cone2->SetYTitle("#it{x}_{#it{E}}");
1830  fhXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1831 
1832  fhXEPosCharged =
1833  new TH2F("hXEPositiveCharged","#it{x}_{#it{E}} for positive charged tracks",
1834  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1835  fhXEPosCharged->SetYTitle("#it{x}_{#it{E}}");
1836  fhXEPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1837 
1838  fhXENegCharged =
1839  new TH2F("hXENegativeCharged","#it{x}_{#it{E}} for negative charged tracks",
1840  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1841  fhXENegCharged->SetYTitle("#it{x}_{#it{E}}");
1842  fhXENegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1843 
1844  fhPtHbpXECharged =
1845  new TH2F("hHbpXECharged","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
1846  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1847  fhPtHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1848  fhPtHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1849 
1850  fhPtHbpXECharged_Cone2 =
1851  new TH2F("hHbpXECharged_Cone2","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons in cone 2 (5#pi/6-7#pi/6)",
1852  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1853  fhPtHbpXECharged_Cone2->SetYTitle("ln(1/#it{x}_{#it{E}})");
1854  fhPtHbpXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1855 
1856  fhZTCharged =
1857  new TH2F("hZTCharged","#it{z}_{T} for charged tracks",
1858  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1859  fhZTCharged->SetYTitle("#it{z}_{T}");
1860  fhZTCharged->SetXTitle("#it{p}_{T trigger}");
1861 
1862  fhZTPosCharged =
1863  new TH2F("hZTPositiveCharged","#it{z}_{T} for positive charged tracks",
1864  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1865  fhZTPosCharged->SetYTitle("#it{z}_{T}");
1866  fhZTPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1867 
1868  fhZTNegCharged =
1869  new TH2F("hZTNegativeCharged","#it{z}_{T} for negative charged tracks",
1870  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1871  fhZTNegCharged->SetYTitle("#it{z}_{T}");
1872  fhZTNegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1873 
1874  fhPtHbpZTCharged =
1875  new TH2F("hHbpZTCharged","#xi = ln(1/#it{z}_{T}) with charged hadrons",
1876  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1877  fhPtHbpZTCharged->SetYTitle("ln(1/#it{z}_{T})");
1878  fhPtHbpZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1879 
1880  fhPtTrigPout =
1881  new TH2F("hPtTrigPout","Pout with triggers",
1882  nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
1883  fhPtTrigPout->SetYTitle("#it{p}_{out} (GeV/#it{c})");
1884  fhPtTrigPout->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1885 
1886  fhPtTrigCharged =
1887  new TH2F("hPtTrigCharged","trigger and charged tracks pt distribution",
1888  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1889  fhPtTrigCharged->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1890  fhPtTrigCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1891 
1892  outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
1893  outputContainer->Add(fhDeltaPhiDeltaEtaChargedPtA3GeV);
1894  outputContainer->Add(fhPhiCharged) ;
1895  outputContainer->Add(fhEtaCharged) ;
1896  outputContainer->Add(fhDeltaPhiCharged) ;
1897  outputContainer->Add(fhDeltaPhiChargedPtA3GeV) ;
1898  outputContainer->Add(fhDeltaEtaCharged) ;
1899  outputContainer->Add(fhDeltaEtaChargedPtA3GeV) ;
1900  outputContainer->Add(fhDeltaPhiChargedPt) ;
1901 
1902  outputContainer->Add(fhXECharged) ;
1903  outputContainer->Add(fhXECharged_Cone2) ;
1904 
1905  if(IsDataMC())
1906  {
1907  for(Int_t i=0; i < fgkNmcTypes; i++)
1908  {
1909 
1910  fhDeltaPhiChargedMC[i] = new TH2F(Form("hDeltaPhiCharged_MC%s",nameMC[i].Data()),
1911  Form("#Delta #phi for charged tracks, trigger origin is %s",nameMC[i].Data()),
1912  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1913  fhDeltaPhiChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1914  fhDeltaPhiChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1915  outputContainer->Add(fhDeltaPhiChargedMC[i]) ;
1916 
1917  fhXEChargedMC[i] = new TH2F(Form("hXECharged_MC%s",nameMC[i].Data()),
1918  Form("#it{x}_{#it{E}} for charged tracks, trigger origin is %s",nameMC[i].Data()),
1919  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1920  fhXEChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1921  fhXEChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1922  outputContainer->Add(fhXEChargedMC[i]) ;
1923 
1924  fhXEUeChargedRightMC[i] = new TH2F(Form("hXEUeChargedRight_MC%s",nameMC[i].Data()),
1925  Form("#it{x}_{#it{E}} for charged tracks in right UE cone, trigger origin is %s",nameMC[i].Data()),
1926  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1927  fhXEUeChargedRightMC[i]->SetYTitle("#it{x}_{#it{E}}");
1928  fhXEUeChargedRightMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1929  outputContainer->Add(fhXEUeChargedRightMC[i]) ;
1930 
1931  fhXEUeChargedLeftMC[i] = new TH2F(Form("hXEUeChargedLeft_MC%s",nameMC[i].Data()),
1932  Form("#it{x}_{#it{E}} for charged tracks in left UE cone, trigger origin is %s",nameMC[i].Data()),
1933  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1934  fhXEUeChargedLeftMC[i]->SetYTitle("#it{x}_{#it{E}}");
1935  fhXEUeChargedLeftMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1936  outputContainer->Add(fhXEUeChargedLeftMC[i]) ;
1937  }
1938  }
1939 
1940  outputContainer->Add(fhXEPosCharged) ;
1941  outputContainer->Add(fhXENegCharged) ;
1942  outputContainer->Add(fhPtHbpXECharged) ;
1943  outputContainer->Add(fhPtHbpXECharged_Cone2) ;
1944 
1945  outputContainer->Add(fhZTCharged) ;
1946  outputContainer->Add(fhZTPosCharged) ;
1947  outputContainer->Add(fhZTNegCharged) ;
1948  outputContainer->Add(fhPtHbpZTCharged) ;
1949 
1950  outputContainer->Add(fhPtTrigPout) ;
1951  outputContainer->Add(fhPtTrigCharged) ;
1952 
1953  TString right = "";
1954  if(fMakeSeveralUE) right = "Right";
1955 
1956  fhUePart = new TH1F("hUePart","UE particles distribution vs pt trig",
1957  nptbins,ptmin,ptmax);
1958  fhUePart->SetYTitle("dNch");
1959  fhUePart->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1960 
1961  fhDeltaPhiUeChargedPt = new TH2F
1962  (Form("hDeltaPhiUe%sChargedPt",right.Data()),"#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}}",
1963  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1964  fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi (rad)");
1965  fhDeltaPhiUeChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1966 
1967  fhXEUeCharged =
1968  new TH2F(Form("hXEUeCharged%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event",
1969  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1970  fhXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
1971  fhXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1972 
1973  fhXEUeChargedSmallCone =
1974  new TH2F(Form("hXEUeChargedSmallCone%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event in a cone [80,100] deg",
1975  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1976  fhXEUeChargedSmallCone->SetYTitle("#it{x}_{#it{E}}");
1977  fhXEUeChargedSmallCone->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1978 
1979  fhXEUeChargedMediumCone =
1980  new TH2F(Form("hXEUeChargedMediumCone%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event in a cone [70,110] deg",
1981  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1982  fhXEUeChargedMediumCone->SetYTitle("#it{x}_{#it{E}}");
1983  fhXEUeChargedMediumCone->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1984 
1985  fhXEUeChargedLargeCone =
1986  new TH2F(Form("hXEUeChargedLargeCone%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event in a cone [60,120] deg",
1987  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1988  fhXEUeChargedLargeCone->SetYTitle("#it{x}_{#it{E}}");
1989  fhXEUeChargedLargeCone->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1990 
1991  fhPtHbpXEUeCharged =
1992  new TH2F(Form("hHbpXEUeCharged%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for Underlying Event",
1993  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1994  fhPtHbpXEUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1995  fhPtHbpXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1996 
1997  fhZTUeCharged =
1998  new TH2F(Form("hZTUeCharged%s",right.Data()),"#it{z}_{T} for Underlying Event",
1999  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2000  fhZTUeCharged->SetYTitle("#it{z}_{T}");
2001  fhZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2002 
2003  fhPtHbpZTUeCharged =
2004  new TH2F(Form("hHbpZTUeCharged%s",right.Data()),"#xi = ln(1/#it{z}_{T}) for Underlying Event",
2005  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2006  fhPtHbpZTUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2007  fhPtHbpZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2008 
2009  outputContainer->Add(fhUePart);
2010  outputContainer->Add(fhDeltaPhiUeChargedPt) ;
2011  outputContainer->Add(fhXEUeCharged) ;
2012  outputContainer->Add(fhXEUeChargedSmallCone) ;
2013  outputContainer->Add(fhXEUeChargedMediumCone) ;
2014  outputContainer->Add(fhXEUeChargedLargeCone) ;
2015  outputContainer->Add(fhPtHbpXEUeCharged) ;
2016  outputContainer->Add(fhZTUeCharged) ;
2017  outputContainer->Add(fhPtHbpZTUeCharged) ;
2018 
2019  if(fMakeSeveralUE)
2020  {
2021  fhDeltaPhiUeLeftCharged = new TH2F
2022  ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left side range of trigger particles",
2023  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2024  fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi (rad)");
2025  fhDeltaPhiUeLeftCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2026  outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
2027 
2028  fhDeltaPhiUeLeftUpCharged = new TH2F
2029  ("hDeltaPhiUeLeftUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Up side range of trigger particles",
2030  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2031  fhDeltaPhiUeLeftUpCharged->SetYTitle("#Delta #phi (rad)");
2032  fhDeltaPhiUeLeftUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2033  outputContainer->Add(fhDeltaPhiUeLeftUpCharged) ;
2034 
2035  fhDeltaPhiUeRightUpCharged = new TH2F
2036  ("hDeltaPhiUeRightUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Up side range of trigger particles",
2037  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2038  fhDeltaPhiUeRightUpCharged->SetYTitle("#Delta #phi (rad)");
2039  fhDeltaPhiUeRightUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2040  outputContainer->Add(fhDeltaPhiUeRightUpCharged) ;
2041 
2042  fhDeltaPhiUeLeftDownCharged = new TH2F
2043  ("hDeltaPhiUeLeftDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Down side range of trigger particles",
2044  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2045  fhDeltaPhiUeLeftDownCharged->SetYTitle("#Delta #phi (rad)");
2046  fhDeltaPhiUeLeftDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2047  outputContainer->Add(fhDeltaPhiUeLeftDownCharged) ;
2048 
2049  fhDeltaPhiUeRightDownCharged = new TH2F
2050  ("hDeltaPhiUeRightDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Down side range of trigger particles",
2051  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2052  fhDeltaPhiUeRightDownCharged->SetYTitle("#Delta #phi (rad)");
2053  fhDeltaPhiUeRightDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2054  outputContainer->Add(fhDeltaPhiUeRightDownCharged) ;
2055 
2056  fhXEUeLeftCharged =
2057  new TH2F("hXEUeChargedLeft","#it{x}_{#it{E}} with UE left side of trigger",
2058  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2059  fhXEUeLeftCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
2060  fhXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2061  outputContainer->Add(fhXEUeLeftCharged) ;
2062 
2063  fhXEUeLeftUpCharged =
2064  new TH2F("hXEUeChargedLeftUp","#it{x}_{#it{E}} with UE left Up side of trigger",
2065  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2066  fhXEUeLeftUpCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
2067  fhXEUeLeftUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2068  outputContainer->Add(fhXEUeLeftUpCharged) ;
2069 
2070  fhXEUeRightUpCharged =
2071  new TH2F("hXEUeChargedRightUp","#it{x}_{#it{E} h^{#pm}} with UE right Up side of trigger",
2072  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2073  fhXEUeRightUpCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
2074  fhXEUeRightUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2075  outputContainer->Add(fhXEUeRightUpCharged) ;
2076 
2077  fhXEUeLeftDownCharged =
2078  new TH2F("hXEUeChargedLeftDown","#it{x}_{#it{E}} with UE left Down side of trigger",
2079  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2080  fhXEUeLeftDownCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
2081  fhXEUeLeftDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2082  outputContainer->Add(fhXEUeLeftDownCharged) ;
2083 
2084  fhXEUeRightDownCharged =
2085  new TH2F("hXEUeChargedRightDown","#it{x}_{#it{E} h^{#pm}} with UE right Down side of trigger",
2086  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2087  fhXEUeRightDownCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
2088  fhXEUeRightDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2089  outputContainer->Add(fhXEUeRightDownCharged) ;
2090 
2091  fhPtHbpXEUeLeftCharged =
2092  new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/#it{x}_{#it{E}}) with charged UE left side of trigger",
2093  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2094  fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2095  fhPtHbpXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2096  outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
2097 
2098  fhZTUeLeftCharged =
2099  new TH2F("hZTUeChargedLeft","#it{z}_{trigger h^{#pm}} = #it{p}_{T Ueh^{#pm}} / #it{p}_{T trigger} with UE left side of trigger",
2100  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2101  fhZTUeLeftCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
2102  fhZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2103  outputContainer->Add(fhZTUeLeftCharged) ;
2104 
2105  fhPtHbpZTUeLeftCharged =
2106  new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/#it{z}_{T}) with charged UE left side of trigger",
2107  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2108  fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/#it{z}_{T})");
2109  fhPtHbpZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2110  outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
2111  }
2112 
2113  if(IsPileUpAnalysisOn())
2114  {
2115  fhDeltaPhiChargedOtherBC = new TH2F
2116  ("hDeltaPhiChargedOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC!=0",
2117  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2118  fhDeltaPhiChargedOtherBC->SetYTitle("#Delta #phi (rad)");
2119  fhDeltaPhiChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2120 
2121  fhDeltaPhiChargedPtA3GeVOtherBC = new TH2F
2122  ("hDeltaPhiChargedPtA3GeVOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC!=0",
2123  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2124  fhDeltaPhiChargedPtA3GeVOtherBC->SetYTitle("#Delta #phi (rad)");
2125  fhDeltaPhiChargedPtA3GeVOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2126 
2127  fhPtTrigChargedOtherBC =
2128  new TH2F("hPtTrigChargedOtherBC","trigger and charged tracks pt distribution, track BC!=0",
2129  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2130  fhPtTrigChargedOtherBC->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2131  fhPtTrigChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2132 
2133  fhXEChargedOtherBC =
2134  new TH2F("hXEChargedOtherBC","#it{x}_{#it{E}} for charged tracks, track BC!=0",
2135  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2136  fhXEChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
2137  fhXEChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2138 
2139  fhXEUeChargedOtherBC =
2140  new TH2F("hXEUeChargedOtherBC","#it{x}_{#it{E}} for Underlying Event, track BC!=0",
2141  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2142  fhXEUeChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
2143  fhXEUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2144 
2145  fhZTChargedOtherBC =
2146  new TH2F("hZTChargedOtherBC","#it{z}_{T} for charged tracks, track BC!=0",
2147  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2148  fhZTChargedOtherBC->SetYTitle("#it{z}_{T}");
2149  fhZTChargedOtherBC->SetXTitle("#it{p}_{T trigger}");
2150 
2151  fhZTUeChargedOtherBC =
2152  new TH2F("hZTUeChargedOtherBC","#it{z}_{T} for Underlying Event, track BC!=0",
2153  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2154  fhZTUeChargedOtherBC->SetYTitle("#it{z}_{T}");
2155  fhZTUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2156 
2157  outputContainer->Add(fhDeltaPhiChargedOtherBC) ;
2158  outputContainer->Add(fhDeltaPhiChargedPtA3GeVOtherBC) ;
2159  outputContainer->Add(fhXEChargedOtherBC) ;
2160  outputContainer->Add(fhXEUeChargedOtherBC) ;
2161  outputContainer->Add(fhZTChargedOtherBC) ;
2162  outputContainer->Add(fhZTUeChargedOtherBC) ;
2163  outputContainer->Add(fhPtTrigChargedOtherBC) ;
2164 
2165  fhDeltaPhiChargedBC0 = new TH2F
2166  ("hDeltaPhiChargedBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
2167  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2168  fhDeltaPhiChargedBC0->SetYTitle("#Delta #phi (rad)");
2169  fhDeltaPhiChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2170 
2171  fhDeltaPhiChargedPtA3GeVBC0 = new TH2F
2172  ("hDeltaPhiChargedPtA3GeVBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
2173  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2174  fhDeltaPhiChargedPtA3GeVBC0->SetYTitle("#Delta #phi (rad)");
2175  fhDeltaPhiChargedPtA3GeVBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2176 
2177  fhPtTrigChargedBC0 =
2178  new TH2F("hPtTrigChargedBC0","trigger and charged tracks pt distribution, track BC==0",
2179  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2180  fhPtTrigChargedBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2181  fhPtTrigChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2182 
2183  fhXEChargedBC0 =
2184  new TH2F("hXEChargedBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
2185  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2186  fhXEChargedBC0->SetYTitle("#it{x}_{#it{E}}");
2187  fhXEChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2188 
2189  fhXEUeChargedBC0 =
2190  new TH2F("hXEUeChargedBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
2191  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2192  fhXEUeChargedBC0->SetYTitle("#it{x}_{#it{E}}");
2193  fhXEUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2194 
2195  fhZTChargedBC0 =
2196  new TH2F("hZTChargedBC0","#it{z}_{T} for charged tracks, track BC==0",
2197  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2198  fhZTChargedBC0->SetYTitle("#it{z}_{T}");
2199  fhZTChargedBC0->SetXTitle("#it{p}_{T trigger}");
2200 
2201  fhZTUeChargedBC0 =
2202  new TH2F("hZTUeChargedBC0","#it{z}_{T} for Underlying Event, track BC==0",
2203  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2204  fhZTUeChargedBC0->SetYTitle("#it{z}_{T}");
2205  fhZTUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2206 
2207  outputContainer->Add(fhDeltaPhiChargedBC0) ;
2208  outputContainer->Add(fhDeltaPhiChargedPtA3GeVBC0) ;
2209  outputContainer->Add(fhXEChargedBC0) ;
2210  outputContainer->Add(fhXEUeChargedBC0) ;
2211  outputContainer->Add(fhZTChargedBC0) ;
2212  outputContainer->Add(fhZTUeChargedBC0) ;
2213  outputContainer->Add(fhPtTrigChargedBC0) ;
2214 
2215  fhPtTriggerVtxBC0 = new TH1F("hPtTriggerVtxBC0","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax);
2216  fhPtTriggerVtxBC0->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2217 
2218  fhDeltaPhiChargedVtxBC0 = new TH2F
2219  ("hDeltaPhiChargedVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
2220  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2221  fhDeltaPhiChargedVtxBC0->SetYTitle("#Delta #phi (rad)");
2222  fhDeltaPhiChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2223 
2224  fhDeltaPhiChargedPtA3GeVVtxBC0 = new TH2F
2225  ("hDeltaPhiChargedPtA3GeVVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
2226  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2227  fhDeltaPhiChargedPtA3GeVVtxBC0->SetYTitle("#Delta #phi (rad)");
2228  fhDeltaPhiChargedPtA3GeVVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2229 
2230  fhPtTrigChargedVtxBC0 =
2231  new TH2F("hPtTrigChargedVtxBC0","trigger and charged tracks pt distribution, track BC==0",
2232  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2233  fhPtTrigChargedVtxBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2234  fhPtTrigChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2235 
2236  fhXEChargedVtxBC0 =
2237  new TH2F("hXEChargedVtxBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
2238  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2239  fhXEChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
2240  fhXEChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2241 
2242  fhXEUeChargedVtxBC0 =
2243  new TH2F("hXEUeChargedVtxBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
2244  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2245  fhXEUeChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
2246  fhXEUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2247 
2248  fhZTChargedVtxBC0 =
2249  new TH2F("hZTChargedVtxBC0","#it{z}_{T} for charged tracks, track BC==0",
2250  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2251  fhZTChargedVtxBC0->SetYTitle("#it{z}_{T}");
2252  fhZTChargedVtxBC0->SetXTitle("#it{p}_{T trigger}");
2253 
2254  fhZTUeChargedVtxBC0 =
2255  new TH2F("hZTUeChargedVtxBC0","#it{z}_{T} for Underlying Event, track BC==0",
2256  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2257  fhZTUeChargedVtxBC0->SetYTitle("#it{z}_{T}");
2258  fhZTUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2259 
2260  outputContainer->Add(fhPtTriggerVtxBC0);
2261  outputContainer->Add(fhDeltaPhiChargedVtxBC0) ;
2262  outputContainer->Add(fhDeltaPhiChargedPtA3GeVVtxBC0) ;
2263  outputContainer->Add(fhXEChargedVtxBC0) ;
2264  outputContainer->Add(fhXEUeChargedVtxBC0) ;
2265  outputContainer->Add(fhZTChargedVtxBC0) ;
2266  outputContainer->Add(fhZTUeChargedVtxBC0) ;
2267  outputContainer->Add(fhPtTrigChargedVtxBC0) ;
2268 
2269  for(Int_t i = 0 ; i < 7 ; i++)
2270  {
2271  fhPtTriggerPileUp[i] = new TH1F(Form("hPtTriggerPileUp%s",pileUpName[i].Data()),
2272  Form("#it{p}_{T} distribution of trigger particles, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2273  fhPtTriggerPileUp[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2274  outputContainer->Add(fhPtTriggerPileUp[i]);
2275 
2276  fhDeltaPhiChargedPileUp[i] = new TH2F(Form("hDeltaPhiChargedPileUp%s",pileUpName[i].Data()),
2277  Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
2278  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2279  fhDeltaPhiChargedPileUp[i]->SetYTitle("#Delta #phi (rad)");
2280  fhDeltaPhiChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2281  outputContainer->Add(fhDeltaPhiChargedPileUp[i]) ;
2282 
2283  fhDeltaPhiChargedPtA3GeVPileUp[i] = new TH2F(Form("hDeltaPhiChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
2284  Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
2285  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2286  fhDeltaPhiChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #phi (rad)");
2287  fhDeltaPhiChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2288  outputContainer->Add(fhDeltaPhiChargedPtA3GeVPileUp[i]) ;
2289 
2290  fhDeltaEtaChargedPileUp[i] = new TH2F(Form("hDeltaEtaChargedPileUp%s",pileUpName[i].Data()),
2291  Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
2292  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
2293  fhDeltaEtaChargedPileUp[i]->SetYTitle("#Delta #eta");
2294  fhDeltaEtaChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2295  outputContainer->Add(fhDeltaEtaChargedPileUp[i]) ;
2296 
2297  fhDeltaEtaChargedPtA3GeVPileUp[i] = new TH2F(Form("hDeltaEtaChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
2298  Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
2299  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
2300  fhDeltaEtaChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #eta");
2301  fhDeltaEtaChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2302  outputContainer->Add(fhDeltaEtaChargedPtA3GeVPileUp[i]) ;
2303 
2304  fhXEChargedPileUp[i] = new TH2F(Form("hXEChargedPileUp%s",pileUpName[i].Data()),
2305  Form("#it{x}_{#it{E}} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
2306  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2307  fhXEChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
2308  fhXEChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2309  outputContainer->Add(fhXEChargedPileUp[i]) ;
2310 
2311  fhXEUeChargedPileUp[i] = new TH2F(Form("hXEUeChargedPileUp%s",pileUpName[i].Data()),
2312  Form("#it{x}_{#it{E}} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
2313  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2314  fhXEUeChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
2315  fhXEUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2316  outputContainer->Add(fhXEUeChargedPileUp[i]) ;
2317 
2318  fhZTChargedPileUp[i] = new TH2F(Form("hZTChargedPileUp%s",pileUpName[i].Data()),
2319  Form("#it{z}_{T} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
2320  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2321  fhZTChargedPileUp[i]->SetYTitle("#it{z}_{T}");
2322  fhZTChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2323  outputContainer->Add(fhZTChargedPileUp[i]) ;
2324 
2325  fhZTUeChargedPileUp[i] = new TH2F(Form("hZTUeChargedPileUp%s",pileUpName[i].Data()),
2326  Form("#it{z}_{T} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
2327  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2328  fhZTUeChargedPileUp[i]->SetYTitle("#it{z}_{T}");
2329  fhZTUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2330  outputContainer->Add(fhZTUeChargedPileUp[i]) ;
2331 
2332  fhPtTrigChargedPileUp[i] = new TH2F(Form("hPtTrigChargedPileUp%s",pileUpName[i].Data()),
2333  Form("trigger and charged tracks pt distribution, %s Pile-Up event",pileUpName[i].Data()),
2334  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2335  fhPtTrigChargedPileUp[i]->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2336  fhPtTrigChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2337  outputContainer->Add(fhPtTrigChargedPileUp[i]) ;
2338  }
2339  }
2340 
2341  if(IsHighMultiplicityAnalysisOn())
2342  {
2343  Int_t nMultiBins = GetNCentrBin();
2344  fhDeltaPhiChargedMult = new TH2F*[nMultiBins] ;
2345  fhDeltaEtaChargedMult = new TH2F*[nMultiBins] ;
2346  fhXEMult = new TH2F*[nMultiBins] ;
2347  fhXEUeMult = new TH2F*[nMultiBins] ;
2348  fhZTMult = new TH2F*[nMultiBins] ;
2349  fhZTUeMult = new TH2F*[nMultiBins] ;
2350 
2351  for(Int_t im=0; im<nMultiBins; im++)
2352  {
2353  fhDeltaPhiChargedMult[im] = new TH2F
2354  (Form("hDeltaPhiCharged_Mult%d",im),Form("#Delta #phi charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2355  fhDeltaPhiChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2356  fhDeltaPhiChargedMult[im]->SetYTitle("#Delta #phi (rad)");
2357 
2358  fhDeltaEtaChargedMult[im] = new TH2F
2359  (Form("hDeltaEtaCharged_Mult%d",im),Form("#Delta #eta charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2360  fhDeltaEtaChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2361  fhDeltaEtaChargedMult[im]->SetYTitle("#Delta #eta");
2362 
2363  fhXEMult[im] = new TH2F
2364  (Form("hXECharged_Mult%d",im),Form("#it{x}_{E} charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2365  fhXEMult[im]->SetYTitle("#it{x}_{E}");
2366  fhXEMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2367 
2368  fhXEUeMult[im] = new TH2F
2369  (Form("hXEUeCharged_Mult%d",im),Form("#it{x}_{E} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2370  fhXEUeMult[im]->SetYTitle("#it{x}_{E}");
2371  fhXEUeMult[im]->SetXTitle("#it{p}_{T trigger}(GeV/#it{c})");
2372 
2373  fhZTMult[im] = new TH2F
2374  (Form("hZTCharged_Mult%d",im),Form("#it{z}_{T} charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2375  fhZTMult[im]->SetYTitle("#it{z}_{T}");
2376  fhZTMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2377 
2378  fhZTUeMult[im] = new TH2F
2379  (Form("hZTUeCharged_Mult%d",im),Form("#it{z}_{T} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2380  fhZTUeMult[im]->SetYTitle("#it{z}_{T}");
2381  fhZTUeMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2382 
2383  outputContainer->Add(fhDeltaPhiChargedMult[im]) ;
2384  outputContainer->Add(fhDeltaEtaChargedMult[im]) ;
2385  outputContainer->Add(fhXEMult [im]);
2386  outputContainer->Add(fhXEUeMult[im]);
2387  outputContainer->Add(fhZTMult [im]);
2388  outputContainer->Add(fhZTUeMult[im]);
2389  }
2390  }
2391 
2392  if(fFillBradHisto)
2393  {
2394  fhAssocPtBkg = new TH2F("hAssocPtBkg", " Trigger #it{p}_{T} vs associated hadron #it{p}_{T} from background",
2395  nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
2396  fhAssocPtBkg->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2397  fhAssocPtBkg->SetYTitle("#it{p}_{T associated} (GeV/#it{c})");
2398  outputContainer->Add(fhAssocPtBkg) ;
2399 
2400  fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} ",
2401  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2402  fhDeltaPhiBrad->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2403  fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2404  outputContainer->Add(fhDeltaPhiBrad) ;
2405  }
2406 
2407  fhDeltaPhiDeltaEtaAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2408  fhDeltaPhiAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2409  fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
2410  fhDeltaPhiAssocPtBinDEta0 = new TH2F*[fNAssocPtBins*nz];
2411  if(fFillMomImbalancePtAssocBinsHisto)
2412  {
2413  fhXEAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2414  fhZTAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2415  }
2416 
2417  if(fCorrelVzBin)
2418  {
2419  fhXEVZ = new TH2F*[nz];
2420  fhZTVZ = new TH2F*[nz];
2421  }
2422 
2423  if(fFillBradHisto)
2424  fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2425 
2426 
2427  fhDeltaPhiAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2428  if(fFillEtaGapsHisto)fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
2429  if(fDecayTrigger) fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2430 
2431 
2432  if(fHMPIDCorrelation)
2433  {
2434  fhDeltaPhiAssocPtBinHMPID = new TH2F*[fNAssocPtBins*nz];
2435  fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins*nz];
2436  }
2437 
2438  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2439  {
2440  for(Int_t z = 0 ; z < nz ; z++)
2441  {
2442  Int_t bin = i*nz+z;
2443 
2444  if(fCorrelVzBin)
2445  {
2446  sz = Form("_vz%d",z);
2447  tz = Form(", #it{v}_{#it{z}} bin %d",z);
2448  }
2449 
2450  //printf("iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2451 
2452  fhDeltaPhiDeltaEtaAssocPtBin[bin] = new TH2F(Form("hDeltaPhiDeltaEtaPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2453  Form("#Delta #phi vs #Delta #eta for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2454  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
2455  fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2456  fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetYTitle("#Delta #eta");
2457 
2458  fhDeltaPhiAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2459  Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2460  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2461  fhDeltaPhiAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2462  fhDeltaPhiAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2463 
2464  outputContainer->Add(fhDeltaPhiDeltaEtaAssocPtBin[bin]) ;
2465  outputContainer->Add(fhDeltaPhiAssocPtBin[bin]) ;
2466 
2467  if(fFillEtaGapsHisto)
2468  {
2469  fhDeltaPhiAssocPtBinDEta08[bin] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2470  Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2471  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2472  fhDeltaPhiAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2473  fhDeltaPhiAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2474 
2475  fhDeltaPhiAssocPtBinDEta0[bin] = new TH2F(Form("hDeltaPhiDeltaEta0PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2476  Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0.", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2477  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2478  fhDeltaPhiAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2479  fhDeltaPhiAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2480 
2481  outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[bin]) ;
2482  outputContainer->Add(fhDeltaPhiAssocPtBinDEta0[bin]) ;
2483  }
2484 
2485  if(fDecayTrigger)
2486  {
2487  fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s_bit%d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data(),fDecayBits[0]),
2488  Form("#Delta #phi vs #it{p}_{T trigger} tagged as decay for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, Bit %d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data(),fDecayBits[0]),
2489  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2490  fhDeltaPhiDecayChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2491  fhDeltaPhiDecayChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2492 
2493  outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[bin]) ;
2494  }
2495 
2496  if(fFillBradHisto)
2497  {
2498  fhDeltaPhiBradAssocPtBin[bin] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2499  Form("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2500  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2501  fhDeltaPhiBradAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2502  fhDeltaPhiBradAssocPtBin[bin]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2503  outputContainer->Add(fhDeltaPhiBradAssocPtBin[bin]) ;
2504  }
2505 
2506  if(fHMPIDCorrelation)
2507  {
2508  fhDeltaPhiAssocPtBinHMPID[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2509  Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track having HMPID signal", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2510  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2511  fhDeltaPhiAssocPtBinHMPID[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})" );
2512  fhDeltaPhiAssocPtBinHMPID[bin]->SetYTitle("#Delta #phi (rad)");
2513 
2514  fhDeltaPhiAssocPtBinHMPIDAcc[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2515  Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track within 5<phi<20 deg", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2516  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2517  fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2518  fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetYTitle("#Delta #phi (rad)");
2519 
2520  outputContainer->Add(fhDeltaPhiAssocPtBinHMPID [bin]) ;
2521  outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[bin]) ;
2522  }
2523  }
2524  }
2525 
2526  if(fFillMomImbalancePtAssocBinsHisto)
2527  {
2528  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2529  {
2530  fhXEAssocPtBin[i] = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2531  Form("#it{x}_{#it{E}} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2532  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2533  fhXEAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2534  fhXEAssocPtBin[i]->SetYTitle("#it{x}_{#it{E}}");
2535 
2536  fhZTAssocPtBin[i] = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2537  Form("#it{z}_{T} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2538  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2539  fhZTAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2540  fhZTAssocPtBin[i]->SetYTitle("#it{z}_{T}");
2541 
2542 
2543  outputContainer->Add(fhXEAssocPtBin[i]);
2544  outputContainer->Add(fhZTAssocPtBin[i]);
2545  }
2546  }
2547 
2548  if(fCorrelVzBin)
2549  {
2550  for(Int_t z = 0 ; z < nz ; z++)
2551  {
2552  sz = Form("_vz%d",z);
2553  tz = Form(", #it{v}_{#it{z}} bin %d",z);
2554 
2555  fhXEVZ[z] = new TH2F(Form("hXE%s", sz.Data()),
2556  Form("#it{x}_{#it{E}} vs #it{p}_{T trigger}%s", tz.Data()),
2557  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2558  fhXEVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2559  fhXEVZ[z]->SetYTitle("#it{x}_{#it{E}}");
2560 
2561  fhZTVZ[z] = new TH2F(Form("hZT%s",sz.Data()),
2562  Form("#it{z}_{T} vs #it{p}_{T trigger}%s", tz.Data()),
2563  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2564  fhZTVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2565  fhZTVZ[z]->SetYTitle("#it{z}_{T}");
2566 
2567  outputContainer->Add(fhXEVZ[z]);
2568  outputContainer->Add(fhZTVZ[z]);
2569  }
2570  }
2571 
2572  if(fPi0Trigger)
2573  {
2574  fhPtPi0DecayRatio = new TH2F
2575  ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
2576  nptbins,ptmin,ptmax, 100,0.,2.);
2577  fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
2578  fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
2579  outputContainer->Add(fhPtPi0DecayRatio) ;
2580 
2581  fhDeltaPhiPi0DecayCharged = new TH2F
2582  ("hDeltaPhiPi0DecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
2583  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2584  fhDeltaPhiPi0DecayCharged->SetYTitle("#Delta #phi (rad)");
2585  fhDeltaPhiPi0DecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2586 
2587  fhXEPi0DecayCharged =
2588  new TH2F("hXEPi0DecayCharged","#it{x}_{#it{E}} Decay",
2589  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2590  fhXEPi0DecayCharged->SetYTitle("#it{x}_{#it{E}}");
2591  fhXEPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2592 
2593  fhZTPi0DecayCharged =
2594  new TH2F("hZTPi0DecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
2595  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2596  fhZTPi0DecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
2597  fhZTPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2598 
2599  outputContainer->Add(fhDeltaPhiPi0DecayCharged) ;
2600  outputContainer->Add(fhXEPi0DecayCharged) ;
2601  outputContainer->Add(fhZTPi0DecayCharged) ;
2602  }
2603 
2604  if(fDecayTrigger)
2605  {
2606  for(Int_t ibit = 0; ibit< fNDecayBits; ibit++)
2607  {
2608  fhDeltaPhiDecayCharged[ibit] = new TH2F
2609  (Form("hDeltaPhiDecayCharged_bit%d",fDecayBits[ibit]),
2610  Form("#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2611  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2612  fhDeltaPhiDecayCharged[ibit]->SetYTitle("#Delta #phi (rad)");
2613  fhDeltaPhiDecayCharged[ibit]->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2614 
2615  fhXEDecayCharged[ibit] =
2616  new TH2F(Form("hXEDecayCharged_bit%d",fDecayBits[ibit]),
2617  Form("#it{x}_{#it{E}} Decay, Bit %d",fDecayBits[ibit]),
2618  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2619  fhXEDecayCharged[ibit]->SetYTitle("#it{x}_{#it{E}}");
2620  fhXEDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2621 
2622  fhZTDecayCharged[ibit] =
2623  new TH2F(Form("hZTDecayCharged_bit%d",fDecayBits[ibit]),
2624  Form("#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2625  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2626  fhZTDecayCharged[ibit]->SetYTitle("#it{z}_{decay h^{#pm}}");
2627  fhZTDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2628 
2629  outputContainer->Add(fhDeltaPhiDecayCharged[ibit]) ;
2630  outputContainer->Add(fhXEDecayCharged[ibit]) ;
2631  outputContainer->Add(fhZTDecayCharged[ibit]) ;
2632  }
2633  }
2634 
2635  //Correlation with neutral hadrons
2636  if(fNeutralCorr)
2637  {
2638  fhDeltaPhiDeltaEtaNeutral = new TH2F
2639  ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
2640  ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
2641  fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi (rad)");
2642  fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
2643 
2644  fhPhiNeutral = new TH2F
2645  ("hPhiNeutral","#phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}",
2646  nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2647  fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
2648  fhPhiNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2649 
2650  fhEtaNeutral = new TH2F
2651  ("hEtaNeutral","#eta_{#pi^{0}} vs #it{p}_{T #pi^{0}}",
2652  nptbins,ptmin,ptmax,200,-1.,1.);
2653  fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
2654  fhEtaNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2655 
2656  fhDeltaPhiNeutral = new TH2F
2657  ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T trigger}",
2658  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2659  fhDeltaPhiNeutral->SetYTitle("#Delta #phi (rad)");
2660  fhDeltaPhiNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2661 
2662  fhDeltaPhiNeutralPt = new TH2F
2663  ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2664  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2665  fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi (rad)");
2666  fhDeltaPhiNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2667 
2668  fhDeltaEtaNeutral = new TH2F
2669  ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs #it{p}_{T trigger}",
2670  nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2671  fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
2672  fhDeltaEtaNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2673 
2674  fhXENeutral =
2675  new TH2F("hXENeutral","#it{x}_{#it{E}} for #pi^{0} associated",
2676  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2677  fhXENeutral->SetYTitle("#it{x}_{#it{E}}");
2678  fhXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2679 
2680  fhPtHbpXENeutral =
2681  new TH2F("hHbpXENeutral","#xi = ln(1/#it{x}_{#it{E}})for #pi^{0} associated",
2682  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2683  fhPtHbpXENeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2684  fhPtHbpXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2685 
2686  fhZTNeutral =
2687  new TH2F("hZTNeutral","#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2688  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2689  fhZTNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2690  fhZTNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2691 
2692  fhPtHbpZTNeutral =
2693  new TH2F("hHbpZTNeutral","#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2694  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2695  fhPtHbpZTNeutral->SetYTitle("ln(1/#it{z}_{T})");
2696  fhPtHbpZTNeutral->SetXTitle("#it{p}_{T trigger}");
2697 
2698  outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
2699  outputContainer->Add(fhPhiNeutral) ;
2700  outputContainer->Add(fhEtaNeutral) ;
2701  outputContainer->Add(fhDeltaPhiNeutral) ;
2702  outputContainer->Add(fhDeltaPhiNeutralPt) ;
2703  outputContainer->Add(fhDeltaEtaNeutral) ;
2704  outputContainer->Add(fhXENeutral) ;
2705  outputContainer->Add(fhPtHbpXENeutral) ;
2706  outputContainer->Add(fhZTNeutral) ;
2707  outputContainer->Add(fhPtHbpZTNeutral) ;
2708 
2709  fhDeltaPhiUeNeutralPt = new TH2F
2710  (Form("hDeltaPhiUe%sNeutralPt",right.Data()),"#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2711  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2712  fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi (rad)");
2713  fhDeltaPhiUeNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2714 
2715  fhXEUeNeutral =
2716  new TH2F(Form("hXEUeNeutral%s",right.Data()),"#it{x}_{#it{E}} for #pi^{0} associated",
2717  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2718  fhXEUeNeutral->SetYTitle("#it{x}_{#it{E}}");
2719  fhXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2720 
2721  fhPtHbpXEUeNeutral =
2722  new TH2F(Form("hHbpXEUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2723  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2724  fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2725  fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2726 
2727  fhZTUeNeutral =
2728  new TH2F(Form("hZTUeNeutral%s",right.Data()),"#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2729  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2730  fhZTUeNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2731  fhZTUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2732 
2733  fhPtHbpZTUeNeutral =
2734  new TH2F(Form("hHbpZTUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2735  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2736  fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{z}_{T})");
2737  fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2738 
2739  outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
2740  outputContainer->Add(fhXEUeNeutral) ;
2741  outputContainer->Add(fhPtHbpXEUeNeutral) ;
2742  outputContainer->Add(fhZTUeNeutral) ;
2743  outputContainer->Add(fhPtHbpZTUeNeutral) ;
2744 
2745  if(fMakeSeveralUE)
2746  {
2747  fhDeltaPhiUeLeftNeutral = new TH2F
2748  ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs #it{p}_{T h^{0}} with neutral UE left side range of trigger particles",
2749  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2750  fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi (rad)");
2751  fhDeltaPhiUeLeftNeutral->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2752  outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
2753 
2754  fhXEUeLeftNeutral =
2755  new TH2F("hXEUeNeutralLeft","#it{x}_{#it{E}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2756  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2757  fhXEUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2758  fhXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2759  outputContainer->Add(fhXEUeLeftNeutral) ;
2760 
2761  fhPtHbpXEUeLeftNeutral =
2762  new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/#it{x}_{#it{E}}) with neutral UE left side of trigger",
2763  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2764  fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2765  fhPtHbpXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2766  outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
2767 
2768  fhZTUeLeftNeutral =
2769  new TH2F("hZTUeNeutralLeft","#it{z}_{trigger h^{0}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2770  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2771  fhZTUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2772  fhZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2773  outputContainer->Add(fhZTUeLeftNeutral) ;
2774 
2775  fhPtHbpZTUeLeftNeutral =
2776  new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/#it{z}_{T}) with neutral UE left side of trigger",
2777  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2778  fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/#it{z}_{T})");
2779  fhPtHbpZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger}");
2780  outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
2781  }
2782 
2783  if(fPi0Trigger)
2784  {
2785  fhDeltaPhiPi0DecayNeutral = new TH2F
2786  ("hDeltaPhiPi0DecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
2787  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2788  fhDeltaPhiPi0DecayNeutral->SetYTitle("#Delta #phi (rad)");
2789  fhDeltaPhiPi0DecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2790 
2791  fhXEPi0DecayNeutral =
2792  new TH2F("hXEPi0DecayNeutral","#it{x}_{#it{E}} for decay trigger",
2793  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2794  fhXEPi0DecayNeutral->SetYTitle("#it{x}_{#it{E}}");
2795  fhXEPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2796 
2797  fhZTPi0DecayNeutral =
2798  new TH2F("hZTPi0DecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
2799  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2800  fhZTPi0DecayNeutral->SetYTitle("#it{z}_{h^{0}}");
2801  fhZTPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2802 
2803  outputContainer->Add(fhDeltaPhiPi0DecayNeutral) ;
2804  outputContainer->Add(fhXEPi0DecayNeutral) ;
2805  outputContainer->Add(fhZTPi0DecayNeutral) ;
2806  }
2807  }//Correlation with neutral hadrons
2808 
2809  // If data is MC, fill more histograms, depending on origin
2810  if(IsDataMC())
2811  {
2812  for(Int_t i= fMCGenTypeMin; i <= fMCGenTypeMax; i++)
2813  {
2814  fhMCPtTrigger[i] = new TH1F (Form("hMCPtTrigger_%s",nameMC[i].Data()),
2815  Form("MC %s: generated trigger #it{p}_{T}",nameMC[i].Data()),
2816  nptbins,ptmin,ptmax);
2817  fhMCPtTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2818 
2819  fhMCPhiTrigger[i] = new TH2F (Form("hMCPhiTrigger_%s",nameMC[i].Data()),
2820  Form("MC %s: generated trigger #phi",nameMC[i].Data()),
2821  nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2822  fhMCPhiTrigger[i]->SetYTitle("#phi (rad)");
2823  fhMCPhiTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2824 
2825  fhMCEtaTrigger[i] = new TH2F (Form("hMCEtaTrigger_%s",nameMC[i].Data()),
2826  Form("MC %s: generated trigger #eta",nameMC[i].Data()),
2827  nptbins,ptmin,ptmax, netabins,etamin,etamax);
2828  fhMCEtaTrigger[i]->SetYTitle("#eta");
2829  fhMCEtaTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2830 
2831  if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2832  {
2833  fhMCPtTriggerNotLeading[i] = new TH1F (Form("hMCPtTriggerNotLeading_%s",nameMC[i].Data()),
2834  Form("MC %s: generated trigger #it{p}_{T}, when not leading of primaries",nameMC[i].Data()),
2835  nptbins,ptmin,ptmax);
2836  fhMCPtTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2837 
2838  fhMCPhiTriggerNotLeading[i] = new TH2F (Form("hMCPhiTriggerNotLeading_%s",nameMC[i].Data()),
2839  Form("MC %s: generated trigger #phi, when not leading of primaries",nameMC[i].Data()),
2840  nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2841  fhMCPhiTriggerNotLeading[i]->SetYTitle("#phi (rad)");
2842  fhMCPhiTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2843 
2844 
2845  fhMCEtaTriggerNotLeading[i] = new TH2F (Form("hMCEtaTriggerNotLeading_%s",nameMC[i].Data()),
2846  Form("MC %s: generated triogger #eta, when not leading of primaries",nameMC[i].Data()),
2847  nptbins,ptmin,ptmax, netabins,etamin,etamax);
2848  fhMCEtaTriggerNotLeading[i]->SetYTitle("#eta ");
2849  fhMCEtaTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2850  }
2851 
2852  fhMCEtaCharged[i] = new TH2F (Form("hMCEtaCharged_%s",nameMC[i].Data()),
2853  Form("MC %s: #eta_{h^{#pm}} vs #it{p}_{T #pm}",nameMC[i].Data()),
2854  nptbins,ptmin,ptmax,100,-1.,1.);
2855  fhMCEtaCharged[i]->SetYTitle("#eta_{h^{#pm}} (rad)");
2856  fhMCEtaCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2857 
2858  fhMCPhiCharged[i] = new TH2F(Form("hMCPhiCharged_%s",nameMC[i].Data()),
2859  Form("MC %s: phi_{h^{#pm}} vs #it{p}_{T #pm}",nameMC[i].Data()),
2860  nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2861  fhMCPhiCharged[i]->SetYTitle("MC #phi_{h^{#pm}} (rad)");
2862  fhMCPhiCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2863 
2864  fhMCDeltaPhiDeltaEtaCharged[i] = new TH2F (Form("hMCDeltaPhiDeltaEtaCharged_%s",nameMC[i].Data()),
2865  Form("MC %s: phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",nameMC[i].Data()),
2866  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2867  fhMCDeltaPhiDeltaEtaCharged[i]->SetXTitle("#Delta #phi (rad)");
2868  fhMCDeltaPhiDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2869 
2870  fhMCDeltaEtaCharged[i] = new TH2F (Form("hMCDeltaEtaCharged_%s",nameMC[i].Data()),
2871  Form("MC %s: #eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger} and #it{p}_{T assoc}",nameMC[i].Data()),
2872  nptbins,ptmin,ptmax,ndeltaetabins ,deltaetamin,deltaetamax);
2873  fhMCDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2874  fhMCDeltaEtaCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2875 
2876  fhMCDeltaPhiCharged[i] = new TH2F (Form("hMCDeltaPhiCharged_%s",nameMC[i].Data()),
2877  Form("MC %s: #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",nameMC[i].Data()),
2878  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2879  fhMCDeltaPhiCharged[i]->SetYTitle("#Delta #phi (rad)");
2880  fhMCDeltaPhiCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2881 
2882  fhMCDeltaPhiChargedPt[i] = new TH2F (Form("hMCDeltaPhiChargedPt_%s",nameMC[i].Data()),
2883  Form("MC %s: #phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",nameMC[i].Data()),
2884  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2885  fhMCDeltaPhiChargedPt[i]->SetYTitle("#Delta #phi (rad)");
2886  fhMCDeltaPhiChargedPt[i]->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2887 
2888  fhMCPtXECharged[i] = new TH2F (Form("hMCPtXECharged_%s",nameMC[i].Data()),
2889  Form("MC %s: #it{x}_{#it{E}} with charged hadrons",nameMC[i].Data()),
2890  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2891  fhMCPtXECharged[i]->SetYTitle("#it{x}_{#it{E}}");
2892  fhMCPtXECharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2893 
2894  fhMCPtHbpXECharged[i] = new TH2F(Form("hMCHbpXECharged_%s",nameMC[i].Data()),
2895  Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",nameMC[i].Data()),
2896  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2897  fhMCPtHbpXECharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2898  fhMCPtHbpXECharged[i]->SetXTitle("#it{p}_{T trigger}");
2899 
2900  fhMCPtZTCharged[i] = new TH2F(Form("hMCPtZTCharged_%s",nameMC[i].Data()),
2901  Form("MC %s: #it{z}_{T} with charged hadrons",nameMC[i].Data()),
2902  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2903  fhMCPtZTCharged[i]->SetYTitle("#it{z}_{T}");
2904  fhMCPtZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2905 
2906  fhMCPtHbpZTCharged[i] = new TH2F(Form("hMCHbpZTCharged_%s",nameMC[i].Data()),
2907  Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons",nameMC[i].Data()),
2908  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2909  fhMCPtHbpZTCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2910  fhMCPtHbpZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2911 
2912  fhMCPtTrigPout[i] = new TH2F(Form("hMCPtTrigPout_%s",nameMC[i].Data()),
2913  Form("MC %s: #it{p}_{out} with triggers",nameMC[i].Data()),
2914  nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
2915  fhMCPtTrigPout[i]->SetYTitle("#it{p}_{out} (GeV/#it{c})");
2916  fhMCPtTrigPout[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2917 
2918  fhMCPtAssocDeltaPhi[i] = new TH2F(Form("hMCPtAssocDeltaPhi_%s",nameMC[i].Data()),
2919  Form("MC %s: #Delta #phi with associated charged hadrons",nameMC[i].Data()),
2920  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2921  fhMCPtAssocDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
2922  fhMCPtAssocDeltaPhi[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2923 
2924  outputContainer->Add(fhMCPtTrigger[i]);
2925  outputContainer->Add(fhMCPhiTrigger[i]);
2926  outputContainer->Add(fhMCEtaTrigger[i]);
2927 
2928  if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2929  {
2930  outputContainer->Add(fhMCPtTriggerNotLeading[i]);
2931  outputContainer->Add(fhMCPhiTriggerNotLeading[i]);
2932  outputContainer->Add(fhMCEtaTriggerNotLeading[i]);
2933  }
2934 
2935  outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged[i]);
2936  outputContainer->Add(fhMCPhiCharged[i]) ;
2937  outputContainer->Add(fhMCEtaCharged[i]) ;
2938  outputContainer->Add(fhMCDeltaEtaCharged[i]) ;
2939  outputContainer->Add(fhMCDeltaPhiCharged[i]) ;
2940 
2941  outputContainer->Add(fhMCDeltaPhiChargedPt[i]) ;
2942  outputContainer->Add(fhMCPtXECharged[i]) ;
2943  outputContainer->Add(fhMCPtZTCharged[i]) ;
2944  outputContainer->Add(fhMCPtHbpXECharged[i]) ;
2945  outputContainer->Add(fhMCPtHbpZTCharged[i]) ;
2946  outputContainer->Add(fhMCPtTrigPout[i]) ;
2947  outputContainer->Add(fhMCPtAssocDeltaPhi[i]) ;
2948 
2949  // Underlying event
2950 
2951  fhMCUePart[i] =
2952  new TH1F(Form("hMCUePart_%s",nameMC[i].Data()),
2953  Form("MC %s: UE particles distribution vs #it{p}_{T trigger}",nameMC[i].Data()),
2954  nptbins,ptmin,ptmax);
2955  fhMCUePart[i]->SetYTitle("#it{dN}^{ch}");
2956  fhMCUePart[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2957 
2958  fhMCPtXEUeCharged[i] =
2959  new TH2F(Form("hMCPtXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2960  Form("MC %s: #it{x}_{#it{E}} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2961  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2962  fhMCPtXEUeCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2963  fhMCPtXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2964 
2965  fhMCPtHbpXEUeCharged[i] =
2966  new TH2F(Form("hMCPtHbpXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2967  Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2968  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2969  fhMCPtHbpXEUeCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2970  fhMCPtHbpXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2971 
2972  fhMCPtZTUeCharged[i] =
2973  new TH2F(Form("hMCPtZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2974  Form("MC %s: #it{z}_{T} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2975  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2976  fhMCPtZTUeCharged[i]->SetYTitle("#it{z}_{T}");
2977  fhMCPtZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2978 
2979  fhMCPtHbpZTUeCharged[i] =
2980  new TH2F(Form("hMCPtHbpZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2981  Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2982  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2983  fhMCPtHbpZTUeCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2984  fhMCPtHbpZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2985 
2986  outputContainer->Add(fhMCUePart[i]);
2987  outputContainer->Add(fhMCPtXEUeCharged[i]) ;
2988  outputContainer->Add(fhMCPtZTUeCharged[i]) ;
2989  outputContainer->Add(fhMCPtHbpZTUeCharged[i]);
2990  outputContainer->Add(fhMCPtHbpXEUeCharged[i]);
2991 
2992  if(fMakeSeveralUE)
2993  {
2994  fhMCPtXEUeLeftCharged[i] = new TH2F(Form("hMCPtXEUeChargedLeft_%s",nameMC[i].Data()),
2995  Form("MC %s: #it{x}_{#it{E}} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2996  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2997  fhMCPtXEUeLeftCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2998  fhMCPtXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2999 
3000  fhMCPtHbpXEUeLeftCharged[i] = new TH2F(Form("hMCPtHbpXEUeChargedLeft_%s",nameMC[i].Data()),
3001  Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
3002  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
3003  fhMCPtHbpXEUeLeftCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
3004  fhMCPtHbpXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3005 
3006  fhMCPtZTUeLeftCharged[i] = new TH2F(Form("hMCPtZTUeChargedLeft_%s",nameMC[i].Data()),
3007  Form("MC %s: #it{z}_{T} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
3008  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
3009  fhMCPtZTUeLeftCharged[i]->SetYTitle("#it{z}_{T}");
3010  fhMCPtZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3011 
3012  fhMCPtHbpZTUeLeftCharged[i] = new TH2F(Form("hMCPtHbpZTUeChargedLeft_%s",nameMC[i].Data()),
3013  Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
3014  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
3015  fhMCPtHbpZTUeLeftCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
3016  fhMCPtHbpZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3017 
3018  outputContainer->Add(fhMCPtXEUeLeftCharged[i]) ;
3019  outputContainer->Add(fhMCPtZTUeLeftCharged[i]) ;
3020  outputContainer->Add(fhMCPtHbpXEUeLeftCharged[i]);
3021  outputContainer->Add(fhMCPtHbpZTUeLeftCharged[i]) ;
3022 
3023  }
3024  }
3025  } // For MC histogram
3026 
3027  if(DoOwnMix())
3028  {
3029  // Create event containers
3030 
3031  if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
3032  {
3033  Int_t nvz = GetNZvertBin();
3034  Int_t nrp = GetNRPBin();
3035  Int_t nce = GetNCentrBin();
3036 
3037  fListMixTrackEvents= new TList*[nvz*nrp*nce] ;
3038 
3039  for( Int_t ice = 0 ; ice < nce ; ice++ )
3040  {
3041  for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
3042  {
3043  for( Int_t irp = 0 ; irp < nrp ; irp++ )
3044  {
3045  Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
3046 
3047  //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
3048  // ic,iz, irp, bin);
3049 
3050  fListMixTrackEvents[bin] = new TList() ;
3051  fListMixTrackEvents[bin]->SetOwner(kFALSE);
3052  }
3053  }
3054  }
3055  }
3056 
3057  fhPtTriggerMixed = new TH1F ("hPtTriggerMixed","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax);
3058  fhPtTriggerMixed->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
3059 
3060  if(fCorrelVzBin)
3061  {
3062  fhPtTriggerMixedVzBin = new TH2F ("hPtTriggerMixedVzBin","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
3063  fhPtTriggerMixedVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
3064  fhPtTriggerMixedVzBin->SetYTitle("#it{v}_{#it{z}} bin");
3065  outputContainer->Add(fhPtTriggerMixedVzBin);
3066  }
3067 
3068  fhPtTriggerMixedBin = new TH2F ("hPtTriggerMixedBin","#it{p}_{T} distribution of trigger particles vs mixing bin", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
3069  fhPtTriggerMixedBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
3070  fhPtTriggerMixedBin->SetYTitle("Bin");
3071 
3072  fhPhiTriggerMixed = new TH2F ("hPhiTriggerMixed","#phi distribution of trigger Particles, used for mixing",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
3073  fhPhiTriggerMixed->SetYTitle("#phi (rad)");
3074 
3075  fhEtaTriggerMixed = new TH2F ("hEtaTriggerMixed","#eta distribution of trigger, used for mixing",nptbins,ptmin,ptmax, netabins,etamin,etamax);
3076  fhEtaTriggerMixed->SetYTitle("#eta ");
3077 
3078  outputContainer->Add(fhPtTriggerMixed);
3079  outputContainer->Add(fhPtTriggerMixedBin);
3080  outputContainer->Add(fhPhiTriggerMixed);
3081  outputContainer->Add(fhEtaTriggerMixed);
3082 
3083  // Fill the cluster pool only in isolation analysis or if requested
3084  if( neutralMix && (!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForCaloExists())))
3085  {
3086  Int_t nvz = GetNZvertBin();
3087  Int_t nrp = GetNRPBin();
3088  Int_t nce = GetNCentrBin();
3089 
3090  fListMixCaloEvents= new TList*[nvz*nrp*nce] ;
3091 
3092  for( Int_t ice = 0 ; ice < nce ; ice++ )
3093  {
3094  for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
3095  {
3096  for( Int_t irp = 0 ; irp < nrp ; irp++ )
3097  {
3098  Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
3099 
3100  //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
3101  // ic,iz, irp, bin);
3102 
3103  fListMixCaloEvents[bin] = new TList() ;
3104  fListMixCaloEvents[bin]->SetOwner(kFALSE);
3105  }
3106  }
3107  }
3108  }
3109 
3110  // Init the list in the reader if not done previously
3111  if(fUseMixStoredInReader)
3112  {
3113  if( !GetReader()->ListWithMixedEventsForTracksExists() )
3114  GetReader()->SetListWithMixedEventsForTracks(fListMixTrackEvents);
3115 
3116  if( !GetReader()->ListWithMixedEventsForCaloExists() )
3117  GetReader()->SetListWithMixedEventsForCalo (fListMixCaloEvents );
3118  }
3119 
3120  fhEventBin=new TH1I("hEventBin","Number of triggers per bin(cen,vz,rp)",
3121  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3122  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3123  fhEventBin->SetXTitle("event bin");
3124  outputContainer->Add(fhEventBin) ;
3125 
3126  fhEventMixBin=new TH1I("hEventMixBin","Number of triggers mixed per event bin(cen,vz,rp)",
3127  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3128  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3129  fhEventMixBin->SetXTitle("event bin");
3130  outputContainer->Add(fhEventMixBin) ;
3131 
3132  fhEventMBBin=new TH1I("hEventMBBin","Number of min bias events per bin(cen,vz,rp)",
3133  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3134  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3135  fhEventMBBin->SetXTitle("event bin");
3136  outputContainer->Add(fhEventMBBin) ;
3137 
3138  fhNtracksMB=new TH2F("hNtracksMBEvent","Number of filtered tracks in MB event per event bin",ntrbins,trmin,trmax,
3139  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3140  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3141  fhNtracksMB->SetYTitle("event bin");
3142  fhNtracksMB->SetXTitle("#it{N}_{track}");
3143  outputContainer->Add(fhNtracksMB);
3144 
3145  if( neutralMix )
3146  {
3147  fhNclustersMB=new TH2F("hNclustersMBEvent","Number of filtered clusters in MB events per event bin",nclbins,clmin,clmax,
3148  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3149  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3150  fhNclustersMB->SetYTitle("event bin");
3151  fhNclustersMB->SetXTitle("#it{N}_{cluster}");
3152  outputContainer->Add(fhNclustersMB);
3153  }
3154 
3155  fhMixDeltaPhiCharged = new TH2F
3156  ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
3157  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
3158  fhMixDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
3159  fhMixDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3160  outputContainer->Add(fhMixDeltaPhiCharged);
3161 
3162  fhMixDeltaPhiDeltaEtaCharged = new TH2F
3163  ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
3164  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
3165  fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
3166  fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
3167  outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
3168 
3169  fhMixXECharged =
3170  new TH2F("hMixXECharged","Mixed event : #it{x}_{#it{E}} for charged tracks",
3171  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
3172  fhMixXECharged->SetYTitle("#it{x}_{#it{E}}");
3173  fhMixXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3174  outputContainer->Add(fhMixXECharged);
3175 
3176  fhMixXEUeCharged =
3177  new TH2F("hMixXEUeCharged","Mixed event : #it{x}_{#it{E}} for charged tracks in Ue region",
3178  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
3179  fhMixXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
3180  fhMixXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3181  outputContainer->Add(fhMixXEUeCharged);
3182 
3183  fhMixHbpXECharged =
3184  new TH2F("hMixHbpXECharged","mixed event : #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
3185  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
3186  fhMixHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
3187  fhMixHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3188  outputContainer->Add(fhMixHbpXECharged);
3189 
3190  fhMixDeltaPhiChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
3191  fhMixDeltaPhiChargedAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
3192  fhMixDeltaPhiChargedAssocPtBinDEta0 = new TH2F*[fNAssocPtBins*nz];
3193  fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
3194 
3195  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3196  {
3197  for(Int_t z = 0 ; z < nz ; z++)
3198  {
3199  Int_t bin = i*nz+z;
3200 
3201  if(fCorrelVzBin)
3202  {
3203  sz = Form("_vz%d",z);
3204  tz = Form(", #it{v}_{#it{z}} bin %d",z);
3205  }
3206 
3207  //printf("MIX : iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
3208 
3209  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3210  Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
3211  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
3212  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
3213  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetYTitle("#Delta #eta");
3214 
3215  outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]);
3216 
3217  fhMixDeltaPhiChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3218  Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
3219  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
3220  fhMixDeltaPhiChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3221  fhMixDeltaPhiChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
3222 
3223  outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[bin]);
3224 
3225  if(fFillEtaGapsHisto)
3226  {
3227  fhMixDeltaPhiChargedAssocPtBinDEta08[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3228  Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
3229  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
3230  fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3231  fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
3232 
3233  fhMixDeltaPhiChargedAssocPtBinDEta0[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3234  Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
3235  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
3236  fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3237  fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
3238 
3239  outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[bin]);
3240  outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta0[bin]);
3241  }
3242  }
3243  }
3244  }
3245 
3246  return outputContainer;
3247 }
3248 
3249 //_____________________________________________________________________________________________________________________
3257 //_____________________________________________________________________________________________________________________
3258 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(Int_t indexPhoton1, Int_t indexPhoton2, Int_t idetector)
3259 {
3260  if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
3261 
3262  AliDebug(1,Form("indexPhoton1 = %d, indexPhoton2 = %d", indexPhoton1, indexPhoton2));
3263 
3264  TObjArray * clusters = 0x0 ;
3265  if(idetector==kEMCAL) clusters = GetEMCALClusters() ;
3266  else clusters = GetPHOSClusters() ;
3267 
3268  for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
3269  {
3270  AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
3271 
3272  if(photon->GetID()==indexPhoton1) photon->GetMomentum(fDecayMom1,GetVertex(0)) ;
3273  if(photon->GetID()==indexPhoton2) photon->GetMomentum(fDecayMom2,GetVertex(0)) ;
3274 
3275  AliDebug(1,Form("Photon1 = %f, Photon2 = %f", fDecayMom1.Pt(), fDecayMom2.Pt()));
3276  } //cluster loop
3277 
3278  return kTRUE;
3279 }
3280 
3281 //________________________________________________________________________
3284 //________________________________________________________________________
3285 Int_t AliAnaParticleHadronCorrelation::GetMCTagHistogramIndex(Int_t mcTag)
3286 {
3287  if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||
3288  GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation) ||
3289  GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCISR) ) return 0;
3290  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ) return 1;
3291  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay) ) return 2;
3292  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ) return 3;
3293  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay) ) return 4;
3294  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) ) return 5; // other decays
3295  else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) ) return 6;
3296  else return 7;
3297 }
3298 
3299 //_________________________________________
3301 //_________________________________________
3302 void AliAnaParticleHadronCorrelation::Init()
3303 {
3304  if(!GetReader()->IsCTSSwitchedOn())
3305  AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!");
3306 }
3307 
3308 //____________________________________________________
3310 //____________________________________________________
3311 void AliAnaParticleHadronCorrelation::InitParameters()
3312 {
3313  SetInputAODName("Particle");
3314  SetAODObjArrayName("Hadrons");
3315  AddToHistogramsName("AnaHadronCorr_");
3316 
3317  SetPtCutRange(0.,300);
3318  fDeltaPhiMinCut = TMath::DegToRad()*120.;
3319  fDeltaPhiMaxCut = TMath::DegToRad()*240. ;
3320  fSelectIsolated = kFALSE;
3321  fMakeSeveralUE = kFALSE;
3322  fUeDeltaPhiMinCut = TMath::DegToRad()*60.;
3323  fUeDeltaPhiMaxCut = TMath::DegToRad()*120 ;
3324 
3325  fNeutralCorr = kFALSE ;
3326  fPi0Trigger = kFALSE ;
3327  fDecayTrigger = kFALSE ;
3328  fHMPIDCorrelation = kFALSE ;
3329 
3330  fMakeAbsoluteLeading = kTRUE;
3331  fMakeNearSideLeading = kFALSE;
3332 
3333  fNAssocPtBins = 9 ;
3334  fAssocPtBinLimit[0] = 0.2 ;
3335  fAssocPtBinLimit[1] = 0.5 ;
3336  fAssocPtBinLimit[2] = 1.0 ;
3337  fAssocPtBinLimit[3] = 2.0 ;
3338  fAssocPtBinLimit[4] = 3.0 ;
3339  fAssocPtBinLimit[5] = 4.0 ;
3340  fAssocPtBinLimit[6] = 5.0 ;
3341  fAssocPtBinLimit[7] = 6.0 ;
3342  fAssocPtBinLimit[8] = 7.0 ;
3343  fAssocPtBinLimit[9] = 8.0 ;
3344  fAssocPtBinLimit[10] = 9.0 ;
3345  fAssocPtBinLimit[11] = 10.0 ;
3346  fAssocPtBinLimit[12] = 12.0 ;
3347  fAssocPtBinLimit[13] = 14.0 ;
3348  fAssocPtBinLimit[14] = 16.0 ;
3349  fAssocPtBinLimit[15] = 20.0 ;
3350  fAssocPtBinLimit[16] = 30.0 ;
3351  fAssocPtBinLimit[17] = 40.0 ;
3352  fAssocPtBinLimit[18] = 50.0 ;
3353  fAssocPtBinLimit[19] = 200.0 ;
3354 
3355  fUseMixStoredInReader = kTRUE;
3356 
3357  fM02MinCut = -1 ;
3358  fM02MaxCut = -1 ;
3359 
3360  fDecayTagsM02Cut = 0.27;
3361 
3362  fSelectLeadingHadronAngle = kFALSE;
3363  fFillLeadHadOppositeHisto = kFALSE;
3364  fMinLeadHadPhi = 150*TMath::DegToRad();
3365  fMaxLeadHadPhi = 210*TMath::DegToRad();
3366 
3367  fMinLeadHadPt = 1;
3368  fMaxLeadHadPt = 100;
3369 
3370  fMCGenTypeMin = 0;
3371  fMCGenTypeMax = fgkNmcTypes-1;
3372 
3373  fNDecayBits = 1;
3374  fDecayBits[0] = AliNeutralMesonSelection::kPi0;
3375  fDecayBits[1] = AliNeutralMesonSelection::kEta;
3376  fDecayBits[2] = AliNeutralMesonSelection::kPi0RightSide;
3377  fDecayBits[3] = AliNeutralMesonSelection::kEtaRightSide;
3378  fDecayBits[4] = AliNeutralMesonSelection::kEtaLeftSide ;
3379  fDecayBits[5] = AliNeutralMesonSelection::kEtaBothSides;
3380  fDecayBits[6] = AliNeutralMesonSelection::kPi0LeftSide ; // Leave it last since likely not used
3381  fDecayBits[7] = AliNeutralMesonSelection::kPi0BothSides; // Leave it last since likely not used
3382 
3383  fNBkgBin = 11;
3384  fBkgBinLimit[ 0] = 00.0; fBkgBinLimit[ 1] = 00.2; fBkgBinLimit[ 2] = 00.3; fBkgBinLimit[ 3] = 00.4; fBkgBinLimit[ 4] = 00.5;
3385  fBkgBinLimit[ 5] = 01.0; fBkgBinLimit[ 6] = 01.5; fBkgBinLimit[ 7] = 02.0; fBkgBinLimit[ 8] = 03.0; fBkgBinLimit[ 9] = 05.0;
3386  fBkgBinLimit[10] = 10.0; fBkgBinLimit[11] = 100.;
3387  for(Int_t ibin = fNBkgBin+1; ibin < 20; ibin++)
3388  {
3389  fBkgBinLimit[ibin] = 00.0;
3390  }
3391 
3392 }
3393 
3394 //________________________________________________________________________________________________________
3399 //________________________________________________________________________________________________________
3400 void AliAnaParticleHadronCorrelation::InvMassHisto(AliAODPWG4ParticleCorrelation * trigger, Int_t mcIndex)
3401 {
3402  Int_t idTrig = trigger->GetCaloLabel(0);
3403  Float_t ptTrig = trigger->Pt();
3404  Float_t tofTrig = trigger->GetTime();
3405 
3406  fMomentum = *(trigger->Momentum());
3407 
3408  // Loop on second cluster to be associated to trigger particle to have the invariant mass
3409  Int_t nphoton = GetInputAODBranch()->GetEntriesFast();
3410 
3411  for(Int_t iphoton = 0; iphoton < nphoton; iphoton++)
3412  {
3413  AliAODPWG4ParticleCorrelation * photon1 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iphoton));
3414 
3415  if(idTrig == photon1->GetCaloLabel(0)) continue; //Do not the invariant mass for the same particle
3416 
3417  fMomentumIM = *(photon1->Momentum());
3418 
3419  // Select secondary photon with proper invariant mass
3420  if(fM02MaxCut > 0 && fM02MinCut > 0) //clID1 > 0 && clID2 < 0 &&
3421  {
3422  Float_t m02 = photon1->GetM02();
3423 
3424  if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3425  }
3426 
3427  // Select clusters with good time window difference
3428  Double_t tdiff = tofTrig - photon1->GetTime();
3429  if(TMath::Abs(tdiff) > GetPairTimeCut()) continue;
3430 
3431  // Add momenta of trigger and possible secondary decay
3432  fMomentumIM += fMomentum;
3433 
3434  Double_t mass = fMomentumIM.M();
3435  // Double_t epair = fMomentumIM.E();
3436  // Float_t ptpair = fMomentumIM.Pt();
3437 
3438  // Mass of all pairs
3439  fhMassPtTrigger->Fill(ptTrig, mass, GetEventWeight());
3440  if ( IsDataMC() ) fhMCMassPtTrigger[mcIndex]->Fill(ptTrig, mass, GetEventWeight());
3441  }
3442 }
3443 
3444 //_________________________________________________________________________
3447 //_________________________________________________________________________
3448 Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
3449 {
3450  Double_t ptTrig = GetMinPt();
3451  Double_t phiTrig = 0 ;
3452  fLeadingTriggerIndex =-1 ;
3453  Int_t index =-1 ;
3454  AliAODPWG4ParticleCorrelation* pLeading = 0;
3455 
3456  // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
3457 
3458  for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3459  {
3460  AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3461  particle->SetLeadingParticle(kFALSE); // set it later
3462 
3463  // Vertex cut in case of mixing
3464  Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3465  if(check == 0) continue;
3466  if(check == -1) return kFALSE; // not sure if it is correct.
3467 
3468  // Find the leading particles with highest momentum
3469  if (particle->Pt() > ptTrig)
3470  {
3471  ptTrig = particle->Pt() ;
3472  phiTrig = particle->Phi();
3473  index = iaod ;
3474  pLeading = particle ;
3475  }
3476  }// Finish search of leading trigger particle on the AOD branch.
3477 
3478  if(index < 0) return kFALSE;
3479 
3480  //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3481 
3482  if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3483 
3484  // Compare if it is the leading of all tracks
3485 
3486  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3487  {
3488  AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3489 
3490  // In case of isolation of single tracks or conversion photon (2 tracks) or pi0 (4 tracks),
3491  // do not count the candidate or the daughters of the candidate
3492  // in the isolation conte
3493  if ( pLeading->GetDetectorTag() == kCTS ) // make sure conversions are tagged as kCTS!!!
3494  {
3495  Int_t trackID = GetReader()->GetTrackID(track) ; // needed instead of track->GetID() since AOD needs some manipulations
3496  Bool_t contained = kFALSE;
3497 
3498  for(Int_t i = 0; i < 4; i++)
3499  {
3500  if( trackID == pLeading->GetTrackLabel(i) ) contained = kTRUE;
3501  }
3502 
3503  if ( contained ) continue ;
3504  }
3505 
3506  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3507  Float_t pt = fTrackVector.Pt();
3508  Float_t phi = fTrackVector.Phi() ;
3509  if(phi < 0) phi+=TMath::TwoPi();
3510 
3511  //jump out this event if near side associated particle pt larger than trigger
3512  if (fMakeNearSideLeading)
3513  {
3514  Float_t deltaPhi = phiTrig-phi;
3515  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3516  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3517 
3518  if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE;
3519  }
3520  //jump out this event if there is any other particle with pt larger than trigger
3521  else
3522  {
3523  if(pt > ptTrig) return kFALSE ;
3524  }
3525  }// track loop
3526 
3527  // Compare if it is leading of all calorimeter clusters
3528 
3529  if(fCheckLeadingWithNeutralClusters)
3530  {
3531  // Select the calorimeter cluster list
3532  TObjArray * nePl = 0x0;
3533  if (pLeading->GetDetectorTag() == kPHOS )
3534  nePl = GetPHOSClusters();
3535  else
3536  nePl = GetEMCALClusters();
3537 
3538  if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3539 
3540  for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3541  {
3542  AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3543 
3544  if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3545 
3546  cluster->GetMomentum(fMomentum,GetVertex(0));
3547 
3548  Float_t pt = fMomentum.Pt();
3549  Float_t phi = fMomentum.Phi() ;
3550  if(phi < 0) phi+=TMath::TwoPi();
3551 
3552  if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3553 
3554  //jump out this event if near side associated particle pt larger than trigger
3555  // not really needed for calorimeter, unless DCal is included
3556  if (fMakeNearSideLeading)
3557  {
3558  Float_t deltaPhi = phiTrig-phi;
3559  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3560  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3561 
3562  if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3563  }
3564  //jump out this event if there is any other particle with pt larger than trigger
3565  else
3566  {
3567  if(pt > ptTrig) return kFALSE ;
3568  }
3569  }// cluster loop
3570  } // check neutral clusters
3571 
3572  fLeadingTriggerIndex = index ;
3573  pLeading->SetLeadingParticle(kTRUE);
3574 
3575  AliDebug(1,Form("\t particle AOD with index %d is leading with pT %2.2f", fLeadingTriggerIndex, pLeading->Pt()));
3576 
3577  return kTRUE;
3578 }
3579 
3580 //_________________________________________________________________
3594 //_________________________________________________________________
3595 void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
3596 {
3597  if(!GetInputAODBranch())
3598  {
3599  AliFatal(Form("No input particles in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3600  return ; // coverity
3601  }
3602 
3603  Int_t naod = GetInputAODBranch()->GetEntriesFast();
3604  if( naod == 0 )
3605  {
3606  AliDebug(1,"No particle AOD found!");
3607  return ; // no trigger particles found.
3608  }
3609 
3610  AliDebug(1,Form("Begin hadron correlation analysis, fill histograms"));
3611  AliDebug(1,Form("n particle branch aod entries %d", naod));
3612  AliDebug(1,Form("In CTS aod entries %d",GetCTSTracks()->GetEntriesFast()));
3613 
3614  //------------------------------------------------------
3615  // Find leading trigger if analysis request only leading,
3616  // if there is no leading trigger, then skip the event
3617 
3618  Int_t iaod = 0 ;
3619  if( fMakeAbsoluteLeading || fMakeNearSideLeading )
3620  {
3621  Bool_t leading = IsTriggerTheEventLeadingParticle();
3622 
3623  AliDebug(1,Form("AOD Leading trigger? %d, with index %d",leading,fLeadingTriggerIndex));
3624 
3625  if(!leading)
3626  {
3627  AliDebug(1,"Leading was requested and not found");
3628  return ;
3629  }
3630  else
3631  {
3632  // Select only the leading in the trigger AOD loop
3633  naod = fLeadingTriggerIndex+1 ;
3634  iaod = fLeadingTriggerIndex ;
3635  }
3636  }
3637 
3638  //------------------------------------------------------
3639  // Get event multiplicity and bins
3640 
3641  Float_t cen = GetEventCentrality();
3642  Float_t ep = GetEventPlaneAngle();
3643  if(IsHighMultiplicityAnalysisOn()) fhTriggerEventPlaneCentrality->Fill(cen, ep, GetEventWeight());
3644 
3645  Int_t mixEventBin = GetEventMixBin();
3646  Int_t vzbin = GetEventVzBin();
3647 
3648  //------------------------------------------------------
3649  // Loop on trigger AOD
3650 
3651  for( iaod = 0; iaod < naod; iaod++ )
3652  {
3653  AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3654 
3655  //
3656  // Trigger particle selection criteria:
3657  //
3658  Float_t pt = particle->Pt();
3659 
3660  if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3661 
3662  fhPtTriggerInput->Fill(pt, GetEventWeight());
3663 
3664  //
3665  // check if it was a calorimeter cluster
3666  // and if the shower shape cut was requested apply it.
3667  // Not needed if already done at the particle identification level,
3668  // but for isolation studies, it is preferred not to remove so we do it here
3669  //
3670 
3671  AliDebug(1,Form("%s Trigger : min %f, max %f, det %d",
3672  GetInputAODName().Data(),fM02MinCut,fM02MaxCut,particle->GetDetectorTag()));
3673 
3674  if(fM02MaxCut > 0 && fM02MinCut > 0) //clID1 > 0 && clID2 < 0 &&
3675  {
3676  // Int_t iclus = -1;
3677  // TObjArray* clusters = 0x0;
3678  // if (particle->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
3679  // else if(particle->GetDetectorTag() == kPHOS ) clusters = GetPHOSClusters();
3680  //
3681  // if(clusters)
3682  // {
3683  // AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
3684  // Float_t m02 = cluster->GetM02();
3685  // if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3686  // }
3687 
3688  Float_t m02 = particle->GetM02();
3689 
3690  if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3691 
3692  fhPtTriggerSSCut->Fill(pt, GetEventWeight());
3693 
3694  AliDebug(1,"Pass the shower shape cut");
3695  }
3696 
3697  //
3698  // Check if the particle is isolated or if we want to take the isolation into account
3699  // This bool is set in AliAnaParticleIsolation
3700  //
3701  if(OnlyIsolated())
3702  {
3703  if( !particle->IsIsolated() ) continue;
3704 
3705  fhPtTriggerIsoCut->Fill(pt, GetEventWeight());
3706 
3707  AliDebug(1,"Pass the isolation cut");
3708  }
3709 
3710  //
3711  // Check if trigger is in fiducial region
3712  //
3713  if(IsFiducialCutOn())
3714  {
3715  Bool_t in = GetFiducialCut()->IsInFiducialCut(particle->Eta(),particle->Phi(),particle->GetDetectorTag()) ;
3716 
3717  if(! in ) continue ;
3718 
3719  AliDebug(1,"Pass the fiducial cut");
3720  }
3721 
3722  fhPtTriggerFidCut->Fill(pt, GetEventWeight());
3723 
3724  //---------------------------------------
3725  // Make correlation
3726 
3727  // Find the leading hadron in the opposite hemisphere to the trigger
3728  // and accept the trigger if leading is in defined window.
3729  Bool_t okLeadHad = kTRUE;
3730  if(fSelectLeadingHadronAngle || fFillLeadHadOppositeHisto)
3731  {
3732  okLeadHad = FindLeadingOppositeHadronInWindow(particle);
3733  if(!okLeadHad && fSelectLeadingHadronAngle) continue;
3734  }
3735 
3736  //
3737  // Charged particles correlation
3738  //
3739  MakeChargedCorrelation(particle);
3740 
3741  // MC
3742  Int_t mcIndex = -1;
3743  Int_t mcTag = particle->GetTag();
3744  Bool_t lostDecayPair = kFALSE;
3745  if(IsDataMC())
3746  {
3747  mcIndex = GetMCTagHistogramIndex(mcTag);
3748  lostDecayPair = GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost);
3749  MakeMCChargedCorrelation(particle->GetLabel(), mcIndex,lostDecayPair);
3750  }
3751 
3752  // Do own mixed event with charged,
3753  // add event and remove previous or fill the mixed histograms
3754  if(DoOwnMix())
3755  MakeChargedMixCorrelation(particle);
3756 
3757  //
3758  // Neutral particles correlation
3759  //
3760  if(fNeutralCorr)
3761  MakeNeutralCorrelation(particle);
3762 
3763  //----------------------------------------------------------------
3764  // Fill trigger pT related histograms if not absolute leading
3765 
3766  //
3767  // pT of the trigger, vs trigger origin if MC
3768  //
3769  fhPtTrigger->Fill(pt, GetEventWeight());
3770 
3771 
3772  if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3773  {
3774  fhPtTriggerMC[mcIndex]->Fill(pt, GetEventWeight());
3775  if( lostDecayPair )
3776  {
3777  // check index in GetMCTagIndexHistogram
3778  if ( mcIndex == 2 ) fhPtTriggerMC[8]->Fill(pt, GetEventWeight()); // pi0 decay
3779  else if( mcIndex == 4 ) fhPtTriggerMC[9]->Fill(pt, GetEventWeight()); // eta decay
3780  }
3781  }
3782 
3783  //
3784  // pT lead cone, pT sum cone background bins distribution
3785  //
3786  if(fFillBkgBinsHisto)
3787  {
3788  Float_t m02 = particle->GetM02();
3789  Float_t pTLeadTrackInCone = 0;
3790  Float_t pTSumTrackInCone = 0;
3791  Float_t pTLeadClusterInCone = 0;
3792  Float_t pTSumClusterInCone = 0;
3793 
3794  pTLeadTrackInCone = particle->GetChargedLeadPtInCone();
3795  pTLeadClusterInCone = particle->GetNeutralLeadPtInCone();
3796 
3797  pTSumTrackInCone = particle->GetChargedPtSumInCone();
3798  pTSumClusterInCone = particle->GetNeutralPtSumInCone();
3799 
3800  Float_t pTLeadInCone = pTLeadTrackInCone;
3801  if(pTLeadClusterInCone > pTLeadInCone) pTLeadInCone = pTLeadClusterInCone;
3802  Float_t pTSumInCone = pTSumTrackInCone + pTSumClusterInCone;
3803 
3804  // Get the background bin for this cone and trigger
3805  Int_t pTSumBin = -1;
3806  Int_t pTLeadBin = -1;
3807 
3808  for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
3809  {
3810  if( pTSumInCone >= fBkgBinLimit[ibin] && pTSumInCone < fBkgBinLimit[ibin+1]) pTSumBin = ibin;
3811  if( pTLeadInCone >= fBkgBinLimit[ibin] && pTLeadInCone < fBkgBinLimit[ibin+1]) pTLeadBin = ibin;
3812  }
3813  if(pTSumBin > 0)
3814  {
3815  fhPtSumInConeBin[pTSumBin]->Fill(pt, GetEventWeight());
3816  }
3817 
3818  if(pTLeadBin > 0)
3819  {
3820  fhPtLeadInConeBin[pTLeadBin]->Fill(pt, GetEventWeight());
3821  }
3822 
3823  // Check if it was a decay
3824  if( fFillTaggedDecayHistograms && m02 < fDecayTagsM02Cut )
3825  {
3826  Int_t decayTag = particle->DecayTag();
3827  if(decayTag < 0) decayTag = 0;
3828 
3829  for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
3830  {
3831  if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3832  {
3833  Int_t pTLeadBinDecay = pTLeadBin+ibit*fNBkgBin;
3834  Int_t pTSumBinDecay = pTSumBin+ibit*fNBkgBin;
3835  if( pTLeadBin >=0 ) fhPtLeadConeBinDecay[pTLeadBinDecay]->Fill(pt, GetEventWeight());
3836  if( pTSumBin >=0 ) fhSumPtConeBinDecay [pTSumBinDecay] ->Fill(pt, GetEventWeight());
3837  }
3838  }
3839  }
3840 
3841  if(IsDataMC())
3842  {
3843  Int_t pTLeadBinMC = pTLeadBin+mcIndex*fNBkgBin;
3844  Int_t pTSumBinMC = pTSumBin+mcIndex*fNBkgBin;
3845 
3846  if( pTLeadBin >=0 )
3847  {
3848  fhPtLeadConeBinMC[pTLeadBinMC]->Fill(pt, GetEventWeight());
3849  }
3850 
3851  if( pTSumBin >=0 )
3852  {
3853  fhSumPtConeBinMC [pTSumBinMC]->Fill(pt, GetEventWeight());
3854  }
3855 
3856  if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3857  {
3858  pTLeadBinMC = pTLeadBin+kmcPhoton*fNBkgBin;
3859  pTSumBinMC = pTSumBin+kmcPhoton*fNBkgBin;
3860  if( pTLeadBin >=0 )
3861  {
3862  fhPtLeadConeBinMC[pTLeadBinMC]->Fill(pt, GetEventWeight());
3863  }
3864 
3865  if( pTSumBin >=0 )
3866  {
3867  fhSumPtConeBinMC [ pTSumBinMC]->Fill(pt, GetEventWeight());
3868  }
3869  }
3870 
3871  // Check if decay and if pair is lost
3872  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
3873  {
3874  if ( mcIndex == kmcPi0Decay )
3875  {
3876  pTLeadBinMC = pTLeadBin+kmcPi0DecayLostPair*fNBkgBin;
3877  pTSumBinMC = pTSumBin+kmcPi0DecayLostPair*fNBkgBin;
3878  }
3879  else if(mcIndex == kmcEtaDecay)
3880  {
3881  pTLeadBinMC = pTLeadBin+kmcEtaDecayLostPair*fNBkgBin;
3882  pTSumBinMC = pTSumBin+kmcEtaDecayLostPair*fNBkgBin;
3883  }
3884  else
3885  AliFatal(Form("Lost decay Bit assigned to bad case, mcIndex %d",mcIndex));
3886 
3887  if( pTLeadBin >=0 )
3888  {
3889  fhPtLeadConeBinMC[pTLeadBinMC]->Fill(pt, GetEventWeight());
3890  }
3891 
3892  if( pTSumBin >=0 )
3893  {
3894  fhSumPtConeBinMC [ pTSumBinMC]->Fill(pt, GetEventWeight());
3895  }
3896 
3897  }
3898 
3899  }
3900  }
3901 
3902  // Invariant mass of trigger particle
3903  if(fFillInvMassHisto) InvMassHisto(particle,mcIndex);
3904 
3905  if(fDecayTrigger)
3906  {
3907  Int_t decayTag = particle->DecayTag();
3908  if(decayTag < 0) decayTag = 0;
3909 
3910  for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
3911  {
3912  if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3913  {
3914  fhPtDecayTrigger[ibit]->Fill(pt, GetEventWeight());
3915 
3916  if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3917  {
3918  fhPtDecayTriggerMC[ibit][mcIndex]->Fill(pt, GetEventWeight());
3919  if( lostDecayPair )
3920  {
3921  // check index in GetMCTagIndexHistogram
3922  if( mcIndex == 2 ) fhPtDecayTriggerMC[ibit][8]->Fill(pt, GetEventWeight()); // pi0 decay
3923  if( mcIndex == 4 ) fhPtDecayTriggerMC[ibit][9]->Fill(pt, GetEventWeight()); // eta decay
3924  }
3925  }
3926  }// check bit
3927  }// bit loop
3928  }
3929 
3930  //
3931  // Acceptance of the trigger
3932  //
3933  Float_t phi = particle->Phi();
3934  if( phi < 0 ) phi+=TMath::TwoPi();
3935  fhPhiTrigger->Fill(pt, phi, GetEventWeight());
3936 
3937  fhEtaTrigger->Fill(pt, particle->Eta(), GetEventWeight());
3938  //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Trigger particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
3939 
3940  //----------------------------------
3941  // Trigger particle pT vs event bins
3942 
3943  fhPtTriggerBin->Fill(pt, mixEventBin, GetEventWeight());
3944  if(fCorrelVzBin)
3945  fhPtTriggerVzBin->Fill(pt, vzbin, GetEventWeight());
3946 
3947  if(IsHighMultiplicityAnalysisOn())
3948  {
3949  fhPtTriggerCentrality->Fill(pt, cen, GetEventWeight());
3950  fhPtTriggerEventPlane->Fill(pt, ep , GetEventWeight());
3951  }
3952 
3953  //----------------------------------
3954  // Trigger particle pT vs pile-up
3955 
3956  if(IsPileUpAnalysisOn())
3957  {
3958  Int_t vtxBC = GetReader()->GetVertexBC();
3959  if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA) fhPtTriggerVtxBC0->Fill(pt, GetEventWeight());
3960 
3961  if(GetReader()->IsPileUpFromSPD()) fhPtTriggerPileUp[0]->Fill(pt, GetEventWeight());
3962  if(GetReader()->IsPileUpFromEMCal()) fhPtTriggerPileUp[1]->Fill(pt, GetEventWeight());
3963  if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtTriggerPileUp[2]->Fill(pt, GetEventWeight());
3964  if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtTriggerPileUp[3]->Fill(pt, GetEventWeight());
3965  if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtTriggerPileUp[4]->Fill(pt, GetEventWeight());
3966  if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtTriggerPileUp[5]->Fill(pt, GetEventWeight());
3967  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTriggerPileUp[6]->Fill(pt, GetEventWeight());
3968  }
3969  } // AOD trigger loop
3970 
3971  // Re-init for next event
3972  fLeadingTriggerIndex = -1;
3973 
3974  AliDebug(1,"End fill histograms");
3975 }
3976 
3977 //_______________________________________________________________________________________________________
3983 //_______________________________________________________________________________________________________
3984 void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3985 {
3986  AliDebug(1,"Make trigger particle - charged hadron correlation");
3987 
3988  Float_t phiTrig = aodParticle->Phi();
3989  Float_t etaTrig = aodParticle->Eta();
3990  Float_t ptTrig = aodParticle->Pt();
3991  Int_t mcTag = aodParticle->GetTag();
3992  Double_t bz = GetReader()->GetInputEvent()->GetMagneticField();
3993 
3994 
3995  Int_t decayTag = 0;
3996  if(fDecayTrigger)
3997  {
3998  //decay = aodParticle->IsTagged();
3999  decayTag = aodParticle->DecayTag();
4000  if(decayTag < 0) decayTag = 0;
4001  // printf("Correlation: pT %2.2f, BTag %d, Tagged %d\n",ptTrig, decayTag, aodParticle->IsTagged());
4002  // printf("\t check bit Pi0 %d, Eta %d, Pi0Side %d, EtaSide %d\n",
4003  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0),
4004  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEta),
4005  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0Side),
4006  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEtaSide));
4007  }
4008 
4009  Float_t pt = -100. ;
4010  Float_t phi = -100. ;
4011  Float_t eta = -100. ;
4012  Float_t deltaPhi = -100. ;
4013 
4014  TObjArray * reftracks = 0x0;
4015  Int_t nrefs = 0;
4016 
4017  // Mixed event settings
4018  Int_t evtIndex11 = -1 ; // cluster trigger or pi0 trigger
4019  Int_t evtIndex12 = -1 ; // pi0 trigger
4020  Int_t evtIndex13 = -1 ; // charged trigger
4021 
4022  if (GetMixedEvent())
4023  {
4024  evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
4025  evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
4026  evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
4027  }
4028 
4029  // Track multiplicity or cent bin
4030  Int_t cenbin = 0;
4031  if(IsHighMultiplicityAnalysisOn()) cenbin = GetEventCentralityBin();
4032 
4033  //
4034  // In case of pi0/eta trigger, we may want to check their decay correlation,
4035  // get their decay children
4036  //
4037 
4038  Bool_t decayFound = kFALSE;
4039  if( fPi0Trigger )
4040  {
4041  decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
4042  if(decayFound)
4043  {
4044  fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom1.Pt()/ptTrig, GetEventWeight());
4045  fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom2.Pt()/ptTrig, GetEventWeight());
4046  }
4047  }
4048 
4049  //-----------------------------------------------------------------------
4050  // Track loop, select tracks with good pt, phi and fill AODs or histograms
4051  //-----------------------------------------------------------------------
4052 
4053  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
4054  {
4055  AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
4056 
4057  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
4058  pt = fTrackVector.Pt();
4059  eta = fTrackVector.Eta();
4060  phi = fTrackVector.Phi() ;
4061  if(phi < 0) phi+=TMath::TwoPi();
4062 
4063  //Select only hadrons in pt range
4064  if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
4065 
4066  // In case of isolation of single tracks or conversion photon (2 tracks) or pi0 (4 tracks),
4067  // do not count the candidate or the daughters of the candidate
4068  // in the isolation conte
4069  if ( aodParticle->GetDetectorTag() == kCTS ) // make sure conversions are tagged as kCTS!!!
4070  {
4071  Int_t trackID = GetReader()->GetTrackID(track) ; // needed instead of track->GetID() since AOD needs some manipulations
4072  Bool_t contained = kFALSE;
4073 
4074  for(Int_t i = 0; i < 4; i++)
4075  {
4076  if( trackID == aodParticle->GetTrackLabel(i) ) contained = kTRUE;
4077  }
4078 
4079  if ( contained ) continue ;
4080  }
4081 
4082  //Only for mixed event frame
4083  Int_t evtIndex2 = 0 ;
4084  if (GetMixedEvent())
4085  {
4086  evtIndex2 = GetMixedEvent()->EventIndex(GetReader()->GetTrackID(track)) ;
4087  if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
4088  continue ;
4089  //vertex cut
4090  if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
4091  continue;
4092  }
4093 
4094  AliDebug(2,Form("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
4095 
4096  // ------------------------------
4097  // Track type bin or bits setting
4098  //
4099 
4100  //
4101  // * Set the pt associated bin for the defined bins *
4102  //
4103  Int_t assocBin = -1;
4104  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
4105  {
4106  if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
4107  }
4108 
4109  //
4110  // * Assign to the histogram array a bin corresponding
4111  // to a combination of pTa and vz bins *
4112  //
4113  Int_t nz = 1;
4114  Int_t vz = 0;
4115 
4116  if(fCorrelVzBin)
4117  {
4118  nz = GetNZvertBin();
4119  vz = GetEventVzBin();
4120  }
4121 
4122  Int_t bin = assocBin*nz+vz;
4123 
4124  //printf("assoc Bin = %d, vZ bin = %d, bin = %d \n", assocBin,GetEventVzBin(),bin);
4125 
4126  //
4127  // * Get the status of the TOF bit *
4128  //
4129  ULong_t status = track->GetStatus();
4130  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
4131  //Double32_t tof = track->GetTOFsignal()*1e-3;
4132  Int_t trackBC = track->GetTOFBunchCrossing(bz);
4133 
4134  Int_t outTOF = -1;
4135  if (okTOF && trackBC!=0) outTOF = 1;
4136  else if(okTOF && trackBC==0) outTOF = 0;
4137 
4138  //----------------
4139  // Fill Histograms
4140 
4141  //
4142  // Azimuthal Angle histograms
4143  //
4144 
4145  deltaPhi = phiTrig-phi;
4146 
4147  //
4148  // Calculate deltaPhi shift so that for the particles on the opposite side
4149  // it is defined between 90 and 270 degrees
4150  // Shift [-360,-90] to [0, 270]
4151  // and [270,360] to [-90,0]
4152  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4153  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4154 
4155  FillChargedAngularCorrelationHistograms(pt, ptTrig, bin, phi, phiTrig, deltaPhi,
4156  eta, etaTrig, decayTag, track->GetHMPIDsignal(),
4157  outTOF, cenbin, mcTag);
4158 
4159  //
4160  // Imbalance zT/xE/pOut histograms
4161  //
4162 
4163  AliESDtrack * esdTrack = dynamic_cast<AliESDtrack*>(track);
4164 
4165  //
4166  // Delta phi cut for momentum imbalance correlation
4167  //
4168  if ( (deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) )
4169  {
4170  // Check track resolution
4171  if ( esdTrack )
4172  {
4173  fhTrackResolution->Fill(pt, TMath::Sqrt(esdTrack->GetCovariance()[14])*pt,
4174  GetEventWeight());
4175  }
4176 
4177  FillChargedMomentumImbalanceHistograms(ptTrig, pt, deltaPhi, cenbin, track->Charge(),
4178  assocBin, decayTag, outTOF, mcTag);
4179  }
4180 
4181  //
4182  // Underlying event, right side, default case
4183  //
4184  if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4185  {
4186  // Check track resolution
4187  if ( esdTrack )
4188  {
4189  fhTrackResolutionUE->Fill(pt, TMath::Sqrt(esdTrack->GetCovariance()[14])*pt,
4190  GetEventWeight());
4191  }
4192 
4193  FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, cenbin, outTOF,mcTag);
4194  }
4195 
4196  //
4197  // Several UE calculation, in different perpendicular regions, up to 6:
4198  // left, right, upper-left, lower left, upper-right, lower-right
4199  //
4200  if(fMakeSeveralUE)
4201  FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,deltaPhi,mcTag);
4202 
4203  //
4204  if(fPi0Trigger && decayFound)
4205  FillDecayPhotonCorrelationHistograms(pt, phi, kTRUE) ;
4206 
4207  //
4208  // Add track reference to array
4209  //
4210  if(fFillAODWithReferences)
4211  {
4212  nrefs++;
4213  if(nrefs==1)
4214  {
4215  reftracks = new TObjArray(0);
4216  TString trackname = Form("%sTracks", GetAODObjArrayName().Data());
4217  reftracks->SetName(trackname.Data());
4218  reftracks->SetOwner(kFALSE);
4219  }
4220 
4221  reftracks->Add(track);
4222  }// reference track to AOD
4223  }// track loop
4224 
4225  //Fill AOD with reference tracks, if not filling histograms
4226  if(fFillAODWithReferences && reftracks)
4227  {
4228  aodParticle->AddObjArray(reftracks);
4229  }
4230 }
4231 
4232 //_________________________________________________________________________________________________________
4234 //_________________________________________________________________________________________________________
4235 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
4236 {
4237  AliDebug(1,Form("Make trigger particle - charged hadron mixed event correlation"));
4238 
4239  if(GetMixedEvent()) return; // This is not the mixed event from general mixing frame
4240 
4241  // Get the event with similar caracteristics
4242  //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
4243 
4244  AliAnalysisManager * manager = AliAnalysisManager::GetAnalysisManager();
4245 
4246  AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
4247 
4248  if(!inputHandler) return;
4249 
4250  if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
4251 
4252  // Get the pool, check if it exits
4253  Int_t eventBin = GetEventMixBin();
4254 
4255  //Check that the bin exists, if not (bad determination of RP, ntrality or vz bin) do nothing
4256  if(eventBin < 0) return;
4257 
4258  fhEventBin->Fill(eventBin, GetEventWeight());
4259 
4260  // Get neutral clusters pool?
4261  Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
4262  Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
4263 
4264  TList * pool = 0;
4265  TList * poolCalo = 0;
4266  if(fUseMixStoredInReader)
4267  {
4268  pool = GetReader()->GetListWithMixedEventsForTracks(eventBin);
4269  if(neutralMix) poolCalo = GetReader()->GetListWithMixedEventsForCalo (eventBin);
4270  }
4271  else
4272  {
4273  pool = fListMixTrackEvents[eventBin];
4274  if(neutralMix) poolCalo = fListMixCaloEvents [eventBin];
4275  }
4276 
4277  if(!pool) return ;
4278 
4279  if( neutralMix && !poolCalo )
4280  AliWarning("Careful, cluster pool not available");
4281 
4282  Double_t ptTrig = aodParticle->Pt();
4283  Double_t etaTrig = aodParticle->Eta();
4284  Double_t phiTrig = aodParticle->Phi();
4285  if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
4286 
4287  AliDebug(1,Form("Pool bin %d size %d, trigger trigger pt=%f, phi=%f, eta=%f",
4288  eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig));
4289 
4290  Double_t ptAssoc = -999.;
4291  Double_t phiAssoc = -999.;
4292  Double_t etaAssoc = -999.;
4293  Double_t deltaPhi = -999.;
4294  Double_t deltaEta = -999.;
4295  Double_t xE = -999.;
4296 
4297  // Start from first event in pool except if in this same event the pool was filled
4298  Int_t ev0 = 0;
4299  if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
4300 
4301  for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
4302  {
4303  //
4304  // Recover the lists of tracks or clusters
4305  //
4306  TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
4307  TObjArray* bgCalo = 0;
4308 
4309  // Recover the clusters list if requested
4310  if( neutralMix && poolCalo )
4311  {
4312  if(pool->GetSize()!=poolCalo->GetSize())
4313  AliWarning("Different size of calo and track pools");
4314 
4315  bgCalo = static_cast<TObjArray*>(poolCalo->At(ev));
4316 
4317  if(!bgCalo) AliDebug(1,Form("Event %d in calo pool not available?",ev));
4318  }
4319 
4320  //
4321  // Isolate the trigger in the mixed event with mixed tracks and clusters
4322  //
4323  if( OnlyIsolated() )
4324  {
4325  Int_t n=0, nfrac = 0;
4326  Bool_t isolated = kFALSE;
4327  Float_t coneptsum = 0, coneptlead = 0;
4328  GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
4329  GetReader(), GetCaloPID(),
4330  kFALSE, aodParticle, "",
4331  n,nfrac,coneptsum,coneptlead,isolated);
4332 
4333  //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
4334  // isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
4335  //if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
4336  //printf("\n");
4337 
4338  if(!isolated) continue ;
4339  }
4340 
4341  //
4342  // Check if the trigger is leading of mixed event
4343  //
4344  Int_t nTracks=bgTracks->GetEntriesFast();
4345 
4346  if(fMakeNearSideLeading || fMakeAbsoluteLeading)
4347  {
4348  Bool_t leading = kTRUE;
4349  for(Int_t jlead = 0;jlead < nTracks; jlead++ )
4350  {
4351  AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(jlead) ;
4352 
4353  ptAssoc = track->Pt();
4354  phiAssoc = track->Phi() ;
4355  if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
4356 
4357  if (fMakeNearSideLeading)
4358  {
4359  deltaPhi = phiTrig-phiAssoc;
4360  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4361  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4362 
4363  if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
4364  {
4365  leading = kFALSE;
4366  break;
4367  }
4368  }
4369  //jump out this event if there is any other particle with pt larger than trigger
4370  else if(fMakeAbsoluteLeading)
4371  {
4372  if(ptAssoc > ptTrig)
4373  {
4374  leading = kFALSE;
4375  break;
4376  }
4377  }
4378  }
4379 
4380  if( !neutralMix && fCheckLeadingWithNeutralClusters )
4381  AliWarning("Leading of clusters requested but no clusters in mixed event");
4382 
4383  if(neutralMix && fCheckLeadingWithNeutralClusters && bgCalo)
4384  {
4385  Int_t nClusters=bgCalo->GetEntriesFast();
4386  for(Int_t jlead = 0;jlead <nClusters; jlead++ )
4387  {
4388  AliAODPWG4Particle *cluster= (AliAODPWG4Particle*) bgCalo->At(jlead) ;
4389 
4390  ptAssoc = cluster->Pt();
4391  phiAssoc = cluster->Phi() ;
4392  if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
4393 
4394  if (fMakeNearSideLeading)
4395  {
4396  deltaPhi = phiTrig-phiAssoc;
4397  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4398  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4399 
4400  if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
4401  {
4402  leading = kFALSE;
4403  break;
4404  }
4405  }
4406  //jump out this event if there is any other particle with pt larger than trigger
4407  else if(fMakeAbsoluteLeading)
4408  {
4409  if(ptAssoc > ptTrig)
4410  {
4411  leading = kFALSE;
4412  break;
4413  }
4414  }
4415  }
4416  }
4417 
4418  if(!leading) continue; // not leading, check the next event in pool
4419  }
4420 
4421  //
4422  // Fill histograms for selected triggers
4423  //
4424 
4425  fhEventMixBin->Fill(eventBin, GetEventWeight());
4426 
4427  //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
4428 
4429  fhPtTriggerMixed ->Fill(ptTrig, GetEventWeight());
4430  fhPhiTriggerMixed ->Fill(ptTrig, phiTrig, GetEventWeight());
4431  fhEtaTriggerMixed ->Fill(ptTrig, etaTrig, GetEventWeight());
4432  fhPtTriggerMixedBin->Fill(ptTrig, eventBin, GetEventWeight());
4433 
4434  if ( fCorrelVzBin ) fhPtTriggerMixedVzBin->Fill(ptTrig, GetEventVzBin(), GetEventWeight());
4435 
4436  //
4437  // Correlation histograms
4438  //
4439  for(Int_t j1 = 0;j1 <nTracks; j1++ )
4440  {
4441  AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
4442 
4443  if(!track) continue;
4444 
4445  ptAssoc = track->Pt();
4446  etaAssoc = track->Eta();
4447  phiAssoc = track->Phi() ;
4448  if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
4449 
4450  deltaPhi = phiTrig-phiAssoc;
4451  if(deltaPhi < -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4452  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4453  deltaEta = etaTrig-etaAssoc;
4454 
4455  AliDebug(1,Form("deltaPhi= %f, deltaEta=%f",deltaPhi, deltaEta));
4456 
4457  // Angular correlation
4458  fhMixDeltaPhiCharged ->Fill(ptTrig , deltaPhi, GetEventWeight());
4459  fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta, GetEventWeight());
4460 
4461  //
4462  // Momentum imbalance
4463  //
4464  if ( (deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) )
4465  {
4466  xE = -ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
4467 
4468  if( xE < 0.)
4469  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4470  xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4471 
4472  fhMixXECharged->Fill(ptTrig, xE, GetEventWeight());
4473  if ( xE > 0 ) fhMixHbpXECharged->Fill(ptTrig, TMath::Log(1./xE), GetEventWeight());
4474  }
4475 
4476  //
4477  // Underlying event momentum imbalance
4478  //
4479  if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4480  {
4481  //Underlying event region
4482  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
4483  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
4484 
4485  if(uexE < 0.)
4486  AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4487  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4488 
4489  fhMixXEUeCharged->Fill(ptTrig, uexE, GetEventWeight());
4490  }
4491 
4492  // Set the pt associated bin for the defined bins
4493  Int_t assocBin = -1;
4494  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
4495  {
4496  if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
4497  }
4498 
4499  //
4500  // Assign to the histogram array a bin corresponding to a combination of pTa and vz bins
4501  //
4502  Int_t nz = 1;
4503  Int_t vz = 0;
4504 
4505  if(fCorrelVzBin)
4506  {
4507  nz = GetNZvertBin();
4508  vz = GetEventVzBin();
4509  }
4510 
4511  Int_t bin = assocBin*nz+vz;
4512 
4513  if(bin < 0) continue ; // this pt bin was not considered
4514 
4515  fhMixDeltaPhiChargedAssocPtBin [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
4516  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->Fill(deltaPhi, deltaEta, GetEventWeight());
4517 
4518  if(fFillEtaGapsHisto)
4519  {
4520  if(TMath::Abs(deltaEta) > 0.8)
4521  fhMixDeltaPhiChargedAssocPtBinDEta08 [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
4522  if(TMath::Abs(deltaEta) < 0.01)
4523  fhMixDeltaPhiChargedAssocPtBinDEta0 [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
4524  }
4525 
4526  } // track loop
4527  } // mixed event loop
4528 }
4529 
4530 //_______________________________________________________________________________________________________
4533 //_______________________________________________________________________________________________________
4534 void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle)
4535 {
4536  TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); // For the future, foresee more possible pi0 lists
4537  if(!pi0list) return ;
4538 
4539  Int_t npi0 = pi0list->GetEntriesFast();
4540  if(npi0 == 0) return ;
4541 
4542  AliDebug(1,Form("Particle - pi0 correlation, %d pi0's",npi0));
4543 
4544  Int_t evtIndex11 = 0 ;
4545  Int_t evtIndex12 = 0 ;
4546  if (GetMixedEvent())
4547  {
4548  evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
4549  evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
4550  }
4551 
4552  Float_t pt = -100. ;
4553  Float_t zT = -100. ;
4554  Float_t phi = -100. ;
4555  Float_t eta = -100. ;
4556  Float_t xE = -100. ;
4557  Float_t hbpXE= -100. ;
4558  Float_t hbpZT= -100. ;
4559 
4560  Float_t ptTrig = aodParticle->Pt();
4561  Float_t phiTrig = aodParticle->Phi();
4562  Float_t etaTrig = aodParticle->Eta();
4563  Float_t deltaPhi= -100. ;
4564  Float_t deltaEta= -100. ;
4565 
4566  // In case of pi0/eta trigger, we may want to check their decay correlation,
4567  // get their decay children
4568 
4569  Bool_t decayFound = kFALSE;
4570  if(fPi0Trigger) decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
4571 
4572  TObjArray * refpi0 = 0x0;
4573  Int_t nrefs = 0;
4574 
4575  // Loop on stored AOD pi0
4576 
4577  for(Int_t iaod = 0; iaod < npi0 ; iaod++)
4578  {
4579  AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (pi0list->At(iaod));
4580 
4581  Int_t evtIndex2 = 0 ;
4582  Int_t evtIndex3 = 0 ;
4583  if (GetMixedEvent())
4584  {
4585  evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
4586  evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
4587 
4588  if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
4589  evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
4590  continue ;
4591  }
4592 
4593  pt = pi0->Pt();
4594 
4595  if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
4596 
4597  // Remove trigger itself for correlation when use charged triggers
4598  if(aodParticle->GetCaloLabel(0) >= 0 &&
4599  (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(0) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(0))) continue ;
4600 
4601  if( aodParticle->GetCaloLabel(1) >= 0 &&
4602  (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(1) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(1))) continue ;
4603 
4604  //
4605  // Angular correlations
4606  //
4607  phi = pi0->Phi() ;
4608  eta = pi0->Eta() ;
4609  deltaEta = etaTrig-eta;
4610  deltaPhi = phiTrig-phi;
4611  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4612  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4613 
4614  fhEtaNeutral ->Fill(pt , eta , GetEventWeight());
4615  fhPhiNeutral ->Fill(pt , phi , GetEventWeight());
4616  fhDeltaEtaNeutral->Fill(ptTrig, deltaEta, GetEventWeight());
4617  fhDeltaPhiNeutral->Fill(ptTrig, deltaPhi, GetEventWeight());
4618 
4619  if ( pt > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta, GetEventWeight());
4620 
4621  //
4622  // Momentum imbalance
4623  //
4624  zT = pt/ptTrig ;
4625 
4626  hbpZT = -100;
4627  hbpXE = -100;
4628 
4629  if(zT > 0 ) hbpZT = TMath::Log(1./zT);
4630 
4631  // Delta phi cut for correlation
4632  if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
4633  {
4634  xE =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
4635 
4636  if( xE < 0. )
4637  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4638  xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4639 
4640  if( xE > 0 ) hbpXE = TMath::Log(1./xE);
4641 
4642  fhDeltaPhiNeutralPt->Fill(pt , deltaPhi, GetEventWeight());
4643  fhXENeutral ->Fill(ptTrig, xE , GetEventWeight());
4644  fhPtHbpXENeutral ->Fill(ptTrig, hbpXE , GetEventWeight());
4645  fhZTNeutral ->Fill(ptTrig, zT , GetEventWeight());
4646  fhPtHbpZTNeutral ->Fill(ptTrig, hbpZT , GetEventWeight());
4647  }
4648  else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4649  {
4650  // Randomize angle for xE calculation
4651  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
4652 
4653  xE = -(pt/ptTrig)*TMath::Cos(randomphi);
4654  if(xE > 0 ) hbpXE = TMath::Log(1./xE);
4655 
4656  fhDeltaPhiUeNeutralPt->Fill(pt , deltaPhi, GetEventWeight());
4657  fhZTUeNeutral ->Fill(ptTrig, zT , GetEventWeight());
4658  fhPtHbpZTUeNeutral ->Fill(ptTrig, hbpZT , GetEventWeight());
4659  fhXEUeNeutral ->Fill(ptTrig, xE , GetEventWeight());
4660  fhPtHbpXEUeNeutral ->Fill(ptTrig, hbpXE , GetEventWeight());
4661  }
4662 
4663  // Several UE calculation, not sure it is useful
4664  // with partical calorimter acceptance
4665  if(fMakeSeveralUE) FillNeutralUnderlyingEventSidesHistograms(ptTrig,pt,zT,hbpZT,deltaPhi);
4666 
4667  //
4668  // Decay photon correlations
4669  //
4670  if(fPi0Trigger && decayFound)
4671  FillDecayPhotonCorrelationHistograms(pt, phi, kFALSE) ;
4672 
4673  if(fFillAODWithReferences)
4674  {
4675  nrefs++;
4676  if(nrefs==1)
4677  {
4678  refpi0 = new TObjArray(0);
4679  refpi0->SetName(GetAODObjArrayName()+"Pi0s");
4680  refpi0->SetOwner(kFALSE);
4681  }
4682  refpi0->Add(pi0);
4683  } // put references in trigger AOD
4684 
4685  AliDebug(1,Form("Selected pi0: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
4686 
4687  } // loop
4688 
4689  // Fill AOD with reference tracks, if not filling histograms
4690  if(fFillAODWithReferences && refpi0)
4691  {
4692  aodParticle->AddObjArray(refpi0);
4693  }
4694 }
4695 
4696 //__________________________________________________________________________________________________________________
4698 //__________________________________________________________________________________________________________________
4699 void AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int_t histoIndex, Bool_t lostDecayPair)
4700 {
4701  AliDebug(1,"Make trigger particle - charged hadron correlation in AOD MC level");
4702 
4703  if( label < 0 )
4704  {
4705  AliDebug(1,Form(" *** bad label ***: label %d", label));
4706  return;
4707  }
4708 
4709  // Do MC correlation for a given particle type range.
4710  // Types defined in GetMCTagHistogramIndex:
4711  // 0 direct gamma; 1 pi0; 2 pi0 decay; 3 eta decay; 4 other decay; 5 electron; 6 other (hadron)
4712  if(histoIndex < fMCGenTypeMin || histoIndex > fMCGenTypeMax) return ;
4713 
4714  AliStack * stack = 0x0 ;
4715  TParticle * primary = 0x0 ;
4716  TClonesArray * mcparticles = 0x0 ;
4717  AliAODMCParticle * aodprimary = 0x0 ;
4718 
4719  Double_t eprim = 0 ;
4720  Double_t ptprim = 0 ;
4721  Double_t phiprim = 0 ;
4722  Double_t etaprim = 0 ;
4723  Int_t nTracks = 0 ;
4724  Int_t iParticle = 0 ;
4725 
4726  Bool_t leadTrig = kTRUE;
4727 
4728  if( GetReader()->ReadStack() )
4729  {
4730  stack = GetMCStack() ;
4731  if(!stack)
4732  {
4733  AliFatal("Stack not available, is the MC handler called? STOP");
4734  return;
4735  }
4736 
4737  //nTracks = stack->GetNtrack() ;
4738  nTracks = stack->GetNprimary();
4739  if( label >= stack->GetNtrack() )
4740  {
4741  if(GetDebug() > 2)
4742  AliInfo(Form("*** large label ***: label %d, n tracks %d", label, stack->GetNtrack()));
4743  return ;
4744  }
4745 
4746  primary = stack->Particle(label);
4747  if ( !primary )
4748  {
4749  AliInfo(Form(" *** no primary ***: label %d", label));
4750  return;
4751  }
4752 
4753  eprim = primary->Energy();
4754  ptprim = primary->Pt();
4755  etaprim = primary->Eta();
4756  phiprim = primary->Phi();
4757  if(phiprim < 0) phiprim+=TMath::TwoPi();
4758 
4759  if(ptprim < 0.01 || eprim < 0.01) return ;
4760 
4761  for (iParticle = 0 ; iParticle < nTracks ; iParticle++)
4762  {
4763  if ( !GetReader()->AcceptParticleMCLabel( iParticle ) ) continue ;
4764 
4765  TParticle * particle = stack->Particle(iParticle);
4766 
4767  //keep only final state particles
4768  if( particle->GetStatusCode() != 1 ) continue ;
4769 
4770  //---------- Charged particles ----------------------
4771  Int_t pdg = particle->GetPdgCode();
4772  Int_t charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
4773  if(charge == 0) continue;
4774 
4775  particle->Momentum(fMomentum);
4776 
4777  //Particles in CTS acceptance, make sure to use the same selection as in the reader
4778  Bool_t inCTS = GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4779  //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4780  if( !inCTS ) continue;
4781 
4782  // Remove conversions
4783  if ( TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode() == 22 ) continue ;
4784 
4785  if ( label == iParticle ) continue; // avoid trigger particle
4786 
4787  Float_t phi = particle->Phi();
4788  if(phi < 0) phi+=TMath::TwoPi();
4789 
4790  Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim,histoIndex,lostDecayPair);
4791 
4792  if ( !lead ) leadTrig = kFALSE;
4793 
4794  // if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading) ) return;
4795  } // track loop
4796  } // ESD MC
4797  else if( GetReader()->ReadAODMCParticles() )
4798  {
4799  // Get the list of MC particles
4800  mcparticles = GetReader()->GetAODMCParticles();
4801  if( !mcparticles ) return;
4802 
4803  nTracks = mcparticles->GetEntriesFast() ;
4804 
4805  if( label >= nTracks )
4806  {
4807  if(GetDebug() > 2)
4808  AliInfo(Form(" *** large label ***: label %d, n tracks %d", label,nTracks));
4809  return;
4810  }
4811 
4812  // Get the particle
4813  aodprimary = (AliAODMCParticle*) mcparticles->At(label);
4814  if( !aodprimary )
4815  {
4816  AliInfo(Form(" *** no AOD primary ***: label %d", label));
4817  return;
4818  }
4819 
4820  eprim = aodprimary->E();
4821  ptprim = aodprimary->Pt();
4822  etaprim = aodprimary->Eta();
4823  phiprim = aodprimary->Phi();
4824  if(phiprim < 0) phiprim+=TMath::TwoPi();
4825 
4826  if(ptprim < 0.01 || eprim < 0.01) return ;
4827 
4828  for (iParticle = 0; iParticle < nTracks; iParticle++)
4829  {
4830  if ( !GetReader()->AcceptParticleMCLabel( iParticle ) ) continue ;
4831 
4832  AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(iParticle);
4833 
4834  if (!part->IsPhysicalPrimary() ) continue; // same as part->GetStatus() !=1
4835 
4836  if ( part->Charge() == 0 ) continue;
4837 
4838  fMomentum.SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->E());
4839 
4840  //Particles in CTS acceptance, make sure to use the same selection as in the reader
4841  Bool_t inCTS = GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4842  //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4843  if( !inCTS ) continue;
4844 
4845  // Remove conversions
4846  Int_t indexmother = part->GetMother();
4847  if ( indexmother > -1 )
4848  {
4849  Int_t pdg = part->GetPdgCode();
4850  Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
4851  if (TMath::Abs(pdg) == 11 && mPdg == 22) continue;
4852  }
4853 
4854  if