AliPhysics  97a96ce (97a96ce)
 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}" , "e^{#pm}" , "hadrons?" , "#gamma_{#pi decay} lost companion", "#gamma_{#eta decay} lost companion"} ;
1497  TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1498 
1499  TString parTitle = Form("#it{R} = %2.2f",GetIsolationCut()->GetConeSize());
1500 
1501  // For vz dependent histograms, if option ON
1502  Int_t nz = 1 ;
1503  if(fCorrelVzBin) nz = GetNZvertBin();
1504  TString sz = "" ;
1505  TString tz = "" ;
1506 
1507  // Fill histograms for neutral clusters in mixing?
1508  Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
1509  Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
1510 
1511  if( GetReader()->GetDataType() == AliCaloTrackReader::kESD )
1512  {
1513  fhTrackResolution = new TH2F ("hTrackResolution","Track resolution: #sigma_{#it{p}_{T}} vs #it{p}_{T}, away side, ESDs",
1514  nptbins,ptmin,ptmax,600,0,0.3);
1515  fhTrackResolution->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1516  fhTrackResolution->SetYTitle("#sigma_{#it{p}_{T}} / #it{p}_{T} (GeV/#it{c})");
1517  outputContainer->Add(fhTrackResolution);
1518 
1519  fhTrackResolutionUE = new TH2F ("hTrackResolutionUE","Track resolution: #sigma_{#it{p}_{T}} vs #it{p}_{T}, UE, ESDs",
1520  nptbins,ptmin,ptmax,600,0,0.3);
1521  fhTrackResolutionUE->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1522  fhTrackResolutionUE->SetYTitle("#sigma_{#it{p}_{T}} / #it{p}_{T} (GeV/#it{c})");
1523  outputContainer->Add(fhTrackResolutionUE);
1524  }
1525 
1526  fhPtTriggerInput = new TH1F("hPtTriggerInput","Input trigger #it{p}_{T}", nptbins,ptmin,ptmax);
1527  fhPtTriggerInput->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1528  outputContainer->Add(fhPtTriggerInput);
1529 
1530  if( fM02MaxCut > 0 && fM02MinCut > 0 )
1531  {
1532  fhPtTriggerSSCut = new TH1F("hPtTriggerSSCut","Trigger #it{p}_{T} after #lambda^{2}_{0} cut", nptbins,ptmin,ptmax);
1533  fhPtTriggerSSCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1534  outputContainer->Add(fhPtTriggerSSCut);
1535  }
1536 
1537  if( OnlyIsolated() )
1538  {
1539  fhPtTriggerIsoCut = new TH1F("hPtTriggerIsoCut","Trigger #it{p}_{T} after isolation (and #lambda^{2}_{0} cut)", nptbins,ptmin,ptmax);
1540  fhPtTriggerIsoCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1541  outputContainer->Add(fhPtTriggerIsoCut);
1542  }
1543 
1544  fhPtTriggerFidCut = new TH1F("hPtTriggerFidCut","Trigger #it{p}_{T} after fiducial (isolation and #lambda^{2}_{0}) cut", nptbins,ptmin,ptmax);
1545  fhPtTriggerFidCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1546  outputContainer->Add(fhPtTriggerFidCut);
1547 
1548  fhPtTrigger = new TH1F("hPtTrigger","#it{p}_{T} distribution of trigger particles (after opposite hadron leading cut and rest)", nptbins,ptmin,ptmax);
1549  fhPtTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1550  outputContainer->Add(fhPtTrigger);
1551 
1552  if(fFillInvMassHisto)
1553  {
1554  fhMassPtTrigger = new TH2F("hMassPtTrigger","2 photons invariant mass vs p_{T}^{trig}",
1555  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1556  fhMassPtTrigger->SetYTitle("M_{#gamma#gamma} (GeV/#it{c}^{2})");
1557  fhMassPtTrigger->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1558  outputContainer->Add(fhMassPtTrigger);
1559  }
1560 
1561  if(fFillBkgBinsHisto)
1562  {
1563 
1564  fhPtLeadInConeBin = new TH1F*[fNBkgBin] ;
1565 
1566  fhPtSumInConeBin = new TH1F*[fNBkgBin] ;
1567 
1568  fhPtLeadConeBinDecay = new TH1F*[fNBkgBin*fNDecayBits] ;
1569 
1570  fhSumPtConeBinDecay = new TH1F*[fNBkgBin*fNDecayBits] ;
1571 
1572  fhPtLeadConeBinMC = new TH1F*[fNBkgBin*fgkNmcTypes] ;
1573 
1574  fhSumPtConeBinMC = new TH1F*[fNBkgBin*fgkNmcTypes] ;
1575 
1576  for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1577  {
1578  fhPtLeadInConeBin[ibin] = new TH1F
1579  (Form("hPtLeadCone_Bin%d",ibin),
1580  Form("cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1581  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1582  fhPtLeadInConeBin[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
1583  fhPtLeadInConeBin[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1584  outputContainer->Add(fhPtLeadInConeBin[ibin]) ;
1585 
1586  fhPtSumInConeBin[ibin] = new TH1F
1587  (Form("hSumPtCone_Bin%d",ibin),
1588  Form("in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, %s",
1589  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1590  fhPtSumInConeBin[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
1591  fhPtSumInConeBin[ibin]->SetXTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
1592  outputContainer->Add(fhPtSumInConeBin[ibin]) ;
1593 
1594  if(fFillTaggedDecayHistograms)
1595  {
1596  for(Int_t idecay = 0; idecay < fNDecayBits; idecay++)
1597  {
1598  Int_t bindecay = ibin+idecay*fNBkgBin;
1599 
1600  fhPtLeadConeBinDecay[bindecay] = new TH1F
1601  (Form("hPtLeadCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
1602  Form("Decay bit %d, cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1603  fDecayBits[idecay],fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1604  fhPtLeadConeBinDecay[bindecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
1605  fhPtLeadConeBinDecay[bindecay]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1606  outputContainer->Add(fhPtLeadConeBinDecay[bindecay]) ;
1607 
1608  fhSumPtConeBinDecay[bindecay] = new TH1F
1609  (Form("hSumPtCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
1610  Form("Decay bit %d, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, %s",
1611  fDecayBits[idecay],fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1612  fhSumPtConeBinDecay[bindecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
1613  fhSumPtConeBinDecay[bindecay]->SetXTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
1614  outputContainer->Add(fhSumPtConeBinDecay[bindecay]) ;
1615  }
1616  }
1617 
1618  if(IsDataMC())
1619  {
1620  for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1621  {
1622  Int_t binmc = ibin+imc*fNBkgBin;
1623  fhPtLeadConeBinMC[binmc] = new TH1F
1624  (Form("hPtLeadCone_Bin%d_MC%s",ibin, nameMC[imc].Data()),
1625  Form("in cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, MC %s, %s",
1626  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
1627  fhPtLeadConeBinMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
1628  fhPtLeadConeBinMC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1629  outputContainer->Add(fhPtLeadConeBinMC[binmc]) ;
1630 
1631  fhSumPtConeBinMC[binmc] = new TH1F
1632  (Form("hSumPtCone_Bin%d_MC%s",ibin,nameMC[imc].Data()),
1633  Form("in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, MC %s, %s",
1634  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
1635  fhSumPtConeBinMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
1636  fhSumPtConeBinMC[binmc]->SetXTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
1637  outputContainer->Add(fhSumPtConeBinMC[binmc]) ;
1638  } // MC particle loop
1639  }
1640  }
1641  }
1642 
1643  if(IsDataMC())
1644  {
1645  for(Int_t i=0; i < fgkNmcTypes; i++)
1646  {
1647  fhPtTriggerMC[i] = new TH1F(Form("hPtTrigger_MC%s",nameMC[i].Data()),
1648  Form("#it{p}_{T} distribution of trigger particles, trigger origin is %s",nameMC[i].Data()),
1649  nptbins,ptmin,ptmax);
1650  fhPtTriggerMC[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1651  outputContainer->Add(fhPtTriggerMC[i]);
1652 
1653  if(fFillInvMassHisto)
1654  {
1655  fhMCMassPtTrigger[i] = new TH2F(Form("hMassPtTrigger_MC%s",nameMC[i].Data()),
1656  Form("2 photons invariant mass, trigger origin is %s",nameMC[i].Data()),
1657  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1658  fhMCMassPtTrigger[i]->SetYTitle("m_{#gamma#gamma}");
1659  fhMCMassPtTrigger[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1660  outputContainer->Add(fhMCMassPtTrigger[i]) ;
1661  }
1662  }
1663  }
1664 
1665  if(fDecayTrigger)
1666  {
1667  for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1668  {
1669  fhPtDecayTrigger[ibit] = new TH1F(Form("hPtDecayTrigger_bit%d",fDecayBits[ibit]),
1670  Form("#it{p}_{T} distribution of trigger particles, decay Bit %d",fDecayBits[ibit]),
1671  nptbins,ptmin,ptmax);
1672  fhPtDecayTrigger[ibit]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1673  outputContainer->Add(fhPtDecayTrigger[ibit]);
1674 
1675  if(IsDataMC())
1676  {
1677  for(Int_t i=0; i < fgkNmcTypes; i++)
1678  {
1679  fhPtDecayTriggerMC[ibit][i] = new TH1F(Form("hPtDecayTrigger_bit%d_MC%s",fDecayBits[ibit], nameMC[i].Data()),
1680  Form("#it{p}_{T} distribution of trigger particles, decay Bit %d, trigger origin is %s",fDecayBits[ibit], nameMC[i].Data()),
1681  nptbins,ptmin,ptmax);
1682  fhPtDecayTriggerMC[ibit][i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1683  outputContainer->Add(fhPtDecayTriggerMC[ibit][i]);
1684  }
1685  }
1686  }
1687  }
1688 
1689  if(fCorrelVzBin)
1690  {
1691  fhPtTriggerVzBin = new TH2F("hPtTriggerVzBin","#it{p}_{T} distribution of trigger particles vs vz bin", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
1692  fhPtTriggerVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1693  fhPtTriggerVzBin->SetYTitle("#it{v}_{#it{z}} bin");
1694  outputContainer->Add(fhPtTriggerVzBin);
1695  }
1696 
1697  fhPtTriggerBin = new TH2F ("hPtTriggerBin","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
1698  fhPtTriggerBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1699  fhPtTriggerBin->SetYTitle("Bin");
1700  outputContainer->Add(fhPtTriggerBin);
1701 
1702  fhPhiTrigger = new TH2F ("hPhiTrigger","#phi distribution of trigger Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1703  fhPhiTrigger->SetYTitle("#phi (rad)");
1704  outputContainer->Add(fhPhiTrigger);
1705 
1706  fhEtaTrigger = new TH2F ("hEtaTrigger","#eta distribution of trigger",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1707  fhEtaTrigger->SetYTitle("#eta ");
1708  outputContainer->Add(fhEtaTrigger);
1709 
1710  if(IsHighMultiplicityAnalysisOn())
1711  {
1712  fhPtTriggerCentrality = new TH2F("hPtTriggerCentrality","Trigger particle #it{p}_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
1713  fhPtTriggerCentrality->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1714  fhPtTriggerCentrality->SetYTitle("Centrality (%)");
1715  outputContainer->Add(fhPtTriggerCentrality) ;
1716 
1717  fhPtTriggerEventPlane = new TH2F("hPtTriggerEventPlane","Trigger particle #it{p}_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
1718  fhPtTriggerEventPlane->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1719  fhPtTriggerEventPlane->SetXTitle("EP angle (rad)");
1720  outputContainer->Add(fhPtTriggerEventPlane) ;
1721 
1722  fhTriggerEventPlaneCentrality = new TH2F("hTriggerEventPlaneCentrality","Trigger particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
1723  fhTriggerEventPlaneCentrality->SetXTitle("Centrality (%)");
1724  fhTriggerEventPlaneCentrality->SetYTitle("EP angle (rad)");
1725  outputContainer->Add(fhTriggerEventPlaneCentrality) ;
1726  }
1727 
1728  // Leading hadron in oposite side
1729  if(fFillLeadHadOppositeHisto)
1730  {
1731  fhPtLeadingOppositeHadron = new TH2F("hPtTriggerPtLeadingOppositeHadron","Leading hadron opposite to trigger vs trigger #it{p}_{T}",
1732  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1733  fhPtLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1734  fhPtLeadingOppositeHadron->SetYTitle("#it{p}_{T}^{lead hadron} (GeV/#it{c})");
1735  outputContainer->Add(fhPtLeadingOppositeHadron);
1736 
1737  fhPtNoLeadingOppositeHadron = new TH1F("hPtTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #it{p}_{T}",
1738  nptbins,ptmin,ptmax);
1739  fhPtNoLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1740  outputContainer->Add(fhPtNoLeadingOppositeHadron);
1741 
1742  fhEtaPhiNoLeadingOppositeHadron = new TH2F("hEtaPhiTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #eta:#phi",
1743  netabins,etamin,etamax,nphibins,phimin,phimax);
1744  fhEtaPhiNoLeadingOppositeHadron->SetXTitle("#eta");
1745  fhEtaPhiNoLeadingOppositeHadron->SetYTitle("#phi");
1746  outputContainer->Add(fhEtaPhiNoLeadingOppositeHadron);
1747 
1748 
1749  fhPtDiffPhiLeadingOppositeHadron = new TH2F("hPtTriggerDiffPhiTriggerLeadingOppositeHadron","#phi_{trigger}-#phi_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1750  nptbins,ptmin,ptmax,ndeltaphibins,deltaphimin,deltaphimax);
1751  fhPtDiffPhiLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1752  fhPtDiffPhiLeadingOppositeHadron->SetYTitle("#phi_{trigger}-#phi_{leading opposite hadron} (rad)");
1753  outputContainer->Add(fhPtDiffPhiLeadingOppositeHadron);
1754 
1755  fhPtDiffEtaLeadingOppositeHadron = new TH2F("hPtTriggerDiffEtaTriggerPhiLeadingOppositeHadron","#eta_{trigger}-#eta_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1756  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1757  fhPtDiffEtaLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1758  fhPtDiffEtaLeadingOppositeHadron->SetYTitle("#eta_{trigger}-#eta_{leading opposite hadron}");
1759  outputContainer->Add(fhPtDiffEtaLeadingOppositeHadron);
1760  }
1761 
1762  // Correlation with charged hadrons
1763 
1764  fhDeltaPhiDeltaEtaCharged = new TH2F
1765  ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
1766  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1767  fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
1768  fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1769 
1770  fhDeltaPhiDeltaEtaChargedPtA3GeV = new TH2F
1771  ("hDeltaPhiDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}, #it{p}_{TA}>3 GeV/#it{c}}",
1772  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1773  fhDeltaPhiDeltaEtaChargedPtA3GeV->SetXTitle("#Delta #phi (rad)");
1774  fhDeltaPhiDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1775 
1776  fhPhiCharged = new TH2F
1777  ("hPhiCharged","#phi_{h^{#pm}} vs #it{p}_{T #pm}",
1778  nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
1779  fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
1780  fhPhiCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1781 
1782  fhEtaCharged = new TH2F
1783  ("hEtaCharged","#eta_{h^{#pm}} vs #it{p}_{T #pm}",
1784  nptbins,ptmin,ptmax,100,-1.,1.);
1785  fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1786  fhEtaCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1787 
1788  fhDeltaPhiCharged = new TH2F
1789  ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
1790  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1791  fhDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
1792  fhDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1793 
1794  fhDeltaPhiChargedPtA3GeV = new TH2F
1795  ("hDeltaPhiChargedPtA3GeV","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1796  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1797  fhDeltaPhiChargedPtA3GeV->SetYTitle("#Delta #phi (rad)");
1798  fhDeltaPhiChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1799 
1800 
1801  fhDeltaPhiChargedPt = new TH2F
1802  ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",
1803  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1804  fhDeltaPhiChargedPt->SetYTitle("#Delta #phi (rad)");
1805  fhDeltaPhiChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1806 
1807  fhDeltaEtaCharged = new TH2F
1808  ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}",
1809  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1810  fhDeltaEtaCharged->SetYTitle("#Delta #eta");
1811  fhDeltaEtaCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1812 
1813  fhDeltaEtaChargedPtA3GeV = new TH2F
1814  ("hDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1815  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1816  fhDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1817  fhDeltaEtaChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1818 
1819  fhXECharged =
1820  new TH2F("hXECharged","#it{x}_{#it{E}} for charged tracks",
1821  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1822  fhXECharged->SetYTitle("#it{x}_{#it{E}}");
1823  fhXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1824 
1825  fhXECharged_Cone2 =
1826  new TH2F("hXECharged_Cone2","#it{x}_{#it{E}} for charged tracks in cone 2 (5#pi/6-7#pi/6)",
1827  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1828  fhXECharged_Cone2->SetYTitle("#it{x}_{#it{E}}");
1829  fhXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1830 
1831  fhXEPosCharged =
1832  new TH2F("hXEPositiveCharged","#it{x}_{#it{E}} for positive charged tracks",
1833  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1834  fhXEPosCharged->SetYTitle("#it{x}_{#it{E}}");
1835  fhXEPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1836 
1837  fhXENegCharged =
1838  new TH2F("hXENegativeCharged","#it{x}_{#it{E}} for negative charged tracks",
1839  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1840  fhXENegCharged->SetYTitle("#it{x}_{#it{E}}");
1841  fhXENegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1842 
1843  fhPtHbpXECharged =
1844  new TH2F("hHbpXECharged","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
1845  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1846  fhPtHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1847  fhPtHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1848 
1849  fhPtHbpXECharged_Cone2 =
1850  new TH2F("hHbpXECharged_Cone2","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons in cone 2 (5#pi/6-7#pi/6)",
1851  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1852  fhPtHbpXECharged_Cone2->SetYTitle("ln(1/#it{x}_{#it{E}})");
1853  fhPtHbpXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1854 
1855  fhZTCharged =
1856  new TH2F("hZTCharged","#it{z}_{T} for charged tracks",
1857  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1858  fhZTCharged->SetYTitle("#it{z}_{T}");
1859  fhZTCharged->SetXTitle("#it{p}_{T trigger}");
1860 
1861  fhZTPosCharged =
1862  new TH2F("hZTPositiveCharged","#it{z}_{T} for positive charged tracks",
1863  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1864  fhZTPosCharged->SetYTitle("#it{z}_{T}");
1865  fhZTPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1866 
1867  fhZTNegCharged =
1868  new TH2F("hZTNegativeCharged","#it{z}_{T} for negative charged tracks",
1869  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1870  fhZTNegCharged->SetYTitle("#it{z}_{T}");
1871  fhZTNegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1872 
1873  fhPtHbpZTCharged =
1874  new TH2F("hHbpZTCharged","#xi = ln(1/#it{z}_{T}) with charged hadrons",
1875  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1876  fhPtHbpZTCharged->SetYTitle("ln(1/#it{z}_{T})");
1877  fhPtHbpZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1878 
1879  fhPtTrigPout =
1880  new TH2F("hPtTrigPout","Pout with triggers",
1881  nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
1882  fhPtTrigPout->SetYTitle("#it{p}_{out} (GeV/#it{c})");
1883  fhPtTrigPout->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1884 
1885  fhPtTrigCharged =
1886  new TH2F("hPtTrigCharged","trigger and charged tracks pt distribution",
1887  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1888  fhPtTrigCharged->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1889  fhPtTrigCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1890 
1891  outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
1892  outputContainer->Add(fhDeltaPhiDeltaEtaChargedPtA3GeV);
1893  outputContainer->Add(fhPhiCharged) ;
1894  outputContainer->Add(fhEtaCharged) ;
1895  outputContainer->Add(fhDeltaPhiCharged) ;
1896  outputContainer->Add(fhDeltaPhiChargedPtA3GeV) ;
1897  outputContainer->Add(fhDeltaEtaCharged) ;
1898  outputContainer->Add(fhDeltaEtaChargedPtA3GeV) ;
1899  outputContainer->Add(fhDeltaPhiChargedPt) ;
1900 
1901  outputContainer->Add(fhXECharged) ;
1902  outputContainer->Add(fhXECharged_Cone2) ;
1903 
1904  if(IsDataMC())
1905  {
1906  for(Int_t i=0; i < fgkNmcTypes; i++)
1907  {
1908 
1909  fhDeltaPhiChargedMC[i] = new TH2F(Form("hDeltaPhiCharged_MC%s",nameMC[i].Data()),
1910  Form("#Delta #phi for charged tracks, trigger origin is %s",nameMC[i].Data()),
1911  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1912  fhDeltaPhiChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1913  fhDeltaPhiChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1914  outputContainer->Add(fhDeltaPhiChargedMC[i]) ;
1915 
1916  fhXEChargedMC[i] = new TH2F(Form("hXECharged_MC%s",nameMC[i].Data()),
1917  Form("#it{x}_{#it{E}} for charged tracks, trigger origin is %s",nameMC[i].Data()),
1918  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1919  fhXEChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1920  fhXEChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1921  outputContainer->Add(fhXEChargedMC[i]) ;
1922 
1923  fhXEUeChargedRightMC[i] = new TH2F(Form("hXEUeChargedRight_MC%s",nameMC[i].Data()),
1924  Form("#it{x}_{#it{E}} for charged tracks in right UE cone, trigger origin is %s",nameMC[i].Data()),
1925  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1926  fhXEUeChargedRightMC[i]->SetYTitle("#it{x}_{#it{E}}");
1927  fhXEUeChargedRightMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1928  outputContainer->Add(fhXEUeChargedRightMC[i]) ;
1929 
1930  fhXEUeChargedLeftMC[i] = new TH2F(Form("hXEUeChargedLeft_MC%s",nameMC[i].Data()),
1931  Form("#it{x}_{#it{E}} for charged tracks in left UE cone, trigger origin is %s",nameMC[i].Data()),
1932  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1933  fhXEUeChargedLeftMC[i]->SetYTitle("#it{x}_{#it{E}}");
1934  fhXEUeChargedLeftMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1935  outputContainer->Add(fhXEUeChargedLeftMC[i]) ;
1936  }
1937  }
1938 
1939  outputContainer->Add(fhXEPosCharged) ;
1940  outputContainer->Add(fhXENegCharged) ;
1941  outputContainer->Add(fhPtHbpXECharged) ;
1942  outputContainer->Add(fhPtHbpXECharged_Cone2) ;
1943 
1944  outputContainer->Add(fhZTCharged) ;
1945  outputContainer->Add(fhZTPosCharged) ;
1946  outputContainer->Add(fhZTNegCharged) ;
1947  outputContainer->Add(fhPtHbpZTCharged) ;
1948 
1949  outputContainer->Add(fhPtTrigPout) ;
1950  outputContainer->Add(fhPtTrigCharged) ;
1951 
1952  TString right = "";
1953  if(fMakeSeveralUE) right = "Right";
1954 
1955  fhUePart = new TH1F("hUePart","UE particles distribution vs pt trig",
1956  nptbins,ptmin,ptmax);
1957  fhUePart->SetYTitle("dNch");
1958  fhUePart->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1959 
1960  fhDeltaPhiUeChargedPt = new TH2F
1961  (Form("hDeltaPhiUe%sChargedPt",right.Data()),"#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}}",
1962  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1963  fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi (rad)");
1964  fhDeltaPhiUeChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1965 
1966  fhXEUeCharged =
1967  new TH2F(Form("hXEUeCharged%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event",
1968  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1969  fhXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
1970  fhXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1971 
1972  fhXEUeChargedSmallCone =
1973  new TH2F(Form("hXEUeChargedSmallCone%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event in a cone [80,100] deg",
1974  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1975  fhXEUeChargedSmallCone->SetYTitle("#it{x}_{#it{E}}");
1976  fhXEUeChargedSmallCone->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1977 
1978  fhXEUeChargedMediumCone =
1979  new TH2F(Form("hXEUeChargedMediumCone%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event in a cone [70,110] deg",
1980  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1981  fhXEUeChargedMediumCone->SetYTitle("#it{x}_{#it{E}}");
1982  fhXEUeChargedMediumCone->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1983 
1984  fhXEUeChargedLargeCone =
1985  new TH2F(Form("hXEUeChargedLargeCone%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event in a cone [60,120] deg",
1986  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1987  fhXEUeChargedLargeCone->SetYTitle("#it{x}_{#it{E}}");
1988  fhXEUeChargedLargeCone->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1989 
1990  fhPtHbpXEUeCharged =
1991  new TH2F(Form("hHbpXEUeCharged%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for Underlying Event",
1992  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1993  fhPtHbpXEUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1994  fhPtHbpXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1995 
1996  fhZTUeCharged =
1997  new TH2F(Form("hZTUeCharged%s",right.Data()),"#it{z}_{T} for Underlying Event",
1998  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1999  fhZTUeCharged->SetYTitle("#it{z}_{T}");
2000  fhZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2001 
2002  fhPtHbpZTUeCharged =
2003  new TH2F(Form("hHbpZTUeCharged%s",right.Data()),"#xi = ln(1/#it{z}_{T}) for Underlying Event",
2004  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2005  fhPtHbpZTUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2006  fhPtHbpZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2007 
2008  outputContainer->Add(fhUePart);
2009  outputContainer->Add(fhDeltaPhiUeChargedPt) ;
2010  outputContainer->Add(fhXEUeCharged) ;
2011  outputContainer->Add(fhXEUeChargedSmallCone) ;
2012  outputContainer->Add(fhXEUeChargedMediumCone) ;
2013  outputContainer->Add(fhXEUeChargedLargeCone) ;
2014  outputContainer->Add(fhPtHbpXEUeCharged) ;
2015  outputContainer->Add(fhZTUeCharged) ;
2016  outputContainer->Add(fhPtHbpZTUeCharged) ;
2017 
2018  if(fMakeSeveralUE)
2019  {
2020  fhDeltaPhiUeLeftCharged = new TH2F
2021  ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left side range of trigger particles",
2022  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2023  fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi (rad)");
2024  fhDeltaPhiUeLeftCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2025  outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
2026 
2027  fhDeltaPhiUeLeftUpCharged = new TH2F
2028  ("hDeltaPhiUeLeftUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Up side range of trigger particles",
2029  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2030  fhDeltaPhiUeLeftUpCharged->SetYTitle("#Delta #phi (rad)");
2031  fhDeltaPhiUeLeftUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2032  outputContainer->Add(fhDeltaPhiUeLeftUpCharged) ;
2033 
2034  fhDeltaPhiUeRightUpCharged = new TH2F
2035  ("hDeltaPhiUeRightUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Up side range of trigger particles",
2036  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2037  fhDeltaPhiUeRightUpCharged->SetYTitle("#Delta #phi (rad)");
2038  fhDeltaPhiUeRightUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2039  outputContainer->Add(fhDeltaPhiUeRightUpCharged) ;
2040 
2041  fhDeltaPhiUeLeftDownCharged = new TH2F
2042  ("hDeltaPhiUeLeftDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Down side range of trigger particles",
2043  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2044  fhDeltaPhiUeLeftDownCharged->SetYTitle("#Delta #phi (rad)");
2045  fhDeltaPhiUeLeftDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2046  outputContainer->Add(fhDeltaPhiUeLeftDownCharged) ;
2047 
2048  fhDeltaPhiUeRightDownCharged = new TH2F
2049  ("hDeltaPhiUeRightDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Down side range of trigger particles",
2050  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2051  fhDeltaPhiUeRightDownCharged->SetYTitle("#Delta #phi (rad)");
2052  fhDeltaPhiUeRightDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2053  outputContainer->Add(fhDeltaPhiUeRightDownCharged) ;
2054 
2055  fhXEUeLeftCharged =
2056  new TH2F("hXEUeChargedLeft","#it{x}_{#it{E}} with UE left side of trigger",
2057  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2058  fhXEUeLeftCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
2059  fhXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2060  outputContainer->Add(fhXEUeLeftCharged) ;
2061 
2062  fhXEUeLeftUpCharged =
2063  new TH2F("hXEUeChargedLeftUp","#it{x}_{#it{E}} with UE left Up side of trigger",
2064  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2065  fhXEUeLeftUpCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
2066  fhXEUeLeftUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2067  outputContainer->Add(fhXEUeLeftUpCharged) ;
2068 
2069  fhXEUeRightUpCharged =
2070  new TH2F("hXEUeChargedRightUp","#it{x}_{#it{E} h^{#pm}} with UE right Up side of trigger",
2071  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2072  fhXEUeRightUpCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
2073  fhXEUeRightUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2074  outputContainer->Add(fhXEUeRightUpCharged) ;
2075 
2076  fhXEUeLeftDownCharged =
2077  new TH2F("hXEUeChargedLeftDown","#it{x}_{#it{E}} with UE left Down side of trigger",
2078  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2079  fhXEUeLeftDownCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
2080  fhXEUeLeftDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2081  outputContainer->Add(fhXEUeLeftDownCharged) ;
2082 
2083  fhXEUeRightDownCharged =
2084  new TH2F("hXEUeChargedRightDown","#it{x}_{#it{E} h^{#pm}} with UE right Down side of trigger",
2085  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2086  fhXEUeRightDownCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
2087  fhXEUeRightDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2088  outputContainer->Add(fhXEUeRightDownCharged) ;
2089 
2090  fhPtHbpXEUeLeftCharged =
2091  new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/#it{x}_{#it{E}}) with charged UE left side of trigger",
2092  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2093  fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2094  fhPtHbpXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2095  outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
2096 
2097  fhZTUeLeftCharged =
2098  new TH2F("hZTUeChargedLeft","#it{z}_{trigger h^{#pm}} = #it{p}_{T Ueh^{#pm}} / #it{p}_{T trigger} with UE left side of trigger",
2099  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2100  fhZTUeLeftCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
2101  fhZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2102  outputContainer->Add(fhZTUeLeftCharged) ;
2103 
2104  fhPtHbpZTUeLeftCharged =
2105  new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/#it{z}_{T}) with charged UE left side of trigger",
2106  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2107  fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/#it{z}_{T})");
2108  fhPtHbpZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2109  outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
2110  }
2111 
2112  if(IsPileUpAnalysisOn())
2113  {
2114  fhDeltaPhiChargedOtherBC = new TH2F
2115  ("hDeltaPhiChargedOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC!=0",
2116  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2117  fhDeltaPhiChargedOtherBC->SetYTitle("#Delta #phi (rad)");
2118  fhDeltaPhiChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2119 
2120  fhDeltaPhiChargedPtA3GeVOtherBC = new TH2F
2121  ("hDeltaPhiChargedPtA3GeVOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC!=0",
2122  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2123  fhDeltaPhiChargedPtA3GeVOtherBC->SetYTitle("#Delta #phi (rad)");
2124  fhDeltaPhiChargedPtA3GeVOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2125 
2126  fhPtTrigChargedOtherBC =
2127  new TH2F("hPtTrigChargedOtherBC","trigger and charged tracks pt distribution, track BC!=0",
2128  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2129  fhPtTrigChargedOtherBC->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2130  fhPtTrigChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2131 
2132  fhXEChargedOtherBC =
2133  new TH2F("hXEChargedOtherBC","#it{x}_{#it{E}} for charged tracks, track BC!=0",
2134  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2135  fhXEChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
2136  fhXEChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2137 
2138  fhXEUeChargedOtherBC =
2139  new TH2F("hXEUeChargedOtherBC","#it{x}_{#it{E}} for Underlying Event, track BC!=0",
2140  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2141  fhXEUeChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
2142  fhXEUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2143 
2144  fhZTChargedOtherBC =
2145  new TH2F("hZTChargedOtherBC","#it{z}_{T} for charged tracks, track BC!=0",
2146  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2147  fhZTChargedOtherBC->SetYTitle("#it{z}_{T}");
2148  fhZTChargedOtherBC->SetXTitle("#it{p}_{T trigger}");
2149 
2150  fhZTUeChargedOtherBC =
2151  new TH2F("hZTUeChargedOtherBC","#it{z}_{T} for Underlying Event, track BC!=0",
2152  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2153  fhZTUeChargedOtherBC->SetYTitle("#it{z}_{T}");
2154  fhZTUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2155 
2156  outputContainer->Add(fhDeltaPhiChargedOtherBC) ;
2157  outputContainer->Add(fhDeltaPhiChargedPtA3GeVOtherBC) ;
2158  outputContainer->Add(fhXEChargedOtherBC) ;
2159  outputContainer->Add(fhXEUeChargedOtherBC) ;
2160  outputContainer->Add(fhZTChargedOtherBC) ;
2161  outputContainer->Add(fhZTUeChargedOtherBC) ;
2162  outputContainer->Add(fhPtTrigChargedOtherBC) ;
2163 
2164  fhDeltaPhiChargedBC0 = new TH2F
2165  ("hDeltaPhiChargedBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
2166  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2167  fhDeltaPhiChargedBC0->SetYTitle("#Delta #phi (rad)");
2168  fhDeltaPhiChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2169 
2170  fhDeltaPhiChargedPtA3GeVBC0 = new TH2F
2171  ("hDeltaPhiChargedPtA3GeVBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
2172  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2173  fhDeltaPhiChargedPtA3GeVBC0->SetYTitle("#Delta #phi (rad)");
2174  fhDeltaPhiChargedPtA3GeVBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2175 
2176  fhPtTrigChargedBC0 =
2177  new TH2F("hPtTrigChargedBC0","trigger and charged tracks pt distribution, track BC==0",
2178  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2179  fhPtTrigChargedBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2180  fhPtTrigChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2181 
2182  fhXEChargedBC0 =
2183  new TH2F("hXEChargedBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
2184  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2185  fhXEChargedBC0->SetYTitle("#it{x}_{#it{E}}");
2186  fhXEChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2187 
2188  fhXEUeChargedBC0 =
2189  new TH2F("hXEUeChargedBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
2190  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2191  fhXEUeChargedBC0->SetYTitle("#it{x}_{#it{E}}");
2192  fhXEUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2193 
2194  fhZTChargedBC0 =
2195  new TH2F("hZTChargedBC0","#it{z}_{T} for charged tracks, track BC==0",
2196  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2197  fhZTChargedBC0->SetYTitle("#it{z}_{T}");
2198  fhZTChargedBC0->SetXTitle("#it{p}_{T trigger}");
2199 
2200  fhZTUeChargedBC0 =
2201  new TH2F("hZTUeChargedBC0","#it{z}_{T} for Underlying Event, track BC==0",
2202  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2203  fhZTUeChargedBC0->SetYTitle("#it{z}_{T}");
2204  fhZTUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2205 
2206  outputContainer->Add(fhDeltaPhiChargedBC0) ;
2207  outputContainer->Add(fhDeltaPhiChargedPtA3GeVBC0) ;
2208  outputContainer->Add(fhXEChargedBC0) ;
2209  outputContainer->Add(fhXEUeChargedBC0) ;
2210  outputContainer->Add(fhZTChargedBC0) ;
2211  outputContainer->Add(fhZTUeChargedBC0) ;
2212  outputContainer->Add(fhPtTrigChargedBC0) ;
2213 
2214  fhPtTriggerVtxBC0 = new TH1F("hPtTriggerVtxBC0","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax);
2215  fhPtTriggerVtxBC0->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2216 
2217  fhDeltaPhiChargedVtxBC0 = new TH2F
2218  ("hDeltaPhiChargedVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
2219  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2220  fhDeltaPhiChargedVtxBC0->SetYTitle("#Delta #phi (rad)");
2221  fhDeltaPhiChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2222 
2223  fhDeltaPhiChargedPtA3GeVVtxBC0 = new TH2F
2224  ("hDeltaPhiChargedPtA3GeVVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
2225  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2226  fhDeltaPhiChargedPtA3GeVVtxBC0->SetYTitle("#Delta #phi (rad)");
2227  fhDeltaPhiChargedPtA3GeVVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2228 
2229  fhPtTrigChargedVtxBC0 =
2230  new TH2F("hPtTrigChargedVtxBC0","trigger and charged tracks pt distribution, track BC==0",
2231  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2232  fhPtTrigChargedVtxBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2233  fhPtTrigChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2234 
2235  fhXEChargedVtxBC0 =
2236  new TH2F("hXEChargedVtxBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
2237  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2238  fhXEChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
2239  fhXEChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2240 
2241  fhXEUeChargedVtxBC0 =
2242  new TH2F("hXEUeChargedVtxBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
2243  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2244  fhXEUeChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
2245  fhXEUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2246 
2247  fhZTChargedVtxBC0 =
2248  new TH2F("hZTChargedVtxBC0","#it{z}_{T} for charged tracks, track BC==0",
2249  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2250  fhZTChargedVtxBC0->SetYTitle("#it{z}_{T}");
2251  fhZTChargedVtxBC0->SetXTitle("#it{p}_{T trigger}");
2252 
2253  fhZTUeChargedVtxBC0 =
2254  new TH2F("hZTUeChargedVtxBC0","#it{z}_{T} for Underlying Event, track BC==0",
2255  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2256  fhZTUeChargedVtxBC0->SetYTitle("#it{z}_{T}");
2257  fhZTUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2258 
2259  outputContainer->Add(fhPtTriggerVtxBC0);
2260  outputContainer->Add(fhDeltaPhiChargedVtxBC0) ;
2261  outputContainer->Add(fhDeltaPhiChargedPtA3GeVVtxBC0) ;
2262  outputContainer->Add(fhXEChargedVtxBC0) ;
2263  outputContainer->Add(fhXEUeChargedVtxBC0) ;
2264  outputContainer->Add(fhZTChargedVtxBC0) ;
2265  outputContainer->Add(fhZTUeChargedVtxBC0) ;
2266  outputContainer->Add(fhPtTrigChargedVtxBC0) ;
2267 
2268  for(Int_t i = 0 ; i < 7 ; i++)
2269  {
2270  fhPtTriggerPileUp[i] = new TH1F(Form("hPtTriggerPileUp%s",pileUpName[i].Data()),
2271  Form("#it{p}_{T} distribution of trigger particles, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2272  fhPtTriggerPileUp[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2273  outputContainer->Add(fhPtTriggerPileUp[i]);
2274 
2275  fhDeltaPhiChargedPileUp[i] = new TH2F(Form("hDeltaPhiChargedPileUp%s",pileUpName[i].Data()),
2276  Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
2277  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2278  fhDeltaPhiChargedPileUp[i]->SetYTitle("#Delta #phi (rad)");
2279  fhDeltaPhiChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2280  outputContainer->Add(fhDeltaPhiChargedPileUp[i]) ;
2281 
2282  fhDeltaPhiChargedPtA3GeVPileUp[i] = new TH2F(Form("hDeltaPhiChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
2283  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()),
2284  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2285  fhDeltaPhiChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #phi (rad)");
2286  fhDeltaPhiChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2287  outputContainer->Add(fhDeltaPhiChargedPtA3GeVPileUp[i]) ;
2288 
2289  fhDeltaEtaChargedPileUp[i] = new TH2F(Form("hDeltaEtaChargedPileUp%s",pileUpName[i].Data()),
2290  Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
2291  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
2292  fhDeltaEtaChargedPileUp[i]->SetYTitle("#Delta #eta");
2293  fhDeltaEtaChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2294  outputContainer->Add(fhDeltaEtaChargedPileUp[i]) ;
2295 
2296  fhDeltaEtaChargedPtA3GeVPileUp[i] = new TH2F(Form("hDeltaEtaChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
2297  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()),
2298  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
2299  fhDeltaEtaChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #eta");
2300  fhDeltaEtaChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2301  outputContainer->Add(fhDeltaEtaChargedPtA3GeVPileUp[i]) ;
2302 
2303  fhXEChargedPileUp[i] = new TH2F(Form("hXEChargedPileUp%s",pileUpName[i].Data()),
2304  Form("#it{x}_{#it{E}} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
2305  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2306  fhXEChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
2307  fhXEChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2308  outputContainer->Add(fhXEChargedPileUp[i]) ;
2309 
2310  fhXEUeChargedPileUp[i] = new TH2F(Form("hXEUeChargedPileUp%s",pileUpName[i].Data()),
2311  Form("#it{x}_{#it{E}} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
2312  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2313  fhXEUeChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
2314  fhXEUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2315  outputContainer->Add(fhXEUeChargedPileUp[i]) ;
2316 
2317  fhZTChargedPileUp[i] = new TH2F(Form("hZTChargedPileUp%s",pileUpName[i].Data()),
2318  Form("#it{z}_{T} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
2319  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2320  fhZTChargedPileUp[i]->SetYTitle("#it{z}_{T}");
2321  fhZTChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2322  outputContainer->Add(fhZTChargedPileUp[i]) ;
2323 
2324  fhZTUeChargedPileUp[i] = new TH2F(Form("hZTUeChargedPileUp%s",pileUpName[i].Data()),
2325  Form("#it{z}_{T} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
2326  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2327  fhZTUeChargedPileUp[i]->SetYTitle("#it{z}_{T}");
2328  fhZTUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2329  outputContainer->Add(fhZTUeChargedPileUp[i]) ;
2330 
2331  fhPtTrigChargedPileUp[i] = new TH2F(Form("hPtTrigChargedPileUp%s",pileUpName[i].Data()),
2332  Form("trigger and charged tracks pt distribution, %s Pile-Up event",pileUpName[i].Data()),
2333  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2334  fhPtTrigChargedPileUp[i]->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2335  fhPtTrigChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2336  outputContainer->Add(fhPtTrigChargedPileUp[i]) ;
2337  }
2338  }
2339 
2340  if(IsHighMultiplicityAnalysisOn())
2341  {
2342  Int_t nMultiBins = GetNCentrBin();
2343  fhDeltaPhiChargedMult = new TH2F*[nMultiBins] ;
2344  fhDeltaEtaChargedMult = new TH2F*[nMultiBins] ;
2345  fhXEMult = new TH2F*[nMultiBins] ;
2346  fhXEUeMult = new TH2F*[nMultiBins] ;
2347  fhZTMult = new TH2F*[nMultiBins] ;
2348  fhZTUeMult = new TH2F*[nMultiBins] ;
2349 
2350  for(Int_t im=0; im<nMultiBins; im++)
2351  {
2352  fhDeltaPhiChargedMult[im] = new TH2F
2353  (Form("hDeltaPhiCharged_Mult%d",im),Form("#Delta #phi charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2354  fhDeltaPhiChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2355  fhDeltaPhiChargedMult[im]->SetYTitle("#Delta #phi (rad)");
2356 
2357  fhDeltaEtaChargedMult[im] = new TH2F
2358  (Form("hDeltaEtaCharged_Mult%d",im),Form("#Delta #eta charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2359  fhDeltaEtaChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2360  fhDeltaEtaChargedMult[im]->SetYTitle("#Delta #eta");
2361 
2362  fhXEMult[im] = new TH2F
2363  (Form("hXECharged_Mult%d",im),Form("#it{x}_{E} charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2364  fhXEMult[im]->SetYTitle("#it{x}_{E}");
2365  fhXEMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2366 
2367  fhXEUeMult[im] = new TH2F
2368  (Form("hXEUeCharged_Mult%d",im),Form("#it{x}_{E} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2369  fhXEUeMult[im]->SetYTitle("#it{x}_{E}");
2370  fhXEUeMult[im]->SetXTitle("#it{p}_{T trigger}(GeV/#it{c})");
2371 
2372  fhZTMult[im] = new TH2F
2373  (Form("hZTCharged_Mult%d",im),Form("#it{z}_{T} charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2374  fhZTMult[im]->SetYTitle("#it{z}_{T}");
2375  fhZTMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2376 
2377  fhZTUeMult[im] = new TH2F
2378  (Form("hZTUeCharged_Mult%d",im),Form("#it{z}_{T} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2379  fhZTUeMult[im]->SetYTitle("#it{z}_{T}");
2380  fhZTUeMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2381 
2382  outputContainer->Add(fhDeltaPhiChargedMult[im]) ;
2383  outputContainer->Add(fhDeltaEtaChargedMult[im]) ;
2384  outputContainer->Add(fhXEMult [im]);
2385  outputContainer->Add(fhXEUeMult[im]);
2386  outputContainer->Add(fhZTMult [im]);
2387  outputContainer->Add(fhZTUeMult[im]);
2388  }
2389  }
2390 
2391  if(fFillBradHisto)
2392  {
2393  fhAssocPtBkg = new TH2F("hAssocPtBkg", " Trigger #it{p}_{T} vs associated hadron #it{p}_{T} from background",
2394  nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
2395  fhAssocPtBkg->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2396  fhAssocPtBkg->SetYTitle("#it{p}_{T associated} (GeV/#it{c})");
2397  outputContainer->Add(fhAssocPtBkg) ;
2398 
2399  fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} ",
2400  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2401  fhDeltaPhiBrad->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2402  fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2403  outputContainer->Add(fhDeltaPhiBrad) ;
2404  }
2405 
2406  fhDeltaPhiDeltaEtaAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2407  fhDeltaPhiAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2408  fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
2409  fhDeltaPhiAssocPtBinDEta0 = new TH2F*[fNAssocPtBins*nz];
2410  if(fFillMomImbalancePtAssocBinsHisto)
2411  {
2412  fhXEAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2413  fhZTAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2414  }
2415 
2416  if(fCorrelVzBin)
2417  {
2418  fhXEVZ = new TH2F*[nz];
2419  fhZTVZ = new TH2F*[nz];
2420  }
2421 
2422  if(fFillBradHisto)
2423  fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2424 
2425 
2426  fhDeltaPhiAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2427  if(fFillEtaGapsHisto)fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
2428  if(fDecayTrigger) fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2429 
2430 
2431  if(fHMPIDCorrelation)
2432  {
2433  fhDeltaPhiAssocPtBinHMPID = new TH2F*[fNAssocPtBins*nz];
2434  fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins*nz];
2435  }
2436 
2437  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2438  {
2439  for(Int_t z = 0 ; z < nz ; z++)
2440  {
2441  Int_t bin = i*nz+z;
2442 
2443  if(fCorrelVzBin)
2444  {
2445  sz = Form("_vz%d",z);
2446  tz = Form(", #it{v}_{#it{z}} bin %d",z);
2447  }
2448 
2449  //printf("iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2450 
2451  fhDeltaPhiDeltaEtaAssocPtBin[bin] = new TH2F(Form("hDeltaPhiDeltaEtaPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2452  Form("#Delta #phi vs #Delta #eta for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2453  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
2454  fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2455  fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetYTitle("#Delta #eta");
2456 
2457  fhDeltaPhiAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2458  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()),
2459  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2460  fhDeltaPhiAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2461  fhDeltaPhiAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2462 
2463  outputContainer->Add(fhDeltaPhiDeltaEtaAssocPtBin[bin]) ;
2464  outputContainer->Add(fhDeltaPhiAssocPtBin[bin]) ;
2465 
2466  if(fFillEtaGapsHisto)
2467  {
2468  fhDeltaPhiAssocPtBinDEta08[bin] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2469  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()),
2470  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2471  fhDeltaPhiAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2472  fhDeltaPhiAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2473 
2474  fhDeltaPhiAssocPtBinDEta0[bin] = new TH2F(Form("hDeltaPhiDeltaEta0PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2475  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()),
2476  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2477  fhDeltaPhiAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2478  fhDeltaPhiAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2479 
2480  outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[bin]) ;
2481  outputContainer->Add(fhDeltaPhiAssocPtBinDEta0[bin]) ;
2482  }
2483 
2484  if(fDecayTrigger)
2485  {
2486  fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s_bit%d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data(),fDecayBits[0]),
2487  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]),
2488  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2489  fhDeltaPhiDecayChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2490  fhDeltaPhiDecayChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2491 
2492  outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[bin]) ;
2493  }
2494 
2495  if(fFillBradHisto)
2496  {
2497  fhDeltaPhiBradAssocPtBin[bin] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2498  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()),
2499  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2500  fhDeltaPhiBradAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2501  fhDeltaPhiBradAssocPtBin[bin]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2502  outputContainer->Add(fhDeltaPhiBradAssocPtBin[bin]) ;
2503  }
2504 
2505  if(fHMPIDCorrelation)
2506  {
2507  fhDeltaPhiAssocPtBinHMPID[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2508  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()),
2509  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2510  fhDeltaPhiAssocPtBinHMPID[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})" );
2511  fhDeltaPhiAssocPtBinHMPID[bin]->SetYTitle("#Delta #phi (rad)");
2512 
2513  fhDeltaPhiAssocPtBinHMPIDAcc[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2514  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()),
2515  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2516  fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2517  fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetYTitle("#Delta #phi (rad)");
2518 
2519  outputContainer->Add(fhDeltaPhiAssocPtBinHMPID [bin]) ;
2520  outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[bin]) ;
2521  }
2522  }
2523  }
2524 
2525  if(fFillMomImbalancePtAssocBinsHisto)
2526  {
2527  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2528  {
2529  fhXEAssocPtBin[i] = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2530  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]),
2531  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2532  fhXEAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2533  fhXEAssocPtBin[i]->SetYTitle("#it{x}_{#it{E}}");
2534 
2535  fhZTAssocPtBin[i] = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2536  Form("#it{z}_{T} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2537  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2538  fhZTAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2539  fhZTAssocPtBin[i]->SetYTitle("#it{z}_{T}");
2540 
2541 
2542  outputContainer->Add(fhXEAssocPtBin[i]);
2543  outputContainer->Add(fhZTAssocPtBin[i]);
2544  }
2545  }
2546 
2547  if(fCorrelVzBin)
2548  {
2549  for(Int_t z = 0 ; z < nz ; z++)
2550  {
2551  sz = Form("_vz%d",z);
2552  tz = Form(", #it{v}_{#it{z}} bin %d",z);
2553 
2554  fhXEVZ[z] = new TH2F(Form("hXE%s", sz.Data()),
2555  Form("#it{x}_{#it{E}} vs #it{p}_{T trigger}%s", tz.Data()),
2556  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2557  fhXEVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2558  fhXEVZ[z]->SetYTitle("#it{x}_{#it{E}}");
2559 
2560  fhZTVZ[z] = new TH2F(Form("hZT%s",sz.Data()),
2561  Form("#it{z}_{T} vs #it{p}_{T trigger}%s", tz.Data()),
2562  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2563  fhZTVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2564  fhZTVZ[z]->SetYTitle("#it{z}_{T}");
2565 
2566  outputContainer->Add(fhXEVZ[z]);
2567  outputContainer->Add(fhZTVZ[z]);
2568  }
2569  }
2570 
2571  if(fPi0Trigger)
2572  {
2573  fhPtPi0DecayRatio = new TH2F
2574  ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
2575  nptbins,ptmin,ptmax, 100,0.,2.);
2576  fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
2577  fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
2578  outputContainer->Add(fhPtPi0DecayRatio) ;
2579 
2580  fhDeltaPhiPi0DecayCharged = new TH2F
2581  ("hDeltaPhiPi0DecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
2582  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2583  fhDeltaPhiPi0DecayCharged->SetYTitle("#Delta #phi (rad)");
2584  fhDeltaPhiPi0DecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2585 
2586  fhXEPi0DecayCharged =
2587  new TH2F("hXEPi0DecayCharged","#it{x}_{#it{E}} Decay",
2588  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2589  fhXEPi0DecayCharged->SetYTitle("#it{x}_{#it{E}}");
2590  fhXEPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2591 
2592  fhZTPi0DecayCharged =
2593  new TH2F("hZTPi0DecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
2594  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2595  fhZTPi0DecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
2596  fhZTPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2597 
2598  outputContainer->Add(fhDeltaPhiPi0DecayCharged) ;
2599  outputContainer->Add(fhXEPi0DecayCharged) ;
2600  outputContainer->Add(fhZTPi0DecayCharged) ;
2601  }
2602 
2603  if(fDecayTrigger)
2604  {
2605  for(Int_t ibit = 0; ibit< fNDecayBits; ibit++)
2606  {
2607  fhDeltaPhiDecayCharged[ibit] = new TH2F
2608  (Form("hDeltaPhiDecayCharged_bit%d",fDecayBits[ibit]),
2609  Form("#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2610  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2611  fhDeltaPhiDecayCharged[ibit]->SetYTitle("#Delta #phi (rad)");
2612  fhDeltaPhiDecayCharged[ibit]->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2613 
2614  fhXEDecayCharged[ibit] =
2615  new TH2F(Form("hXEDecayCharged_bit%d",fDecayBits[ibit]),
2616  Form("#it{x}_{#it{E}} Decay, Bit %d",fDecayBits[ibit]),
2617  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2618  fhXEDecayCharged[ibit]->SetYTitle("#it{x}_{#it{E}}");
2619  fhXEDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2620 
2621  fhZTDecayCharged[ibit] =
2622  new TH2F(Form("hZTDecayCharged_bit%d",fDecayBits[ibit]),
2623  Form("#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2624  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2625  fhZTDecayCharged[ibit]->SetYTitle("#it{z}_{decay h^{#pm}}");
2626  fhZTDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2627 
2628  outputContainer->Add(fhDeltaPhiDecayCharged[ibit]) ;
2629  outputContainer->Add(fhXEDecayCharged[ibit]) ;
2630  outputContainer->Add(fhZTDecayCharged[ibit]) ;
2631  }
2632  }
2633 
2634  //Correlation with neutral hadrons
2635  if(fNeutralCorr)
2636  {
2637  fhDeltaPhiDeltaEtaNeutral = new TH2F
2638  ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
2639  ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
2640  fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi (rad)");
2641  fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
2642 
2643  fhPhiNeutral = new TH2F
2644  ("hPhiNeutral","#phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}",
2645  nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2646  fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
2647  fhPhiNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2648 
2649  fhEtaNeutral = new TH2F
2650  ("hEtaNeutral","#eta_{#pi^{0}} vs #it{p}_{T #pi^{0}}",
2651  nptbins,ptmin,ptmax,200,-1.,1.);
2652  fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
2653  fhEtaNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2654 
2655  fhDeltaPhiNeutral = new TH2F
2656  ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T trigger}",
2657  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2658  fhDeltaPhiNeutral->SetYTitle("#Delta #phi (rad)");
2659  fhDeltaPhiNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2660 
2661  fhDeltaPhiNeutralPt = new TH2F
2662  ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2663  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2664  fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi (rad)");
2665  fhDeltaPhiNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2666 
2667  fhDeltaEtaNeutral = new TH2F
2668  ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs #it{p}_{T trigger}",
2669  nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2670  fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
2671  fhDeltaEtaNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2672 
2673  fhXENeutral =
2674  new TH2F("hXENeutral","#it{x}_{#it{E}} for #pi^{0} associated",
2675  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2676  fhXENeutral->SetYTitle("#it{x}_{#it{E}}");
2677  fhXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2678 
2679  fhPtHbpXENeutral =
2680  new TH2F("hHbpXENeutral","#xi = ln(1/#it{x}_{#it{E}})for #pi^{0} associated",
2681  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2682  fhPtHbpXENeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2683  fhPtHbpXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2684 
2685  fhZTNeutral =
2686  new TH2F("hZTNeutral","#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2687  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2688  fhZTNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2689  fhZTNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2690 
2691  fhPtHbpZTNeutral =
2692  new TH2F("hHbpZTNeutral","#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2693  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2694  fhPtHbpZTNeutral->SetYTitle("ln(1/#it{z}_{T})");
2695  fhPtHbpZTNeutral->SetXTitle("#it{p}_{T trigger}");
2696 
2697  outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
2698  outputContainer->Add(fhPhiNeutral) ;
2699  outputContainer->Add(fhEtaNeutral) ;
2700  outputContainer->Add(fhDeltaPhiNeutral) ;
2701  outputContainer->Add(fhDeltaPhiNeutralPt) ;
2702  outputContainer->Add(fhDeltaEtaNeutral) ;
2703  outputContainer->Add(fhXENeutral) ;
2704  outputContainer->Add(fhPtHbpXENeutral) ;
2705  outputContainer->Add(fhZTNeutral) ;
2706  outputContainer->Add(fhPtHbpZTNeutral) ;
2707 
2708  fhDeltaPhiUeNeutralPt = new TH2F
2709  (Form("hDeltaPhiUe%sNeutralPt",right.Data()),"#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2710  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2711  fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi (rad)");
2712  fhDeltaPhiUeNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2713 
2714  fhXEUeNeutral =
2715  new TH2F(Form("hXEUeNeutral%s",right.Data()),"#it{x}_{#it{E}} for #pi^{0} associated",
2716  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2717  fhXEUeNeutral->SetYTitle("#it{x}_{#it{E}}");
2718  fhXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2719 
2720  fhPtHbpXEUeNeutral =
2721  new TH2F(Form("hHbpXEUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2722  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2723  fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2724  fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2725 
2726  fhZTUeNeutral =
2727  new TH2F(Form("hZTUeNeutral%s",right.Data()),"#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2728  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2729  fhZTUeNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2730  fhZTUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2731 
2732  fhPtHbpZTUeNeutral =
2733  new TH2F(Form("hHbpZTUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2734  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2735  fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{z}_{T})");
2736  fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2737 
2738  outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
2739  outputContainer->Add(fhXEUeNeutral) ;
2740  outputContainer->Add(fhPtHbpXEUeNeutral) ;
2741  outputContainer->Add(fhZTUeNeutral) ;
2742  outputContainer->Add(fhPtHbpZTUeNeutral) ;
2743 
2744  if(fMakeSeveralUE)
2745  {
2746  fhDeltaPhiUeLeftNeutral = new TH2F
2747  ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs #it{p}_{T h^{0}} with neutral UE left side range of trigger particles",
2748  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2749  fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi (rad)");
2750  fhDeltaPhiUeLeftNeutral->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2751  outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
2752 
2753  fhXEUeLeftNeutral =
2754  new TH2F("hXEUeNeutralLeft","#it{x}_{#it{E}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2755  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2756  fhXEUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2757  fhXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2758  outputContainer->Add(fhXEUeLeftNeutral) ;
2759 
2760  fhPtHbpXEUeLeftNeutral =
2761  new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/#it{x}_{#it{E}}) with neutral UE left side of trigger",
2762  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2763  fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2764  fhPtHbpXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2765  outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
2766 
2767  fhZTUeLeftNeutral =
2768  new TH2F("hZTUeNeutralLeft","#it{z}_{trigger h^{0}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2769  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2770  fhZTUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2771  fhZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2772  outputContainer->Add(fhZTUeLeftNeutral) ;
2773 
2774  fhPtHbpZTUeLeftNeutral =
2775  new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/#it{z}_{T}) with neutral UE left side of trigger",
2776  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2777  fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/#it{z}_{T})");
2778  fhPtHbpZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger}");
2779  outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
2780  }
2781 
2782  if(fPi0Trigger)
2783  {
2784  fhDeltaPhiPi0DecayNeutral = new TH2F
2785  ("hDeltaPhiPi0DecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
2786  nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2787  fhDeltaPhiPi0DecayNeutral->SetYTitle("#Delta #phi (rad)");
2788  fhDeltaPhiPi0DecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2789 
2790  fhXEPi0DecayNeutral =
2791  new TH2F("hXEPi0DecayNeutral","#it{x}_{#it{E}} for decay trigger",
2792  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2793  fhXEPi0DecayNeutral->SetYTitle("#it{x}_{#it{E}}");
2794  fhXEPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2795 
2796  fhZTPi0DecayNeutral =
2797  new TH2F("hZTPi0DecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
2798  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2799  fhZTPi0DecayNeutral->SetYTitle("#it{z}_{h^{0}}");
2800  fhZTPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2801 
2802  outputContainer->Add(fhDeltaPhiPi0DecayNeutral) ;
2803  outputContainer->Add(fhXEPi0DecayNeutral) ;
2804  outputContainer->Add(fhZTPi0DecayNeutral) ;
2805  }
2806  }//Correlation with neutral hadrons
2807 
2808  // If data is MC, fill more histograms, depending on origin
2809  if(IsDataMC())
2810  {
2811  for(Int_t i= fMCGenTypeMin; i <= fMCGenTypeMax; i++)
2812  {
2813  fhMCPtTrigger[i] = new TH1F (Form("hMCPtTrigger_%s",nameMC[i].Data()),
2814  Form("MC %s: generated trigger #it{p}_{T}",nameMC[i].Data()),
2815  nptbins,ptmin,ptmax);
2816  fhMCPtTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2817 
2818  fhMCPhiTrigger[i] = new TH2F (Form("hMCPhiTrigger_%s",nameMC[i].Data()),
2819  Form("MC %s: generated trigger #phi",nameMC[i].Data()),
2820  nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2821  fhMCPhiTrigger[i]->SetYTitle("#phi (rad)");
2822  fhMCPhiTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2823 
2824  fhMCEtaTrigger[i] = new TH2F (Form("hMCEtaTrigger_%s",nameMC[i].Data()),
2825  Form("MC %s: generated trigger #eta",nameMC[i].Data()),
2826  nptbins,ptmin,ptmax, netabins,etamin,etamax);
2827  fhMCEtaTrigger[i]->SetYTitle("#eta");
2828  fhMCEtaTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2829 
2830  if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2831  {
2832  fhMCPtTriggerNotLeading[i] = new TH1F (Form("hMCPtTriggerNotLeading_%s",nameMC[i].Data()),
2833  Form("MC %s: generated trigger #it{p}_{T}, when not leading of primaries",nameMC[i].Data()),
2834  nptbins,ptmin,ptmax);
2835  fhMCPtTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2836 
2837  fhMCPhiTriggerNotLeading[i] = new TH2F (Form("hMCPhiTriggerNotLeading_%s",nameMC[i].Data()),
2838  Form("MC %s: generated trigger #phi, when not leading of primaries",nameMC[i].Data()),
2839  nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2840  fhMCPhiTriggerNotLeading[i]->SetYTitle("#phi (rad)");
2841  fhMCPhiTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2842 
2843 
2844  fhMCEtaTriggerNotLeading[i] = new TH2F (Form("hMCEtaTriggerNotLeading_%s",nameMC[i].Data()),
2845  Form("MC %s: generated triogger #eta, when not leading of primaries",nameMC[i].Data()),
2846  nptbins,ptmin,ptmax, netabins,etamin,etamax);
2847  fhMCEtaTriggerNotLeading[i]->SetYTitle("#eta ");
2848  fhMCEtaTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2849  }
2850 
2851  fhMCEtaCharged[i] = new TH2F (Form("hMCEtaCharged_%s",nameMC[i].Data()),
2852  Form("MC %s: #eta_{h^{#pm}} vs #it{p}_{T #pm}",nameMC[i].Data()),
2853  nptbins,ptmin,ptmax,100,-1.,1.);
2854  fhMCEtaCharged[i]->SetYTitle("#eta_{h^{#pm}} (rad)");
2855  fhMCEtaCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2856 
2857  fhMCPhiCharged[i] = new TH2F(Form("hMCPhiCharged_%s",nameMC[i].Data()),
2858  Form("MC %s: phi_{h^{#pm}} vs #it{p}_{T #pm}",nameMC[i].Data()),
2859  nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2860  fhMCPhiCharged[i]->SetYTitle("MC #phi_{h^{#pm}} (rad)");
2861  fhMCPhiCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2862 
2863  fhMCDeltaPhiDeltaEtaCharged[i] = new TH2F (Form("hMCDeltaPhiDeltaEtaCharged_%s",nameMC[i].Data()),
2864  Form("MC %s: phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",nameMC[i].Data()),
2865  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2866  fhMCDeltaPhiDeltaEtaCharged[i]->SetXTitle("#Delta #phi (rad)");
2867  fhMCDeltaPhiDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2868 
2869  fhMCDeltaEtaCharged[i] = new TH2F (Form("hMCDeltaEtaCharged_%s",nameMC[i].Data()),
2870  Form("MC %s: #eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger} and #it{p}_{T assoc}",nameMC[i].Data()),
2871  nptbins,ptmin,ptmax,ndeltaetabins ,deltaetamin,deltaetamax);
2872  fhMCDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2873  fhMCDeltaEtaCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2874 
2875  fhMCDeltaPhiCharged[i] = new TH2F (Form("hMCDeltaPhiCharged_%s",nameMC[i].Data()),
2876  Form("MC %s: #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",nameMC[i].Data()),
2877  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2878  fhMCDeltaPhiCharged[i]->SetYTitle("#Delta #phi (rad)");
2879  fhMCDeltaPhiCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2880 
2881  fhMCDeltaPhiChargedPt[i] = new TH2F (Form("hMCDeltaPhiChargedPt_%s",nameMC[i].Data()),
2882  Form("MC %s: #phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",nameMC[i].Data()),
2883  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2884  fhMCDeltaPhiChargedPt[i]->SetYTitle("#Delta #phi (rad)");
2885  fhMCDeltaPhiChargedPt[i]->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2886 
2887  fhMCPtXECharged[i] = new TH2F (Form("hMCPtXECharged_%s",nameMC[i].Data()),
2888  Form("MC %s: #it{x}_{#it{E}} with charged hadrons",nameMC[i].Data()),
2889  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2890  fhMCPtXECharged[i]->SetYTitle("#it{x}_{#it{E}}");
2891  fhMCPtXECharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2892 
2893  fhMCPtHbpXECharged[i] = new TH2F(Form("hMCHbpXECharged_%s",nameMC[i].Data()),
2894  Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",nameMC[i].Data()),
2895  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2896  fhMCPtHbpXECharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2897  fhMCPtHbpXECharged[i]->SetXTitle("#it{p}_{T trigger}");
2898 
2899  fhMCPtZTCharged[i] = new TH2F(Form("hMCPtZTCharged_%s",nameMC[i].Data()),
2900  Form("MC %s: #it{z}_{T} with charged hadrons",nameMC[i].Data()),
2901  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2902  fhMCPtZTCharged[i]->SetYTitle("#it{z}_{T}");
2903  fhMCPtZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2904 
2905  fhMCPtHbpZTCharged[i] = new TH2F(Form("hMCHbpZTCharged_%s",nameMC[i].Data()),
2906  Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons",nameMC[i].Data()),
2907  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2908  fhMCPtHbpZTCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2909  fhMCPtHbpZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2910 
2911  fhMCPtTrigPout[i] = new TH2F(Form("hMCPtTrigPout_%s",nameMC[i].Data()),
2912  Form("MC %s: #it{p}_{out} with triggers",nameMC[i].Data()),
2913  nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
2914  fhMCPtTrigPout[i]->SetYTitle("#it{p}_{out} (GeV/#it{c})");
2915  fhMCPtTrigPout[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2916 
2917  fhMCPtAssocDeltaPhi[i] = new TH2F(Form("hMCPtAssocDeltaPhi_%s",nameMC[i].Data()),
2918  Form("MC %s: #Delta #phi with associated charged hadrons",nameMC[i].Data()),
2919  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2920  fhMCPtAssocDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
2921  fhMCPtAssocDeltaPhi[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2922 
2923  outputContainer->Add(fhMCPtTrigger[i]);
2924  outputContainer->Add(fhMCPhiTrigger[i]);
2925  outputContainer->Add(fhMCEtaTrigger[i]);
2926 
2927  if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2928  {
2929  outputContainer->Add(fhMCPtTriggerNotLeading[i]);
2930  outputContainer->Add(fhMCPhiTriggerNotLeading[i]);
2931  outputContainer->Add(fhMCEtaTriggerNotLeading[i]);
2932  }
2933 
2934  outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged[i]);
2935  outputContainer->Add(fhMCPhiCharged[i]) ;
2936  outputContainer->Add(fhMCEtaCharged[i]) ;
2937  outputContainer->Add(fhMCDeltaEtaCharged[i]) ;
2938  outputContainer->Add(fhMCDeltaPhiCharged[i]) ;
2939 
2940  outputContainer->Add(fhMCDeltaPhiChargedPt[i]) ;
2941  outputContainer->Add(fhMCPtXECharged[i]) ;
2942  outputContainer->Add(fhMCPtZTCharged[i]) ;
2943  outputContainer->Add(fhMCPtHbpXECharged[i]) ;
2944  outputContainer->Add(fhMCPtHbpZTCharged[i]) ;
2945  outputContainer->Add(fhMCPtTrigPout[i]) ;
2946  outputContainer->Add(fhMCPtAssocDeltaPhi[i]) ;
2947 
2948  // Underlying event
2949 
2950  fhMCUePart[i] =
2951  new TH1F(Form("hMCUePart_%s",nameMC[i].Data()),
2952  Form("MC %s: UE particles distribution vs #it{p}_{T trigger}",nameMC[i].Data()),
2953  nptbins,ptmin,ptmax);
2954  fhMCUePart[i]->SetYTitle("#it{dN}^{ch}");
2955  fhMCUePart[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2956 
2957  fhMCPtXEUeCharged[i] =
2958  new TH2F(Form("hMCPtXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2959  Form("MC %s: #it{x}_{#it{E}} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2960  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2961  fhMCPtXEUeCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2962  fhMCPtXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2963 
2964  fhMCPtHbpXEUeCharged[i] =
2965  new TH2F(Form("hMCPtHbpXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2966  Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2967  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2968  fhMCPtHbpXEUeCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2969  fhMCPtHbpXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2970 
2971  fhMCPtZTUeCharged[i] =
2972  new TH2F(Form("hMCPtZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2973  Form("MC %s: #it{z}_{T} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2974  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2975  fhMCPtZTUeCharged[i]->SetYTitle("#it{z}_{T}");
2976  fhMCPtZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2977 
2978  fhMCPtHbpZTUeCharged[i] =
2979  new TH2F(Form("hMCPtHbpZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2980  Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2981  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2982  fhMCPtHbpZTUeCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2983  fhMCPtHbpZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2984 
2985  outputContainer->Add(fhMCUePart[i]);
2986  outputContainer->Add(fhMCPtXEUeCharged[i]) ;
2987  outputContainer->Add(fhMCPtZTUeCharged[i]) ;
2988  outputContainer->Add(fhMCPtHbpZTUeCharged[i]);
2989  outputContainer->Add(fhMCPtHbpXEUeCharged[i]);
2990 
2991  if(fMakeSeveralUE)
2992  {
2993  fhMCPtXEUeLeftCharged[i] = new TH2F(Form("hMCPtXEUeChargedLeft_%s",nameMC[i].Data()),
2994  Form("MC %s: #it{x}_{#it{E}} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2995  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2996  fhMCPtXEUeLeftCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2997  fhMCPtXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2998 
2999  fhMCPtHbpXEUeLeftCharged[i] = new TH2F(Form("hMCPtHbpXEUeChargedLeft_%s",nameMC[i].Data()),
3000  Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
3001  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
3002  fhMCPtHbpXEUeLeftCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
3003  fhMCPtHbpXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3004 
3005  fhMCPtZTUeLeftCharged[i] = new TH2F(Form("hMCPtZTUeChargedLeft_%s",nameMC[i].Data()),
3006  Form("MC %s: #it{z}_{T} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
3007  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
3008  fhMCPtZTUeLeftCharged[i]->SetYTitle("#it{z}_{T}");
3009  fhMCPtZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3010 
3011  fhMCPtHbpZTUeLeftCharged[i] = new TH2F(Form("hMCPtHbpZTUeChargedLeft_%s",nameMC[i].Data()),
3012  Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
3013  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
3014  fhMCPtHbpZTUeLeftCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
3015  fhMCPtHbpZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3016 
3017  outputContainer->Add(fhMCPtXEUeLeftCharged[i]) ;
3018  outputContainer->Add(fhMCPtZTUeLeftCharged[i]) ;
3019  outputContainer->Add(fhMCPtHbpXEUeLeftCharged[i]);
3020  outputContainer->Add(fhMCPtHbpZTUeLeftCharged[i]) ;
3021 
3022  }
3023  }
3024  } // For MC histogram
3025 
3026  if(DoOwnMix())
3027  {
3028  // Create event containers
3029 
3030  if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
3031  {
3032  Int_t nvz = GetNZvertBin();
3033  Int_t nrp = GetNRPBin();
3034  Int_t nce = GetNCentrBin();
3035 
3036  fListMixTrackEvents= new TList*[nvz*nrp*nce] ;
3037 
3038  for( Int_t ice = 0 ; ice < nce ; ice++ )
3039  {
3040  for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
3041  {
3042  for( Int_t irp = 0 ; irp < nrp ; irp++ )
3043  {
3044  Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
3045 
3046  //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
3047  // ic,iz, irp, bin);
3048 
3049  fListMixTrackEvents[bin] = new TList() ;
3050  fListMixTrackEvents[bin]->SetOwner(kFALSE);
3051  }
3052  }
3053  }
3054  }
3055 
3056  fhPtTriggerMixed = new TH1F ("hPtTriggerMixed","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax);
3057  fhPtTriggerMixed->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
3058 
3059  if(fCorrelVzBin)
3060  {
3061  fhPtTriggerMixedVzBin = new TH2F ("hPtTriggerMixedVzBin","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
3062  fhPtTriggerMixedVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
3063  fhPtTriggerMixedVzBin->SetYTitle("#it{v}_{#it{z}} bin");
3064  outputContainer->Add(fhPtTriggerMixedVzBin);
3065  }
3066 
3067  fhPtTriggerMixedBin = new TH2F ("hPtTriggerMixedBin","#it{p}_{T} distribution of trigger particles vs mixing bin", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
3068  fhPtTriggerMixedBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
3069  fhPtTriggerMixedBin->SetYTitle("Bin");
3070 
3071  fhPhiTriggerMixed = new TH2F ("hPhiTriggerMixed","#phi distribution of trigger Particles, used for mixing",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
3072  fhPhiTriggerMixed->SetYTitle("#phi (rad)");
3073 
3074  fhEtaTriggerMixed = new TH2F ("hEtaTriggerMixed","#eta distribution of trigger, used for mixing",nptbins,ptmin,ptmax, netabins,etamin,etamax);
3075  fhEtaTriggerMixed->SetYTitle("#eta ");
3076 
3077  outputContainer->Add(fhPtTriggerMixed);
3078  outputContainer->Add(fhPtTriggerMixedBin);
3079  outputContainer->Add(fhPhiTriggerMixed);
3080  outputContainer->Add(fhEtaTriggerMixed);
3081 
3082  // Fill the cluster pool only in isolation analysis or if requested
3083  if( neutralMix && (!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForCaloExists())))
3084  {
3085  Int_t nvz = GetNZvertBin();
3086  Int_t nrp = GetNRPBin();
3087  Int_t nce = GetNCentrBin();
3088 
3089  fListMixCaloEvents= new TList*[nvz*nrp*nce] ;
3090 
3091  for( Int_t ice = 0 ; ice < nce ; ice++ )
3092  {
3093  for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
3094  {
3095  for( Int_t irp = 0 ; irp < nrp ; irp++ )
3096  {
3097  Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
3098 
3099  //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
3100  // ic,iz, irp, bin);
3101 
3102  fListMixCaloEvents[bin] = new TList() ;
3103  fListMixCaloEvents[bin]->SetOwner(kFALSE);
3104  }
3105  }
3106  }
3107  }
3108 
3109  // Init the list in the reader if not done previously
3110  if(fUseMixStoredInReader)
3111  {
3112  if( !GetReader()->ListWithMixedEventsForTracksExists() )
3113  GetReader()->SetListWithMixedEventsForTracks(fListMixTrackEvents);
3114 
3115  if( !GetReader()->ListWithMixedEventsForCaloExists() )
3116  GetReader()->SetListWithMixedEventsForCalo (fListMixCaloEvents );
3117  }
3118 
3119  fhEventBin=new TH1I("hEventBin","Number of triggers per bin(cen,vz,rp)",
3120  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3121  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3122  fhEventBin->SetXTitle("event bin");
3123  outputContainer->Add(fhEventBin) ;
3124 
3125  fhEventMixBin=new TH1I("hEventMixBin","Number of triggers mixed per event bin(cen,vz,rp)",
3126  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3127  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3128  fhEventMixBin->SetXTitle("event bin");
3129  outputContainer->Add(fhEventMixBin) ;
3130 
3131  fhEventMBBin=new TH1I("hEventMBBin","Number of min bias events per bin(cen,vz,rp)",
3132  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3133  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3134  fhEventMBBin->SetXTitle("event bin");
3135  outputContainer->Add(fhEventMBBin) ;
3136 
3137  fhNtracksMB=new TH2F("hNtracksMBEvent","Number of filtered tracks in MB event per event bin",ntrbins,trmin,trmax,
3138  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3139  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3140  fhNtracksMB->SetYTitle("event bin");
3141  fhNtracksMB->SetXTitle("#it{N}_{track}");
3142  outputContainer->Add(fhNtracksMB);
3143 
3144  if( neutralMix )
3145  {
3146  fhNclustersMB=new TH2F("hNclustersMBEvent","Number of filtered clusters in MB events per event bin",nclbins,clmin,clmax,
3147  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
3148  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
3149  fhNclustersMB->SetYTitle("event bin");
3150  fhNclustersMB->SetXTitle("#it{N}_{cluster}");
3151  outputContainer->Add(fhNclustersMB);
3152  }
3153 
3154  fhMixDeltaPhiCharged = new TH2F
3155  ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
3156  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
3157  fhMixDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
3158  fhMixDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3159  outputContainer->Add(fhMixDeltaPhiCharged);
3160 
3161  fhMixDeltaPhiDeltaEtaCharged = new TH2F
3162  ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
3163  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
3164  fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
3165  fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
3166  outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
3167 
3168  fhMixXECharged =
3169  new TH2F("hMixXECharged","Mixed event : #it{x}_{#it{E}} for charged tracks",
3170  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
3171  fhMixXECharged->SetYTitle("#it{x}_{#it{E}}");
3172  fhMixXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3173  outputContainer->Add(fhMixXECharged);
3174 
3175  fhMixXEUeCharged =
3176  new TH2F("hMixXEUeCharged","Mixed event : #it{x}_{#it{E}} for charged tracks in Ue region",
3177  nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
3178  fhMixXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
3179  fhMixXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3180  outputContainer->Add(fhMixXEUeCharged);
3181 
3182  fhMixHbpXECharged =
3183  new TH2F("hMixHbpXECharged","mixed event : #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
3184  nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
3185  fhMixHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
3186  fhMixHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3187  outputContainer->Add(fhMixHbpXECharged);
3188 
3189  fhMixDeltaPhiChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
3190  fhMixDeltaPhiChargedAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
3191  fhMixDeltaPhiChargedAssocPtBinDEta0 = new TH2F*[fNAssocPtBins*nz];
3192  fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
3193 
3194  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3195  {
3196  for(Int_t z = 0 ; z < nz ; z++)
3197  {
3198  Int_t bin = i*nz+z;
3199 
3200  if(fCorrelVzBin)
3201  {
3202  sz = Form("_vz%d",z);
3203  tz = Form(", #it{v}_{#it{z}} bin %d",z);
3204  }
3205 
3206  //printf("MIX : iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
3207 
3208  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3209  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()),
3210  ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
3211  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
3212  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetYTitle("#Delta #eta");
3213 
3214  outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]);
3215 
3216  fhMixDeltaPhiChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3217  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()),
3218  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
3219  fhMixDeltaPhiChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3220  fhMixDeltaPhiChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
3221 
3222  outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[bin]);
3223 
3224  if(fFillEtaGapsHisto)
3225  {
3226  fhMixDeltaPhiChargedAssocPtBinDEta08[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3227  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()),
3228  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
3229  fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3230  fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
3231 
3232  fhMixDeltaPhiChargedAssocPtBinDEta0[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
3233  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()),
3234  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
3235  fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
3236  fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
3237 
3238  outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[bin]);
3239  outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta0[bin]);
3240  }
3241  }
3242  }
3243  }
3244 
3245  return outputContainer;
3246 }
3247 
3248 //_____________________________________________________________________________________________________________________
3256 //_____________________________________________________________________________________________________________________
3257 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(Int_t indexPhoton1, Int_t indexPhoton2, Int_t idetector)
3258 {
3259  if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
3260 
3261  AliDebug(1,Form("indexPhoton1 = %d, indexPhoton2 = %d", indexPhoton1, indexPhoton2));
3262 
3263  TObjArray * clusters = 0x0 ;
3264  if(idetector==kEMCAL) clusters = GetEMCALClusters() ;
3265  else clusters = GetPHOSClusters() ;
3266 
3267  for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
3268  {
3269  AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
3270 
3271  if(photon->GetID()==indexPhoton1) photon->GetMomentum(fDecayMom1,GetVertex(0)) ;
3272  if(photon->GetID()==indexPhoton2) photon->GetMomentum(fDecayMom2,GetVertex(0)) ;
3273 
3274  AliDebug(1,Form("Photon1 = %f, Photon2 = %f", fDecayMom1.Pt(), fDecayMom2.Pt()));
3275  } //cluster loop
3276 
3277  return kTRUE;
3278 }
3279 
3280 //________________________________________________________________________
3283 //________________________________________________________________________
3284 Int_t AliAnaParticleHadronCorrelation::GetMCTagHistogramIndex(Int_t mcTag)
3285 {
3286  if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||
3287  GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation) ||
3288  GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCISR) ) return 0;
3289  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ) return 1;
3290  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay) ) return 2;
3291  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ) return 3;
3292  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay) ) return 4;
3293  else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) ) return 5; // other decays
3294  else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) ) return 6;
3295  else return 7;
3296 }
3297 
3298 //_________________________________________
3300 //_________________________________________
3301 void AliAnaParticleHadronCorrelation::Init()
3302 {
3303  if(!GetReader()->IsCTSSwitchedOn())
3304  AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!");
3305 }
3306 
3307 //____________________________________________________
3309 //____________________________________________________
3310 void AliAnaParticleHadronCorrelation::InitParameters()
3311 {
3312  SetInputAODName("Particle");
3313  SetAODObjArrayName("Hadrons");
3314  AddToHistogramsName("AnaHadronCorr_");
3315 
3316  SetPtCutRange(0.,300);
3317  fDeltaPhiMinCut = TMath::DegToRad()*120.;
3318  fDeltaPhiMaxCut = TMath::DegToRad()*240. ;
3319  fSelectIsolated = kFALSE;
3320  fMakeSeveralUE = kFALSE;
3321  fUeDeltaPhiMinCut = TMath::DegToRad()*60.;
3322  fUeDeltaPhiMaxCut = TMath::DegToRad()*120 ;
3323 
3324  fNeutralCorr = kFALSE ;
3325  fPi0Trigger = kFALSE ;
3326  fDecayTrigger = kFALSE ;
3327  fHMPIDCorrelation = kFALSE ;
3328 
3329  fMakeAbsoluteLeading = kTRUE;
3330  fMakeNearSideLeading = kFALSE;
3331 
3332  fNAssocPtBins = 9 ;
3333  fAssocPtBinLimit[0] = 0.2 ;
3334  fAssocPtBinLimit[1] = 0.5 ;
3335  fAssocPtBinLimit[2] = 1.0 ;
3336  fAssocPtBinLimit[3] = 2.0 ;
3337  fAssocPtBinLimit[4] = 3.0 ;
3338  fAssocPtBinLimit[5] = 4.0 ;
3339  fAssocPtBinLimit[6] = 5.0 ;
3340  fAssocPtBinLimit[7] = 6.0 ;
3341  fAssocPtBinLimit[8] = 7.0 ;
3342  fAssocPtBinLimit[9] = 8.0 ;
3343  fAssocPtBinLimit[10] = 9.0 ;
3344  fAssocPtBinLimit[11] = 10.0 ;
3345  fAssocPtBinLimit[12] = 12.0 ;
3346  fAssocPtBinLimit[13] = 14.0 ;
3347  fAssocPtBinLimit[14] = 16.0 ;
3348  fAssocPtBinLimit[15] = 20.0 ;
3349  fAssocPtBinLimit[16] = 30.0 ;
3350  fAssocPtBinLimit[17] = 40.0 ;
3351  fAssocPtBinLimit[18] = 50.0 ;
3352  fAssocPtBinLimit[19] = 200.0 ;
3353 
3354  fUseMixStoredInReader = kTRUE;
3355 
3356  fM02MinCut = -1 ;
3357  fM02MaxCut = -1 ;
3358 
3359  fDecayTagsM02Cut = 0.27;
3360 
3361  fSelectLeadingHadronAngle = kFALSE;
3362  fFillLeadHadOppositeHisto = kFALSE;
3363  fMinLeadHadPhi = 150*TMath::DegToRad();
3364  fMaxLeadHadPhi = 210*TMath::DegToRad();
3365 
3366  fMinLeadHadPt = 1;
3367  fMaxLeadHadPt = 100;
3368 
3369  fMCGenTypeMin = 0;
3370  fMCGenTypeMax = 10;
3371 
3372  fNDecayBits = 1;
3373  fDecayBits[0] = AliNeutralMesonSelection::kPi0;
3374  fDecayBits[1] = AliNeutralMesonSelection::kEta;
3375  fDecayBits[2] = AliNeutralMesonSelection::kPi0RightSide;
3376  fDecayBits[3] = AliNeutralMesonSelection::kEtaRightSide;
3377  fDecayBits[4] = AliNeutralMesonSelection::kEtaLeftSide ;
3378  fDecayBits[5] = AliNeutralMesonSelection::kEtaBothSides;
3379  fDecayBits[6] = AliNeutralMesonSelection::kPi0LeftSide ; // Leave it last since likely not used
3380  fDecayBits[7] = AliNeutralMesonSelection::kPi0BothSides; // Leave it last since likely not used
3381 
3382  fNBkgBin = 11;
3383  fBkgBinLimit[ 0] = 00.0; fBkgBinLimit[ 1] = 00.2; fBkgBinLimit[ 2] = 00.3; fBkgBinLimit[ 3] = 00.4; fBkgBinLimit[ 4] = 00.5;
3384  fBkgBinLimit[ 5] = 01.0; fBkgBinLimit[ 6] = 01.5; fBkgBinLimit[ 7] = 02.0; fBkgBinLimit[ 8] = 03.0; fBkgBinLimit[ 9] = 05.0;
3385  fBkgBinLimit[10] = 10.0; fBkgBinLimit[11] = 100.;
3386  for(Int_t ibin = fNBkgBin+1; ibin < 20; ibin++)
3387  {
3388  fBkgBinLimit[ibin] = 00.0;
3389  }
3390 
3391 }
3392 
3393 //________________________________________________________________________________________________________
3398 //________________________________________________________________________________________________________
3399 void AliAnaParticleHadronCorrelation::InvMassHisto(AliAODPWG4ParticleCorrelation * trigger, Int_t mcIndex)
3400 {
3401  Int_t idTrig = trigger->GetCaloLabel(0);
3402  Float_t ptTrig = trigger->Pt();
3403  Float_t tofTrig = trigger->GetTime();
3404 
3405  fMomentum = *(trigger->Momentum());
3406 
3407  // Loop on second cluster to be associated to trigger particle to have the invariant mass
3408  Int_t nphoton = GetInputAODBranch()->GetEntriesFast();
3409 
3410  for(Int_t iphoton = 0; iphoton < nphoton; iphoton++)
3411  {
3412  AliAODPWG4ParticleCorrelation * photon1 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iphoton));
3413 
3414  if(idTrig == photon1->GetCaloLabel(0)) continue; //Do not the invariant mass for the same particle
3415 
3416  fMomentumIM = *(photon1->Momentum());
3417 
3418  // Select secondary photon with proper invariant mass
3419  if(fM02MaxCut > 0 && fM02MinCut > 0) //clID1 > 0 && clID2 < 0 &&
3420  {
3421  Float_t m02 = photon1->GetM02();
3422 
3423  if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3424  }
3425 
3426  // Select clusters with good time window difference
3427  Double_t tdiff = tofTrig - photon1->GetTime();
3428  if(TMath::Abs(tdiff) > GetPairTimeCut()) continue;
3429 
3430  // Add momenta of trigger and possible secondary decay
3431  fMomentumIM += fMomentum;
3432 
3433  Double_t mass = fMomentumIM.M();
3434  // Double_t epair = fMomentumIM.E();
3435  // Float_t ptpair = fMomentumIM.Pt();
3436 
3437  // Mass of all pairs
3438  fhMassPtTrigger->Fill(ptTrig, mass, GetEventWeight());
3439  if ( IsDataMC() ) fhMCMassPtTrigger[mcIndex]->Fill(ptTrig, mass, GetEventWeight());
3440  }
3441 }
3442 
3443 //_________________________________________________________________________
3446 //_________________________________________________________________________
3447 Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
3448 {
3449  Double_t ptTrig = GetMinPt();
3450  Double_t phiTrig = 0 ;
3451  fLeadingTriggerIndex =-1 ;
3452  Int_t index =-1 ;
3453  AliAODPWG4ParticleCorrelation* pLeading = 0;
3454 
3455  // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
3456 
3457  for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3458  {
3459  AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3460  particle->SetLeadingParticle(kFALSE); // set it later
3461 
3462  // Vertex cut in case of mixing
3463  Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3464  if(check == 0) continue;
3465  if(check == -1) return kFALSE; // not sure if it is correct.
3466 
3467  // Find the leading particles with highest momentum
3468  if (particle->Pt() > ptTrig)
3469  {
3470  ptTrig = particle->Pt() ;
3471  phiTrig = particle->Phi();
3472  index = iaod ;
3473  pLeading = particle ;
3474  }
3475  }// Finish search of leading trigger particle on the AOD branch.
3476 
3477  if(index < 0) return kFALSE;
3478 
3479  //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3480 
3481  if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3482 
3483  // Compare if it is the leading of all tracks
3484 
3485  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3486  {
3487  AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3488 
3489  // In case of isolation of single tracks or conversion photon (2 tracks) or pi0 (4 tracks),
3490  // do not count the candidate or the daughters of the candidate
3491  // in the isolation conte
3492  if ( pLeading->GetDetectorTag() == kCTS ) // make sure conversions are tagged as kCTS!!!
3493  {
3494  Int_t trackID = GetReader()->GetTrackID(track) ; // needed instead of track->GetID() since AOD needs some manipulations
3495  Bool_t contained = kFALSE;
3496 
3497  for(Int_t i = 0; i < 4; i++)
3498  {
3499  if( trackID == pLeading->GetTrackLabel(i) ) contained = kTRUE;
3500  }
3501 
3502  if ( contained ) continue ;
3503  }
3504 
3505  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3506  Float_t pt = fTrackVector.Pt();
3507  Float_t phi = fTrackVector.Phi() ;
3508  if(phi < 0) phi+=TMath::TwoPi();
3509 
3510  //jump out this event if near side associated particle pt larger than trigger
3511  if (fMakeNearSideLeading)
3512  {
3513  Float_t deltaPhi = phiTrig-phi;
3514  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3515  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3516 
3517  if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE;
3518  }
3519  //jump out this event if there is any other particle with pt larger than trigger
3520  else
3521  {
3522  if(pt > ptTrig) return kFALSE ;
3523  }
3524  }// track loop
3525 
3526  // Compare if it is leading of all calorimeter clusters
3527 
3528  if(fCheckLeadingWithNeutralClusters)
3529  {
3530  // Select the calorimeter cluster list
3531  TObjArray * nePl = 0x0;
3532  if (pLeading->GetDetectorTag() == kPHOS )
3533  nePl = GetPHOSClusters();
3534  else
3535  nePl = GetEMCALClusters();
3536 
3537  if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3538 
3539  for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3540  {
3541  AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3542 
3543  if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3544 
3545  cluster->GetMomentum(fMomentum,GetVertex(0));
3546 
3547  Float_t pt = fMomentum.Pt();
3548  Float_t phi = fMomentum.Phi() ;
3549  if(phi < 0) phi+=TMath::TwoPi();
3550 
3551  if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3552 
3553  //jump out this event if near side associated particle pt larger than trigger
3554  // not really needed for calorimeter, unless DCal is included
3555  if (fMakeNearSideLeading)
3556  {
3557  Float_t deltaPhi = phiTrig-phi;
3558  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3559  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3560 
3561  if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3562  }
3563  //jump out this event if there is any other particle with pt larger than trigger
3564  else
3565  {
3566  if(pt > ptTrig) return kFALSE ;
3567  }
3568  }// cluster loop
3569  } // check neutral clusters
3570 
3571  fLeadingTriggerIndex = index ;
3572  pLeading->SetLeadingParticle(kTRUE);
3573 
3574  AliDebug(1,Form("\t particle AOD with index %d is leading with pT %2.2f", fLeadingTriggerIndex, pLeading->Pt()));
3575 
3576  return kTRUE;
3577 }
3578 
3579 //_________________________________________________________________
3593 //_________________________________________________________________
3594 void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
3595 {
3596  if(!GetInputAODBranch())
3597  {
3598  AliFatal(Form("No input particles in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3599  return ; // coverity
3600  }
3601 
3602  Int_t naod = GetInputAODBranch()->GetEntriesFast();
3603  if( naod == 0 )
3604  {
3605  AliDebug(1,"No particle AOD found!");
3606  return ; // no trigger particles found.
3607  }
3608 
3609  AliDebug(1,Form("Begin hadron correlation analysis, fill histograms"));
3610  AliDebug(1,Form("n particle branch aod entries %d", naod));
3611  AliDebug(1,Form("In CTS aod entries %d",GetCTSTracks()->GetEntriesFast()));
3612 
3613  //------------------------------------------------------
3614  // Find leading trigger if analysis request only leading,
3615  // if there is no leading trigger, then skip the event
3616 
3617  Int_t iaod = 0 ;
3618  if( fMakeAbsoluteLeading || fMakeNearSideLeading )
3619  {
3620  Bool_t leading = IsTriggerTheEventLeadingParticle();
3621 
3622  AliDebug(1,Form("AOD Leading trigger? %d, with index %d",leading,fLeadingTriggerIndex));
3623 
3624  if(!leading)
3625  {
3626  AliDebug(1,"Leading was requested and not found");
3627  return ;
3628  }
3629  else
3630  {
3631  // Select only the leading in the trigger AOD loop
3632  naod = fLeadingTriggerIndex+1 ;
3633  iaod = fLeadingTriggerIndex ;
3634  }
3635  }
3636 
3637  //------------------------------------------------------
3638  // Get event multiplicity and bins
3639 
3640  Float_t cen = GetEventCentrality();
3641  Float_t ep = GetEventPlaneAngle();
3642  if(IsHighMultiplicityAnalysisOn()) fhTriggerEventPlaneCentrality->Fill(cen, ep, GetEventWeight());
3643 
3644  Int_t mixEventBin = GetEventMixBin();
3645  Int_t vzbin = GetEventVzBin();
3646 
3647  //------------------------------------------------------
3648  // Loop on trigger AOD
3649 
3650  for( iaod = 0; iaod < naod; iaod++ )
3651  {
3652  AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3653 
3654  //
3655  // Trigger particle selection criteria:
3656  //
3657  Float_t pt = particle->Pt();
3658 
3659  if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3660 
3661  fhPtTriggerInput->Fill(pt, GetEventWeight());
3662 
3663  //
3664  // check if it was a calorimeter cluster
3665  // and if the shower shape cut was requested apply it.
3666  // Not needed if already done at the particle identification level,
3667  // but for isolation studies, it is preferred not to remove so we do it here
3668  //
3669 
3670  AliDebug(1,Form("%s Trigger : min %f, max %f, det %d",
3671  GetInputAODName().Data(),fM02MinCut,fM02MaxCut,particle->GetDetectorTag()));
3672 
3673  if(fM02MaxCut > 0 && fM02MinCut > 0) //clID1 > 0 && clID2 < 0 &&
3674  {
3675  // Int_t iclus = -1;
3676  // TObjArray* clusters = 0x0;
3677  // if (particle->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
3678  // else if(particle->GetDetectorTag() == kPHOS ) clusters = GetPHOSClusters();
3679  //
3680  // if(clusters)
3681  // {
3682  // AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
3683  // Float_t m02 = cluster->GetM02();
3684  // if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3685  // }
3686 
3687  Float_t m02 = particle->GetM02();
3688 
3689  if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3690 
3691  fhPtTriggerSSCut->Fill(pt, GetEventWeight());
3692 
3693  AliDebug(1,"Pass the shower shape cut");
3694  }
3695 
3696  //
3697  // Check if the particle is isolated or if we want to take the isolation into account
3698  // This bool is set in AliAnaParticleIsolation
3699  //
3700  if(OnlyIsolated())
3701  {
3702  if( !particle->IsIsolated() ) continue;
3703 
3704  fhPtTriggerIsoCut->Fill(pt, GetEventWeight());
3705 
3706  AliDebug(1,"Pass the isolation cut");
3707  }
3708 
3709  //
3710  // Check if trigger is in fiducial region
3711  //
3712  if(IsFiducialCutOn())
3713  {
3714  Bool_t in = GetFiducialCut()->IsInFiducialCut(particle->Eta(),particle->Phi(),particle->GetDetectorTag()) ;
3715 
3716  if(! in ) continue ;
3717 
3718  AliDebug(1,"Pass the fiducial cut");
3719  }
3720 
3721  fhPtTriggerFidCut->Fill(pt, GetEventWeight());
3722 
3723  //---------------------------------------
3724  // Make correlation
3725 
3726  // Find the leading hadron in the opposite hemisphere to the trigger
3727  // and accept the trigger if leading is in defined window.
3728  Bool_t okLeadHad = kTRUE;
3729  if(fSelectLeadingHadronAngle || fFillLeadHadOppositeHisto)
3730  {
3731  okLeadHad = FindLeadingOppositeHadronInWindow(particle);
3732  if(!okLeadHad && fSelectLeadingHadronAngle) continue;
3733  }
3734 
3735  //
3736  // Charged particles correlation
3737  //
3738  MakeChargedCorrelation(particle);
3739 
3740  // MC
3741  Int_t mcIndex = -1;
3742  Int_t mcTag = particle->GetTag();
3743  Bool_t lostDecayPair = kFALSE;
3744  if(IsDataMC())
3745  {
3746  mcIndex = GetMCTagHistogramIndex(mcTag);
3747  lostDecayPair = GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost);
3748  MakeMCChargedCorrelation(particle->GetLabel(), mcIndex,lostDecayPair);
3749  }
3750 
3751  // Do own mixed event with charged,
3752  // add event and remove previous or fill the mixed histograms
3753  if(DoOwnMix())
3754  MakeChargedMixCorrelation(particle);
3755 
3756  //
3757  // Neutral particles correlation
3758  //
3759  if(fNeutralCorr)
3760  MakeNeutralCorrelation(particle);
3761 
3762  //----------------------------------------------------------------
3763  // Fill trigger pT related histograms if not absolute leading
3764 
3765  //
3766  // pT of the trigger, vs trigger origin if MC
3767  //
3768  fhPtTrigger->Fill(pt, GetEventWeight());
3769 
3770 
3771  if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3772  {
3773  fhPtTriggerMC[mcIndex]->Fill(pt, GetEventWeight());
3774  if( lostDecayPair )
3775  {
3776  // check index in GetMCTagIndexHistogram
3777  if ( mcIndex == 2 ) fhPtTriggerMC[8]->Fill(pt, GetEventWeight()); // pi0 decay
3778  else if( mcIndex == 4 ) fhPtTriggerMC[9]->Fill(pt, GetEventWeight()); // eta decay
3779  }
3780  }
3781 
3782  //
3783  // pT lead cone, pT sum cone background bins distribution
3784  //
3785  if(fFillBkgBinsHisto)
3786  {
3787  Float_t m02 = particle->GetM02();
3788  Float_t pTLeadTrackInCone = 0;
3789  Float_t pTSumTrackInCone = 0;
3790  Float_t pTLeadClusterInCone = 0;
3791  Float_t pTSumClusterInCone = 0;
3792 
3793  pTLeadTrackInCone = particle->GetChargedLeadPtInCone();
3794  pTLeadClusterInCone = particle->GetNeutralLeadPtInCone();
3795 
3796  pTSumTrackInCone = particle->GetChargedPtSumInCone();
3797  pTSumClusterInCone = particle->GetNeutralPtSumInCone();
3798 
3799  Float_t pTLeadInCone = pTLeadTrackInCone;
3800  if(pTLeadClusterInCone > pTLeadInCone) pTLeadInCone = pTLeadClusterInCone;
3801  Float_t pTSumInCone = pTSumTrackInCone + pTSumClusterInCone;
3802 
3803  // Get the background bin for this cone and trigger
3804  Int_t pTSumBin = -1;
3805  Int_t pTLeadBin = -1;
3806 
3807  for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
3808  {
3809  if( pTSumInCone >= fBkgBinLimit[ibin] && pTSumInCone < fBkgBinLimit[ibin+1]) pTSumBin = ibin;
3810  if( pTLeadInCone >= fBkgBinLimit[ibin] && pTLeadInCone < fBkgBinLimit[ibin+1]) pTLeadBin = ibin;
3811  }
3812  if(pTSumBin > 0)
3813  {
3814  fhPtSumInConeBin[pTSumBin]->Fill(pt, GetEventWeight());
3815  }
3816 
3817  if(pTLeadBin > 0)
3818  {
3819  fhPtLeadInConeBin[pTLeadBin]->Fill(pt, GetEventWeight());
3820  }
3821 
3822  // Check if it was a decay
3823  if( fFillTaggedDecayHistograms && m02 < fDecayTagsM02Cut )
3824  {
3825  Int_t decayTag = particle->DecayTag();
3826  if(decayTag < 0) decayTag = 0;
3827 
3828  for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
3829  {
3830  if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3831  {
3832  Int_t pTLeadBinDecay = pTLeadBin+ibit*fNBkgBin;
3833  Int_t pTSumBinDecay = pTSumBin+ibit*fNBkgBin;
3834  if( pTLeadBin >=0 ) fhPtLeadConeBinDecay[pTLeadBinDecay]->Fill(pt, GetEventWeight());
3835  if( pTSumBin >=0 ) fhSumPtConeBinDecay [pTSumBinDecay] ->Fill(pt, GetEventWeight());
3836  }
3837  }
3838  }
3839 
3840  if(IsDataMC())
3841  {
3842  Int_t pTLeadBinMC = pTLeadBin+mcIndex*fNBkgBin;
3843  Int_t pTSumBinMC = pTSumBin+mcIndex*fNBkgBin;
3844 
3845  if( pTLeadBin >=0 )
3846  {
3847  fhPtLeadConeBinMC[pTLeadBinMC]->Fill(pt, GetEventWeight());
3848  }
3849 
3850  if( pTSumBin >=0 )
3851  {
3852  fhSumPtConeBinMC [pTSumBinMC]->Fill(pt, GetEventWeight());
3853  }
3854 
3855  if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3856  {
3857  pTLeadBinMC = pTLeadBin+kmcPhoton*fNBkgBin;
3858  pTSumBinMC = pTSumBin+kmcPhoton*fNBkgBin;
3859  if( pTLeadBin >=0 )
3860  {
3861  fhPtLeadConeBinMC[pTLeadBinMC]->Fill(pt, GetEventWeight());
3862  }
3863 
3864  if( pTSumBin >=0 )
3865  {
3866  fhSumPtConeBinMC [ pTSumBinMC]->Fill(pt, GetEventWeight());
3867  }
3868  }
3869 
3870  // Check if decay and if pair is lost
3871  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
3872  {
3873  if ( mcIndex == kmcPi0Decay )
3874  {
3875  pTLeadBinMC = pTLeadBin+kmcPi0DecayLostPair*fNBkgBin;
3876  pTSumBinMC = pTSumBin+kmcPi0DecayLostPair*fNBkgBin;
3877  }
3878  else if(mcIndex == kmcEtaDecay)
3879  {
3880  pTLeadBinMC = pTLeadBin+kmcEtaDecayLostPair*fNBkgBin;
3881  pTSumBinMC = pTSumBin+kmcEtaDecayLostPair*fNBkgBin;
3882  }
3883  else
3884  AliFatal(Form("Lost decay Bit assigned to bad case, mcIndex %d",mcIndex));
3885 
3886  if( pTLeadBin >=0 )
3887  {
3888  fhPtLeadConeBinMC[pTLeadBinMC]->Fill(pt, GetEventWeight());
3889  }
3890 
3891  if( pTSumBin >=0 )
3892  {
3893  fhSumPtConeBinMC [ pTSumBinMC]->Fill(pt, GetEventWeight());
3894  }
3895 
3896  }
3897 
3898  }
3899  }
3900 
3901  // Invariant mass of trigger particle
3902  if(fFillInvMassHisto) InvMassHisto(particle,mcIndex);
3903 
3904  if(fDecayTrigger)
3905  {
3906  Int_t decayTag = particle->DecayTag();
3907  if(decayTag < 0) decayTag = 0;
3908 
3909  for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
3910  {
3911  if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3912  {
3913  fhPtDecayTrigger[ibit]->Fill(pt, GetEventWeight());
3914 
3915  if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3916  {
3917  fhPtDecayTriggerMC[ibit][mcIndex]->Fill(pt, GetEventWeight());
3918  if( lostDecayPair )
3919  {
3920  // check index in GetMCTagIndexHistogram
3921  if( mcIndex == 2 ) fhPtDecayTriggerMC[ibit][8]->Fill(pt, GetEventWeight()); // pi0 decay
3922  if( mcIndex == 4 ) fhPtDecayTriggerMC[ibit][9]->Fill(pt, GetEventWeight()); // eta decay
3923  }
3924  }
3925  }// check bit
3926  }// bit loop
3927  }
3928 
3929  //
3930  // Acceptance of the trigger
3931  //
3932  Float_t phi = particle->Phi();
3933  if( phi < 0 ) phi+=TMath::TwoPi();
3934  fhPhiTrigger->Fill(pt, phi, GetEventWeight());
3935 
3936  fhEtaTrigger->Fill(pt, particle->Eta(), GetEventWeight());
3937  //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Trigger particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
3938 
3939  //----------------------------------
3940  // Trigger particle pT vs event bins
3941 
3942  fhPtTriggerBin->Fill(pt, mixEventBin, GetEventWeight());
3943  if(fCorrelVzBin)
3944  fhPtTriggerVzBin->Fill(pt, vzbin, GetEventWeight());
3945 
3946  if(IsHighMultiplicityAnalysisOn())
3947  {
3948  fhPtTriggerCentrality->Fill(pt, cen, GetEventWeight());
3949  fhPtTriggerEventPlane->Fill(pt, ep , GetEventWeight());
3950  }
3951 
3952  //----------------------------------
3953  // Trigger particle pT vs pile-up
3954 
3955  if(IsPileUpAnalysisOn())
3956  {
3957  Int_t vtxBC = GetReader()->GetVertexBC();
3958  if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA) fhPtTriggerVtxBC0->Fill(pt, GetEventWeight());
3959 
3960  if(GetReader()->IsPileUpFromSPD()) fhPtTriggerPileUp[0]->Fill(pt, GetEventWeight());
3961  if(GetReader()->IsPileUpFromEMCal()) fhPtTriggerPileUp[1]->Fill(pt, GetEventWeight());
3962  if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtTriggerPileUp[2]->Fill(pt, GetEventWeight());
3963  if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtTriggerPileUp[3]->Fill(pt, GetEventWeight());
3964  if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtTriggerPileUp[4]->Fill(pt, GetEventWeight());
3965  if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtTriggerPileUp[5]->Fill(pt, GetEventWeight());
3966  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTriggerPileUp[6]->Fill(pt, GetEventWeight());
3967  }
3968  } // AOD trigger loop
3969 
3970  // Re-init for next event
3971  fLeadingTriggerIndex = -1;
3972 
3973  AliDebug(1,"End fill histograms");
3974 }
3975 
3976 //_______________________________________________________________________________________________________
3982 //_______________________________________________________________________________________________________
3983 void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3984 {
3985  AliDebug(1,"Make trigger particle - charged hadron correlation");
3986 
3987  Float_t phiTrig = aodParticle->Phi();
3988  Float_t etaTrig = aodParticle->Eta();
3989  Float_t ptTrig = aodParticle->Pt();
3990  Int_t mcTag = aodParticle->GetTag();
3991  Double_t bz = GetReader()->GetInputEvent()->GetMagneticField();
3992 
3993 
3994  Int_t decayTag = 0;
3995  if(fDecayTrigger)
3996  {
3997  //decay = aodParticle->IsTagged();
3998  decayTag = aodParticle->DecayTag();
3999  if(decayTag < 0) decayTag = 0;
4000  // printf("Correlation: pT %2.2f, BTag %d, Tagged %d\n",ptTrig, decayTag, aodParticle->IsTagged());
4001  // printf("\t check bit Pi0 %d, Eta %d, Pi0Side %d, EtaSide %d\n",
4002  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0),
4003  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEta),
4004  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0Side),
4005  // GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEtaSide));
4006  }
4007 
4008  Float_t pt = -100. ;
4009  Float_t phi = -100. ;
4010  Float_t eta = -100. ;
4011  Float_t deltaPhi = -100. ;
4012 
4013  TObjArray * reftracks = 0x0;
4014  Int_t nrefs = 0;
4015 
4016  // Mixed event settings
4017  Int_t evtIndex11 = -1 ; // cluster trigger or pi0 trigger
4018  Int_t evtIndex12 = -1 ; // pi0 trigger
4019  Int_t evtIndex13 = -1 ; // charged trigger
4020 
4021  if (GetMixedEvent())
4022  {
4023  evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
4024  evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
4025  evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
4026  }
4027 
4028  // Track multiplicity or cent bin
4029  Int_t cenbin = 0;
4030  if(IsHighMultiplicityAnalysisOn()) cenbin = GetEventCentralityBin();
4031 
4032  //
4033  // In case of pi0/eta trigger, we may want to check their decay correlation,
4034  // get their decay children
4035  //
4036 
4037  Bool_t decayFound = kFALSE;
4038  if( fPi0Trigger )
4039  {
4040  decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
4041  if(decayFound)
4042  {
4043  fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom1.Pt()/ptTrig, GetEventWeight());
4044  fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom2.Pt()/ptTrig, GetEventWeight());
4045  }
4046  }
4047 
4048  //-----------------------------------------------------------------------
4049  // Track loop, select tracks with good pt, phi and fill AODs or histograms
4050  //-----------------------------------------------------------------------
4051 
4052  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
4053  {
4054  AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
4055 
4056  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
4057  pt = fTrackVector.Pt();
4058  eta = fTrackVector.Eta();
4059  phi = fTrackVector.Phi() ;
4060  if(phi < 0) phi+=TMath::TwoPi();
4061 
4062  //Select only hadrons in pt range
4063  if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
4064 
4065  // In case of isolation of single tracks or conversion photon (2 tracks) or pi0 (4 tracks),
4066  // do not count the candidate or the daughters of the candidate
4067  // in the isolation conte
4068  if ( aodParticle->GetDetectorTag() == kCTS ) // make sure conversions are tagged as kCTS!!!
4069  {
4070  Int_t trackID = GetReader()->GetTrackID(track) ; // needed instead of track->GetID() since AOD needs some manipulations
4071  Bool_t contained = kFALSE;
4072 
4073  for(Int_t i = 0; i < 4; i++)
4074  {
4075  if( trackID == aodParticle->GetTrackLabel(i) ) contained = kTRUE;
4076  }
4077 
4078  if ( contained ) continue ;
4079  }
4080 
4081  //Only for mixed event frame
4082  Int_t evtIndex2 = 0 ;
4083  if (GetMixedEvent())
4084  {
4085  evtIndex2 = GetMixedEvent()->EventIndex(GetReader()->GetTrackID(track)) ;
4086  if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
4087  continue ;
4088  //vertex cut
4089  if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
4090  continue;
4091  }
4092 
4093  AliDebug(2,Form("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
4094 
4095  // ------------------------------
4096  // Track type bin or bits setting
4097  //
4098 
4099  //
4100  // * Set the pt associated bin for the defined bins *
4101  //
4102  Int_t assocBin = -1;
4103  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
4104  {
4105  if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
4106  }
4107 
4108  //
4109  // * Assign to the histogram array a bin corresponding
4110  // to a combination of pTa and vz bins *
4111  //
4112  Int_t nz = 1;
4113  Int_t vz = 0;
4114 
4115  if(fCorrelVzBin)
4116  {
4117  nz = GetNZvertBin();
4118  vz = GetEventVzBin();
4119  }
4120 
4121  Int_t bin = assocBin*nz+vz;
4122 
4123  //printf("assoc Bin = %d, vZ bin = %d, bin = %d \n", assocBin,GetEventVzBin(),bin);
4124 
4125  //
4126  // * Get the status of the TOF bit *
4127  //
4128  ULong_t status = track->GetStatus();
4129  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
4130  //Double32_t tof = track->GetTOFsignal()*1e-3;
4131  Int_t trackBC = track->GetTOFBunchCrossing(bz);
4132 
4133  Int_t outTOF = -1;
4134  if (okTOF && trackBC!=0) outTOF = 1;
4135  else if(okTOF && trackBC==0) outTOF = 0;
4136 
4137  //----------------
4138  // Fill Histograms
4139 
4140  //
4141  // Azimuthal Angle histograms
4142  //
4143 
4144  deltaPhi = phiTrig-phi;
4145 
4146  //
4147  // Calculate deltaPhi shift so that for the particles on the opposite side
4148  // it is defined between 90 and 270 degrees
4149  // Shift [-360,-90] to [0, 270]
4150  // and [270,360] to [-90,0]
4151  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4152  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4153 
4154  FillChargedAngularCorrelationHistograms(pt, ptTrig, bin, phi, phiTrig, deltaPhi,
4155  eta, etaTrig, decayTag, track->GetHMPIDsignal(),
4156  outTOF, cenbin, mcTag);
4157 
4158  //
4159  // Imbalance zT/xE/pOut histograms
4160  //
4161 
4162  AliESDtrack * esdTrack = dynamic_cast<AliESDtrack*>(track);
4163 
4164  //
4165  // Delta phi cut for momentum imbalance correlation
4166  //
4167  if ( (deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) )
4168  {
4169  // Check track resolution
4170  if ( esdTrack )
4171  {
4172  fhTrackResolution->Fill(pt, TMath::Sqrt(esdTrack->GetCovariance()[14])*pt,
4173  GetEventWeight());
4174  }
4175 
4176  FillChargedMomentumImbalanceHistograms(ptTrig, pt, deltaPhi, cenbin, track->Charge(),
4177  assocBin, decayTag, outTOF, mcTag);
4178  }
4179 
4180  //
4181  // Underlying event, right side, default case
4182  //
4183  if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4184  {
4185  // Check track resolution
4186  if ( esdTrack )
4187  {
4188  fhTrackResolutionUE->Fill(pt, TMath::Sqrt(esdTrack->GetCovariance()[14])*pt,
4189  GetEventWeight());
4190  }
4191 
4192  FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, cenbin, outTOF,mcTag);
4193  }
4194 
4195  //
4196  // Several UE calculation, in different perpendicular regions, up to 6:
4197  // left, right, upper-left, lower left, upper-right, lower-right
4198  //
4199  if(fMakeSeveralUE)
4200  FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,deltaPhi,mcTag);
4201 
4202  //
4203  if(fPi0Trigger && decayFound)
4204  FillDecayPhotonCorrelationHistograms(pt, phi, kTRUE) ;
4205 
4206  //
4207  // Add track reference to array
4208  //
4209  if(fFillAODWithReferences)
4210  {
4211  nrefs++;
4212  if(nrefs==1)
4213  {
4214  reftracks = new TObjArray(0);
4215  TString trackname = Form("%sTracks", GetAODObjArrayName().Data());
4216  reftracks->SetName(trackname.Data());
4217  reftracks->SetOwner(kFALSE);
4218  }
4219 
4220  reftracks->Add(track);
4221  }// reference track to AOD
4222  }// track loop
4223 
4224  //Fill AOD with reference tracks, if not filling histograms
4225  if(fFillAODWithReferences && reftracks)
4226  {
4227  aodParticle->AddObjArray(reftracks);
4228  }
4229 }
4230 
4231 //_________________________________________________________________________________________________________
4233 //_________________________________________________________________________________________________________
4234 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
4235 {
4236  AliDebug(1,Form("Make trigger particle - charged hadron mixed event correlation"));
4237 
4238  if(GetMixedEvent()) return; // This is not the mixed event from general mixing frame
4239 
4240  // Get the event with similar caracteristics
4241  //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
4242 
4243  AliAnalysisManager * manager = AliAnalysisManager::GetAnalysisManager();
4244 
4245  AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
4246 
4247  if(!inputHandler) return;
4248 
4249  if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
4250 
4251  // Get the pool, check if it exits
4252  Int_t eventBin = GetEventMixBin();
4253 
4254  //Check that the bin exists, if not (bad determination of RP, ntrality or vz bin) do nothing
4255  if(eventBin < 0) return;
4256 
4257  fhEventBin->Fill(eventBin, GetEventWeight());
4258 
4259  // Get neutral clusters pool?
4260  Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
4261  Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
4262 
4263  TList * pool = 0;
4264  TList * poolCalo = 0;
4265  if(fUseMixStoredInReader)
4266  {
4267  pool = GetReader()->GetListWithMixedEventsForTracks(eventBin);
4268  if(neutralMix) poolCalo = GetReader()->GetListWithMixedEventsForCalo (eventBin);
4269  }
4270  else
4271  {
4272  pool = fListMixTrackEvents[eventBin];
4273  if(neutralMix) poolCalo = fListMixCaloEvents [eventBin];
4274  }
4275 
4276  if(!pool) return ;
4277 
4278  if( neutralMix && !poolCalo )
4279  AliWarning("Careful, cluster pool not available");
4280 
4281  Double_t ptTrig = aodParticle->Pt();
4282  Double_t etaTrig = aodParticle->Eta();
4283  Double_t phiTrig = aodParticle->Phi();
4284  if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
4285 
4286  AliDebug(1,Form("Pool bin %d size %d, trigger trigger pt=%f, phi=%f, eta=%f",
4287  eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig));
4288 
4289  Double_t ptAssoc = -999.;
4290  Double_t phiAssoc = -999.;
4291  Double_t etaAssoc = -999.;
4292  Double_t deltaPhi = -999.;
4293  Double_t deltaEta = -999.;
4294  Double_t xE = -999.;
4295 
4296  // Start from first event in pool except if in this same event the pool was filled
4297  Int_t ev0 = 0;
4298  if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
4299 
4300  for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
4301  {
4302  //
4303  // Recover the lists of tracks or clusters
4304  //
4305  TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
4306  TObjArray* bgCalo = 0;
4307 
4308  // Recover the clusters list if requested
4309  if( neutralMix && poolCalo )
4310  {
4311  if(pool->GetSize()!=poolCalo->GetSize())
4312  AliWarning("Different size of calo and track pools");
4313 
4314  bgCalo = static_cast<TObjArray*>(poolCalo->At(ev));
4315 
4316  if(!bgCalo) AliDebug(1,Form("Event %d in calo pool not available?",ev));
4317  }
4318 
4319  //
4320  // Isolate the trigger in the mixed event with mixed tracks and clusters
4321  //
4322  if( OnlyIsolated() )
4323  {
4324  Int_t n=0, nfrac = 0;
4325  Bool_t isolated = kFALSE;
4326  Float_t coneptsum = 0, coneptlead = 0;
4327  GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
4328  GetReader(), GetCaloPID(),
4329  kFALSE, aodParticle, "",
4330  n,nfrac,coneptsum,coneptlead,isolated);
4331 
4332  //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
4333  // isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
4334  //if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
4335  //printf("\n");
4336 
4337  if(!isolated) continue ;
4338  }
4339 
4340  //
4341  // Check if the trigger is leading of mixed event
4342  //
4343  Int_t nTracks=bgTracks->GetEntriesFast();
4344 
4345  if(fMakeNearSideLeading || fMakeAbsoluteLeading)
4346  {
4347  Bool_t leading = kTRUE;
4348  for(Int_t jlead = 0;jlead < nTracks; jlead++ )
4349  {
4350  AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(jlead) ;
4351 
4352  ptAssoc = track->Pt();
4353  phiAssoc = track->Phi() ;
4354  if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
4355 
4356  if (fMakeNearSideLeading)
4357  {
4358  deltaPhi = phiTrig-phiAssoc;
4359  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4360  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4361 
4362  if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
4363  {
4364  leading = kFALSE;
4365  break;
4366  }
4367  }
4368  //jump out this event if there is any other particle with pt larger than trigger
4369  else if(fMakeAbsoluteLeading)
4370  {
4371  if(ptAssoc > ptTrig)
4372  {
4373  leading = kFALSE;
4374  break;
4375  }
4376  }
4377  }
4378 
4379  if( !neutralMix && fCheckLeadingWithNeutralClusters )
4380  AliWarning("Leading of clusters requested but no clusters in mixed event");
4381 
4382  if(neutralMix && fCheckLeadingWithNeutralClusters && bgCalo)
4383  {
4384  Int_t nClusters=bgCalo->GetEntriesFast();
4385  for(Int_t jlead = 0;jlead <nClusters; jlead++ )
4386  {
4387  AliAODPWG4Particle *cluster= (AliAODPWG4Particle*) bgCalo->At(jlead) ;
4388 
4389  ptAssoc = cluster->Pt();
4390  phiAssoc = cluster->Phi() ;
4391  if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
4392 
4393  if (fMakeNearSideLeading)
4394  {
4395  deltaPhi = phiTrig-phiAssoc;
4396  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4397  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4398 
4399  if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
4400  {
4401  leading = kFALSE;
4402  break;
4403  }
4404  }
4405  //jump out this event if there is any other particle with pt larger than trigger
4406  else if(fMakeAbsoluteLeading)
4407  {
4408  if(ptAssoc > ptTrig)
4409  {
4410  leading = kFALSE;
4411  break;
4412  }
4413  }
4414  }
4415  }
4416 
4417  if(!leading) continue; // not leading, check the next event in pool
4418  }
4419 
4420  //
4421  // Fill histograms for selected triggers
4422  //
4423 
4424  fhEventMixBin->Fill(eventBin, GetEventWeight());
4425 
4426  //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
4427 
4428  fhPtTriggerMixed ->Fill(ptTrig, GetEventWeight());
4429  fhPhiTriggerMixed ->Fill(ptTrig, phiTrig, GetEventWeight());
4430  fhEtaTriggerMixed ->Fill(ptTrig, etaTrig, GetEventWeight());
4431  fhPtTriggerMixedBin->Fill(ptTrig, eventBin, GetEventWeight());
4432 
4433  if ( fCorrelVzBin ) fhPtTriggerMixedVzBin->Fill(ptTrig, GetEventVzBin(), GetEventWeight());
4434 
4435  //
4436  // Correlation histograms
4437  //
4438  for(Int_t j1 = 0;j1 <nTracks; j1++ )
4439  {
4440  AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
4441 
4442  if(!track) continue;
4443 
4444  ptAssoc = track->Pt();
4445  etaAssoc = track->Eta();
4446  phiAssoc = track->Phi() ;
4447  if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
4448 
4449  deltaPhi = phiTrig-phiAssoc;
4450  if(deltaPhi < -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4451  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4452  deltaEta = etaTrig-etaAssoc;
4453 
4454  AliDebug(1,Form("deltaPhi= %f, deltaEta=%f",deltaPhi, deltaEta));
4455 
4456  // Angular correlation
4457  fhMixDeltaPhiCharged ->Fill(ptTrig , deltaPhi, GetEventWeight());
4458  fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta, GetEventWeight());
4459 
4460  //
4461  // Momentum imbalance
4462  //
4463  if ( (deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) )
4464  {
4465  xE = -ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
4466 
4467  if( xE < 0.)
4468  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4469  xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4470 
4471  fhMixXECharged->Fill(ptTrig, xE, GetEventWeight());
4472  if ( xE > 0 ) fhMixHbpXECharged->Fill(ptTrig, TMath::Log(1./xE), GetEventWeight());
4473  }
4474 
4475  //
4476  // Underlying event momentum imbalance
4477  //
4478  if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4479  {
4480  //Underlying event region
4481  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
4482  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
4483 
4484  if(uexE < 0.)
4485  AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4486  uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4487 
4488  fhMixXEUeCharged->Fill(ptTrig, uexE, GetEventWeight());
4489  }
4490 
4491  // Set the pt associated bin for the defined bins
4492  Int_t assocBin = -1;
4493  for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
4494  {
4495  if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
4496  }
4497 
4498  //
4499  // Assign to the histogram array a bin corresponding to a combination of pTa and vz bins
4500  //
4501  Int_t nz = 1;
4502  Int_t vz = 0;
4503 
4504  if(fCorrelVzBin)
4505  {
4506  nz = GetNZvertBin();
4507  vz = GetEventVzBin();
4508  }
4509 
4510  Int_t bin = assocBin*nz+vz;
4511 
4512  if(bin < 0) continue ; // this pt bin was not considered
4513 
4514  fhMixDeltaPhiChargedAssocPtBin [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
4515  fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->Fill(deltaPhi, deltaEta, GetEventWeight());
4516 
4517  if(fFillEtaGapsHisto)
4518  {
4519  if(TMath::Abs(deltaEta) > 0.8)
4520  fhMixDeltaPhiChargedAssocPtBinDEta08 [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
4521  if(TMath::Abs(deltaEta) < 0.01)
4522  fhMixDeltaPhiChargedAssocPtBinDEta0 [bin]->Fill(ptTrig, deltaPhi, GetEventWeight());
4523  }
4524 
4525  } // track loop
4526  } // mixed event loop
4527 }
4528 
4529 //_______________________________________________________________________________________________________
4532 //_______________________________________________________________________________________________________
4533 void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle)
4534 {
4535  TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); // For the future, foresee more possible pi0 lists
4536  if(!pi0list) return ;
4537 
4538  Int_t npi0 = pi0list->GetEntriesFast();
4539  if(npi0 == 0) return ;
4540 
4541  AliDebug(1,Form("Particle - pi0 correlation, %d pi0's",npi0));
4542 
4543  Int_t evtIndex11 = 0 ;
4544  Int_t evtIndex12 = 0 ;
4545  if (GetMixedEvent())
4546  {
4547  evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
4548  evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
4549  }
4550 
4551  Float_t pt = -100. ;
4552  Float_t zT = -100. ;
4553  Float_t phi = -100. ;
4554  Float_t eta = -100. ;
4555  Float_t xE = -100. ;
4556  Float_t hbpXE= -100. ;
4557  Float_t hbpZT= -100. ;
4558 
4559  Float_t ptTrig = aodParticle->Pt();
4560  Float_t phiTrig = aodParticle->Phi();
4561  Float_t etaTrig = aodParticle->Eta();
4562  Float_t deltaPhi= -100. ;
4563  Float_t deltaEta= -100. ;
4564 
4565  // In case of pi0/eta trigger, we may want to check their decay correlation,
4566  // get their decay children
4567 
4568  Bool_t decayFound = kFALSE;
4569  if(fPi0Trigger) decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
4570 
4571  TObjArray * refpi0 = 0x0;
4572  Int_t nrefs = 0;
4573 
4574  // Loop on stored AOD pi0
4575 
4576  for(Int_t iaod = 0; iaod < npi0 ; iaod++)
4577  {
4578  AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (pi0list->At(iaod));
4579 
4580  Int_t evtIndex2 = 0 ;
4581  Int_t evtIndex3 = 0 ;
4582  if (GetMixedEvent())
4583  {
4584  evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
4585  evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
4586 
4587  if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
4588  evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
4589  continue ;
4590  }
4591 
4592  pt = pi0->Pt();
4593 
4594  if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
4595 
4596  // Remove trigger itself for correlation when use charged triggers
4597  if(aodParticle->GetCaloLabel(0) >= 0 &&
4598  (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(0) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(0))) continue ;
4599 
4600  if( aodParticle->GetCaloLabel(1) >= 0 &&
4601  (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(1) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(1))) continue ;
4602 
4603  //
4604  // Angular correlations
4605  //
4606  phi = pi0->Phi() ;
4607  eta = pi0->Eta() ;
4608  deltaEta = etaTrig-eta;
4609  deltaPhi = phiTrig-phi;
4610  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4611  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4612 
4613  fhEtaNeutral ->Fill(pt , eta , GetEventWeight());
4614  fhPhiNeutral ->Fill(pt , phi , GetEventWeight());
4615  fhDeltaEtaNeutral->Fill(ptTrig, deltaEta, GetEventWeight());
4616  fhDeltaPhiNeutral->Fill(ptTrig, deltaPhi, GetEventWeight());
4617 
4618  if ( pt > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta, GetEventWeight());
4619 
4620  //
4621  // Momentum imbalance
4622  //
4623  zT = pt/ptTrig ;
4624 
4625  hbpZT = -100;
4626  hbpXE = -100;
4627 
4628  if(zT > 0 ) hbpZT = TMath::Log(1./zT);
4629 
4630  // Delta phi cut for correlation
4631  if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
4632  {
4633  xE =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
4634 
4635  if( xE < 0. )
4636  AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4637  xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4638 
4639  if( xE > 0 ) hbpXE = TMath::Log(1./xE);
4640 
4641  fhDeltaPhiNeutralPt->Fill(pt , deltaPhi, GetEventWeight());
4642  fhXENeutral ->Fill(ptTrig, xE , GetEventWeight());
4643  fhPtHbpXENeutral ->Fill(ptTrig, hbpXE , GetEventWeight());
4644  fhZTNeutral ->Fill(ptTrig, zT , GetEventWeight());
4645  fhPtHbpZTNeutral ->Fill(ptTrig, hbpZT , GetEventWeight());
4646  }
4647  else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4648  {
4649  // Randomize angle for xE calculation
4650  Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
4651 
4652  xE = -(pt/ptTrig)*TMath::Cos(randomphi);
4653  if(xE > 0 ) hbpXE = TMath::Log(1./xE);
4654 
4655  fhDeltaPhiUeNeutralPt->Fill(pt , deltaPhi, GetEventWeight());
4656  fhZTUeNeutral ->Fill(ptTrig, zT , GetEventWeight());
4657  fhPtHbpZTUeNeutral ->Fill(ptTrig, hbpZT , GetEventWeight());
4658  fhXEUeNeutral ->Fill(ptTrig, xE , GetEventWeight());
4659  fhPtHbpXEUeNeutral ->Fill(ptTrig, hbpXE , GetEventWeight());
4660  }
4661 
4662  // Several UE calculation, not sure it is useful
4663  // with partical calorimter acceptance
4664  if(fMakeSeveralUE) FillNeutralUnderlyingEventSidesHistograms(ptTrig,pt,zT,hbpZT,deltaPhi);
4665 
4666  //
4667  // Decay photon correlations
4668  //
4669  if(fPi0Trigger && decayFound)
4670  FillDecayPhotonCorrelationHistograms(pt, phi, kFALSE) ;
4671 
4672  if(fFillAODWithReferences)
4673  {
4674  nrefs++;
4675  if(nrefs==1)
4676  {
4677  refpi0 = new TObjArray(0);
4678  refpi0->SetName(GetAODObjArrayName()+"Pi0s");
4679  refpi0->SetOwner(kFALSE);
4680  }
4681  refpi0->Add(pi0);
4682  } // put references in trigger AOD
4683 
4684  AliDebug(1,Form("Selected pi0: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
4685 
4686  } // loop
4687 
4688  // Fill AOD with reference tracks, if not filling histograms
4689  if(fFillAODWithReferences && refpi0)
4690  {
4691  aodParticle->AddObjArray(refpi0);
4692  }
4693 }
4694 
4695 //__________________________________________________________________________________________________________________
4697 //__________________________________________________________________________________________________________________
4698 void AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int_t histoIndex, Bool_t lostDecayPair)
4699 {
4700  AliDebug(1,"Make trigger particle - charged hadron correlation in AOD MC level");
4701 
4702  if( label < 0 )
4703  {
4704  AliDebug(1,Form(" *** bad label ***: label %d", label));
4705  return;
4706  }
4707 
4708  // Do MC correlation for a given particle type range.
4709  // Types defined in GetMCTagHistogramIndex:
4710  // 0 direct gamma; 1 pi0; 2 pi0 decay; 3 eta decay; 4 other decay; 5 electron; 6 other (hadron)
4711  if(histoIndex < fMCGenTypeMin || histoIndex > fMCGenTypeMax) return ;
4712 
4713  AliStack * stack = 0x0 ;
4714  TParticle * primary = 0x0 ;
4715  TClonesArray * mcparticles = 0x0 ;
4716  AliAODMCParticle * aodprimary = 0x0 ;
4717 
4718  Double_t eprim = 0 ;
4719  Double_t ptprim = 0 ;
4720  Double_t phiprim = 0 ;
4721  Double_t etaprim = 0 ;
4722  Int_t nTracks = 0 ;
4723  Int_t iParticle = 0 ;
4724 
4725  Bool_t leadTrig = kTRUE;
4726 
4727  if( GetReader()->ReadStack() )
4728  {
4729  stack = GetMCStack() ;
4730  if(!stack)
4731  {
4732  AliFatal("Stack not available, is the MC handler called? STOP");
4733  return;
4734  }
4735 
4736  //nTracks = stack->GetNtrack() ;
4737  nTracks = stack->GetNprimary();
4738  if( label >= stack->GetNtrack() )
4739  {
4740  if(GetDebug() > 2)
4741  AliInfo(Form("*** large label ***: label %d, n tracks %d", label, stack->GetNtrack()));
4742  return ;
4743  }
4744 
4745  primary = stack->Particle(label);
4746  if ( !primary )
4747  {
4748  AliInfo(Form(" *** no primary ***: label %d", label));
4749  return;
4750  }
4751 
4752  eprim = primary->Energy();
4753  ptprim = primary->Pt();
4754  etaprim = primary->Eta();
4755  phiprim = primary->Phi();
4756  if(phiprim < 0) phiprim+=TMath::TwoPi();
4757 
4758  if(ptprim < 0.01 || eprim < 0.01) return ;
4759 
4760  for (iParticle = 0 ; iParticle < nTracks ; iParticle++)
4761  {
4762  TParticle * particle = stack->Particle(iParticle);
4763 
4764  //keep only final state particles
4765  if( particle->GetStatusCode() != 1 ) continue ;
4766 
4767  //---------- Charged particles ----------------------
4768  Int_t pdg = particle->GetPdgCode();
4769  Int_t charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
4770  if(charge == 0) continue;
4771 
4772  particle->Momentum(fMomentum);
4773 
4774  //Particles in CTS acceptance, make sure to use the same selection as in the reader
4775  Bool_t inCTS = GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4776  //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4777  if( !inCTS ) continue;
4778 
4779  // Remove conversions
4780  if ( TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode() == 22 ) continue ;
4781 
4782  if ( label == iParticle ) continue; // avoid trigger particle
4783 
4784  Float_t phi = particle->Phi();
4785  if(phi < 0) phi+=TMath::TwoPi();
4786 
4787  Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim,histoIndex,lostDecayPair);
4788 
4789  if ( !lead ) leadTrig = kFALSE;
4790 
4791  // if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading) ) return;
4792  } // track loop
4793  } // ESD MC
4794  else if( GetReader()->ReadAODMCParticles() )
4795  {
4796  // Get the list of MC particles
4797  mcparticles = GetReader()->GetAODMCParticles();
4798  if( !mcparticles ) return;
4799 
4800  nTracks = mcparticles->GetEntriesFast() ;
4801 
4802  if( label >= nTracks )
4803  {
4804  if(GetDebug() > 2)
4805  AliInfo(Form(" *** large label ***: label %d, n tracks %d", label,nTracks));
4806  return;
4807  }
4808 
4809  // Get the particle
4810  aodprimary = (AliAODMCParticle*) mcparticles->At(label);
4811  if( !aodprimary )
4812  {
4813  AliInfo(Form(" *** no AOD primary ***: label %d", label));
4814  return;
4815  }
4816 
4817  eprim = aodprimary->E();
4818  ptprim = aodprimary->Pt();
4819  etaprim = aodprimary->Eta();
4820  phiprim = aodprimary->Phi();
4821  if(phiprim < 0) phiprim+=TMath::TwoPi();
4822 
4823  if(ptprim < 0.01 || eprim < 0.01) return ;
4824 
4825  for (iParticle = 0; iParticle < nTracks; iParticle++)
4826  {
4827  AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(iParticle);
4828 
4829  if (!part->IsPhysicalPrimary() ) continue; // same as part->GetStatus() !=1
4830 
4831  if ( part->Charge() == 0 ) continue;
4832 
4833  fMomentum.SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->E());
4834 
4835  //Particles in CTS acceptance, make sure to use the same selection as in the reader
4836  Bool_t inCTS = GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4837  //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4838  if( !inCTS ) continue;
4839 
4840  // Remove conversions
4841  Int_t indexmother = part->GetMother();
4842  if ( indexmother > -1 )
4843  {
4844  Int_t pdg = part->GetPdgCode();
4845  Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
4846  if (TMath::Abs(pdg) == 11 && mPdg == 22) continue;
4847  }
4848 
4849  if ( label == iParticle ) continue; // avoid trigger particle
4850 
4851  Float_t phi = part->Phi();
4852  if(phi < 0) phi+=TMath::TwoPi();
4853 
4854  Bool_t lead = Fil