60c4d8428372cc90c0e2ad05a37f848166e5882b
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleIsolation.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes hereby granted      *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // Class for analysis of particle isolation
18 // Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)
19 //
20 //  Class created from old AliPHOSGammaJet
21 //  (see AliRoot versions previous Release 4-09)
22 //
23 // -- Author: Gustavo Conesa (LNF-INFN)
24
25 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
26 //////////////////////////////////////////////////////////////////////////////
27
28 // --- ROOT system ---
29 #include <TClonesArray.h>
30 #include <TList.h>
31 #include <TObjString.h>
32 #include <TH2F.h>
33 #include <TClass.h>
34 #include <TH2F.h>
35 #include "TParticle.h"
36 #include "TDatabasePDG.h"
37
38 // --- Analysis system ---
39 #include "AliAnaParticleIsolation.h"
40 #include "AliCaloTrackReader.h"
41 #include "AliStack.h"
42 #include "AliIsolationCut.h"
43 #include "AliFiducialCut.h"
44 #include "AliMCAnalysisUtils.h"
45 #include "AliNeutralMesonSelection.h"
46 #include "AliAODMCParticle.h"
47 #include "AliAODPWG4ParticleCorrelation.h"
48 #include "AliMCAnalysisUtils.h"
49 #include "AliVTrack.h"
50 #include "AliVCluster.h"
51 #include "AliESDEvent.h"
52 #include "AliAODEvent.h"
53 // --- Detectors ---
54 #include "AliEMCALGeometry.h"
55 #include "AliPHOSGeoUtils.h"
56
57 ClassImp(AliAnaParticleIsolation)
58
59 //______________________________________________________________________________
60 AliAnaParticleIsolation::AliAnaParticleIsolation() :
61 AliAnaCaloTrackCorrBaseClass(),
62 fIsoDetector(-1),                 fIsoDetectorString(""),
63 fReMakeIC(0),                     fMakeSeveralIC(0),
64 fFillTMHisto(0),                  fFillSSHisto(1),
65 fFillUEBandSubtractHistograms(1), fFillCellHistograms(0),
66 fFillTaggedDecayHistograms(0),
67 fNDecayBits(0),                   fDecayBits(),
68 fFillNLMHistograms(0),
69 fLeadingOnly(0),                  fCheckLeadingWithNeutralClusters(0),
70 fSelectPrimariesInCone(0),        fMakePrimaryPi0DecayStudy(0),
71 fFillBackgroundBinHistograms(0),  fNBkgBin(0),
72 fFillPtTrigBinHistograms(0),      fNPtTrigBin(0),
73 fMinCellsAngleOverlap(0),
74 // Several IC
75 fNCones(0),                       fNPtThresFrac(0),
76 fConeSizes(),                     fPtThresholds(),
77 fPtFractions(),                   fSumPtThresholds(),
78 fMomentum(),                      fMomIso(),
79 fMomDaugh1(),                     fMomDaugh2(),
80 fTrackVector(),
81 // Histograms
82 fhEIso(0),                        fhPtIso(0),
83 fhPtCentralityIso(0),             fhPtEventPlaneIso(0),
84 fhPtNLocMaxIso(0),
85 fhPhiIso(0),                      fhEtaIso(0),                              fhEtaPhiIso(0),
86 fhEtaPhiNoIso(0),
87 fhENoIso(0),                      fhPtNoIso(0),                             fhPtNLocMaxNoIso(0),
88 fhPtInCone(0),
89 fhPtClusterInCone(0),             fhPtCellInCone(0),                        fhPtTrackInCone(0),
90 fhPtTrackInConeOtherBC(0),        fhPtTrackInConeOtherBCPileUpSPD(0),
91 fhPtTrackInConeBC0(0),            fhPtTrackInConeVtxBC0(0),
92 fhPtTrackInConeBC0PileUpSPD(0),
93 fhPtInConePileUp(),               fhPtInConeCent(0),
94 fhPerpConeSumPt(0),               fhPtInPerpCone(0),
95 fhEtaPhiInConeCluster(0),         fhEtaPhiCluster(0),
96 fhEtaPhiInConeTrack(0),           fhEtaPhiTrack(0),
97 fhEtaBandCluster(0),              fhPhiBandCluster(0),
98 fhEtaBandTrack(0),                fhPhiBandTrack(0),
99 fhEtaBandCell(0),                 fhPhiBandCell(0),
100 fhConePtLead(0),                  fhConePtLeadCluster(0),                   fhConePtLeadTrack(0),
101 fhConePtLeadClustervsTrack(0),    fhConePtLeadClusterTrackFrac(0),
102 fhConeSumPt(0),                   fhConeSumPtCellTrack(0),
103 fhConeSumPtCell(0),               fhConeSumPtCluster(0),                    fhConeSumPtTrack(0),
104 fhConeSumPtEtaBandUECluster(0),             fhConeSumPtPhiBandUECluster(0),
105 fhConeSumPtEtaBandUETrack(0),               fhConeSumPtPhiBandUETrack(0),
106 fhConeSumPtEtaBandUECell(0),                fhConeSumPtPhiBandUECell(0),
107 fhConeSumPtTrigEtaPhi(0),
108 fhConeSumPtCellTrackTrigEtaPhi(0),
109 fhConeSumPtEtaBandUEClusterTrigEtaPhi(0),   fhConeSumPtPhiBandUEClusterTrigEtaPhi(0),
110 fhConeSumPtEtaBandUETrackTrigEtaPhi(0),     fhConeSumPtPhiBandUETrackTrigEtaPhi(0),
111 fhConeSumPtEtaBandUECellTrigEtaPhi(0),      fhConeSumPtPhiBandUECellTrigEtaPhi(0),
112 fhConeSumPtEtaUESub(0),                     fhConeSumPtPhiUESub(0),
113 fhConeSumPtEtaUESubTrigEtaPhi(0),           fhConeSumPtPhiUESubTrigEtaPhi(0),
114 fhConeSumPtEtaUESubTrackCell(0),            fhConeSumPtPhiUESubTrackCell(0),
115 fhConeSumPtEtaUESubTrackCellTrigEtaPhi(0),  fhConeSumPtPhiUESubTrackCellTrigEtaPhi(0),
116 fhConeSumPtEtaUESubCluster(0),              fhConeSumPtPhiUESubCluster(0),
117 fhConeSumPtEtaUESubClusterTrigEtaPhi(0),    fhConeSumPtPhiUESubClusterTrigEtaPhi(0),
118 fhConeSumPtEtaUESubCell(0),                 fhConeSumPtPhiUESubCell(0),
119 fhConeSumPtEtaUESubCellTrigEtaPhi(0),       fhConeSumPtPhiUESubCellTrigEtaPhi(0),
120 fhConeSumPtEtaUESubTrack(0),                fhConeSumPtPhiUESubTrack(0),
121 fhConeSumPtEtaUESubTrackTrigEtaPhi(0),      fhConeSumPtPhiUESubTrackTrigEtaPhi(0),
122 fhFractionTrackOutConeEta(0),               fhFractionTrackOutConeEtaTrigEtaPhi(0),
123 fhFractionClusterOutConeEta(0),             fhFractionClusterOutConeEtaTrigEtaPhi(0),
124 fhFractionClusterOutConePhi(0),             fhFractionClusterOutConePhiTrigEtaPhi(0),
125 fhFractionCellOutConeEta(0),                fhFractionCellOutConeEtaTrigEtaPhi(0),
126 fhFractionCellOutConePhi(0),                fhFractionCellOutConePhiTrigEtaPhi(0),
127 fhConeSumPtClustervsTrack(0),               fhConeSumPtClusterTrackFrac(0),
128 fhConeSumPtEtaUESubClustervsTrack(0),       fhConeSumPtPhiUESubClustervsTrack(0),
129 fhConeSumPtCellvsTrack(0),
130 fhConeSumPtEtaUESubCellvsTrack(0),          fhConeSumPtPhiUESubCellvsTrack(0),
131 fhEtaBandClustervsTrack(0),                 fhPhiBandClustervsTrack(0),
132 fhEtaBandNormClustervsTrack(0),             fhPhiBandNormClustervsTrack(0),
133 fhEtaBandCellvsTrack(0),                    fhPhiBandCellvsTrack(0),
134 fhEtaBandNormCellvsTrack(0),                fhPhiBandNormCellvsTrack(0),
135 fhConeSumPtSubvsConeSumPtTotPhiTrack(0),    fhConeSumPtSubNormvsConeSumPtTotPhiTrack(0),
136 fhConeSumPtSubvsConeSumPtTotEtaTrack(0),    fhConeSumPtSubNormvsConeSumPtTotEtaTrack(0),
137 fhConeSumPtSubvsConeSumPtTotPhiCluster(0),  fhConeSumPtSubNormvsConeSumPtTotPhiCluster(0),
138 fhConeSumPtSubvsConeSumPtTotEtaCluster(0),  fhConeSumPtSubNormvsConeSumPtTotEtaCluster(0),
139 fhConeSumPtSubvsConeSumPtTotPhiCell(0),     fhConeSumPtSubNormvsConeSumPtTotPhiCell(0),
140 fhConeSumPtSubvsConeSumPtTotEtaCell(0),     fhConeSumPtSubNormvsConeSumPtTotEtaCell(0),
141 fhConeSumPtVSUETracksEtaBand(0),            fhConeSumPtVSUETracksPhiBand(0),
142 fhConeSumPtVSUEClusterEtaBand(0),           fhConeSumPtVSUEClusterPhiBand(0),
143 fhPtPrimMCPi0DecayPairOutOfCone(0),
144 fhPtPrimMCPi0DecayPairOutOfAcceptance(0),
145 fhPtPrimMCPi0DecayPairOutOfAcceptanceNoOverlap(0),
146 fhPtPrimMCPi0DecayPairAcceptInConeLowPt(0),
147 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlap(0),
148 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlapCaloE(0),
149 fhPtPrimMCPi0DecayPairNoOverlap(0),
150 fhPtPrimMCPi0DecayIsoPairOutOfCone(0),
151 fhPtPrimMCPi0DecayIsoPairOutOfAcceptance(0),
152 fhPtPrimMCPi0DecayIsoPairOutOfAcceptanceNoOverlap(0),
153 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPt(0),
154 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlap(0),
155 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlapCaloE(0),
156 fhPtPrimMCPi0DecayIsoPairNoOverlap(0),
157 fhPtPrimMCPi0Overlap(0),                    fhPtPrimMCPi0IsoOverlap(0),
158 fhPtLeadConeBin(0),                         fhSumPtConeBin(0),
159 fhPtLeadConeBinMC(0),                       fhSumPtConeBinMC(0),
160 fhPtLeadConeBinDecay(0),                    fhSumPtConeBinDecay(0),
161 fhPtLeadConeBinLambda0(0),                  fhSumPtConeBinLambda0(0),
162 fhPtLeadConeBinLambda0MC(0),                fhSumPtConeBinLambda0MC(0),
163 fhPtTrigBinPtLeadCone(0),                   fhPtTrigBinSumPtCone(0),
164 fhPtTrigBinPtLeadConeMC(0),                 fhPtTrigBinSumPtConeMC(0),
165 fhPtTrigBinPtLeadConeDecay(0),              fhPtTrigBinSumPtConeDecay(0),
166 fhPtTrigBinLambda0vsPtLeadCone(0),          fhPtTrigBinLambda0vsSumPtCone(0),
167 fhPtTrigBinLambda0vsPtLeadConeMC(0),        fhPtTrigBinLambda0vsSumPtConeMC(0),
168 // Number of local maxima in cluster
169 fhNLocMax(),
170 fhELambda0LocMax1(),              fhELambda1LocMax1(),
171 fhELambda0LocMax2(),              fhELambda1LocMax2(),
172 fhELambda0LocMaxN(),              fhELambda1LocMaxN(),
173 // PileUp
174 fhEIsoPileUp(),                   fhPtIsoPileUp(),
175 fhENoIsoPileUp(),                 fhPtNoIsoPileUp(),
176 fhTimeENoCut(0),                  fhTimeESPD(0),                  fhTimeESPDMulti(0),
177 fhTimeNPileUpVertSPD(0),          fhTimeNPileUpVertTrack(0),
178 fhTimeNPileUpVertContributors(0),
179 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
180 {
181   //default ctor
182   
183   //Initialize parameters
184   InitParameters();
185   
186   for(Int_t i = 0; i < 5 ; i++)
187   {
188     fConeSizes[i]      = 0 ;
189     
190     for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
191       fhSumPtLeadingPtMC[imc][i] = 0 ;
192     
193     for(Int_t j = 0; j < 5 ; j++)
194     {
195       fhPtThresIsolated             [i][j] = 0 ;
196       fhPtFracIsolated              [i][j] = 0 ;
197       fhSumPtIsolated               [i][j] = 0 ;
198       
199       fhEtaPhiPtThresIso            [i][j] = 0 ;
200       fhEtaPhiPtThresDecayIso       [i][j] = 0 ;
201       fhPtPtThresDecayIso           [i][j] = 0 ;
202       
203       fhEtaPhiPtFracIso             [i][j] = 0 ;
204       fhEtaPhiPtFracDecayIso        [i][j] = 0 ;
205       fhPtPtFracDecayIso            [i][j] = 0 ;
206       fhPtPtSumDecayIso             [i][j] = 0 ;
207       fhPtSumDensityIso             [i][j] = 0 ;
208       fhPtSumDensityDecayIso        [i][j] = 0 ;
209       fhEtaPhiSumDensityIso         [i][j] = 0 ;
210       fhEtaPhiSumDensityDecayIso    [i][j] = 0 ;
211       fhPtFracPtSumIso              [i][j] = 0 ;
212       fhPtFracPtSumDecayIso         [i][j] = 0 ;
213       fhEtaPhiFracPtSumIso          [i][j] = 0 ;
214       fhEtaPhiFracPtSumDecayIso     [i][j] = 0 ;
215       
216       for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
217       {
218         fhPtThresIsolatedMC[imc][i][j] = 0 ;
219         fhPtFracIsolatedMC [imc][i][j] = 0 ;
220         fhSumPtIsolatedMC  [imc][i][j] = 0 ;
221         
222       }
223     }
224   }
225   
226   for(Int_t ibit =0; ibit < 4; ibit++)
227   {
228     for(Int_t iso =0; iso < 2; iso++)
229     {
230       fhPtDecay       [iso][ibit] = 0;
231       fhEtaPhiDecay   [iso][ibit] = 0;
232       fhPtLambda0Decay[iso][ibit] = 0;
233       for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
234         fhPtDecayMC[iso][ibit][imc]    = 0;
235     }
236   }
237   
238   for(Int_t i = 0; i < 5 ; i++)
239   {
240     fPtFractions    [i] = 0 ;
241     fPtThresholds   [i] = 0 ;
242     fSumPtThresholds[i] = 0 ;
243     
244     fhSumPtLeadingPt    [i] = 0 ;
245     fhPtLeadingPt       [i] = 0 ;
246     fhPerpSumPtLeadingPt[i] = 0 ;
247     fhPerpPtLeadingPt   [i] = 0 ;
248   }
249   
250   for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
251   {
252     fhPtNoIsoMC  [imc]    = 0;
253     fhPtIsoMC    [imc]    = 0;
254     fhPhiIsoMC   [imc]    = 0;
255     fhEtaIsoMC   [imc]    = 0;
256     fhPtLambda0MC[imc][0] = 0;
257     fhPtLambda0MC[imc][1] = 0;
258   }
259   
260   for(Int_t i = 0; i < 2 ; i++)
261   {
262     fhTrackMatchedDEta[i] = 0 ;             fhTrackMatchedDPhi[i] = 0 ;   fhTrackMatchedDEtaDPhi  [i] = 0 ;
263     fhdEdx            [i] = 0 ;             fhEOverP          [i] = 0 ;   fhTrackMatchedMCParticle[i] = 0 ;
264     fhELambda0        [i] = 0 ;             fhPtLambda0       [i] = 0 ;   //fhELambda1        [i] = 0 ;
265     fhELambda0TRD     [i] = 0 ;             fhPtLambda0TRD    [i] = 0 ;   //fhELambda1TRD     [i] = 0 ;
266     
267     // Number of local maxima in cluster
268     fhNLocMax        [i] = 0 ;
269     fhELambda0LocMax1[i] = 0 ;              fhELambda1LocMax1[i] = 0 ;
270     fhELambda0LocMax2[i] = 0 ;              fhELambda1LocMax2[i] = 0 ;
271     fhELambda0LocMaxN[i] = 0 ;              fhELambda1LocMaxN[i] = 0 ;
272   }
273   
274   // Acceptance
275   for(Int_t i = 0; i < fgkNmcPrimTypes; i++)
276   {
277     fhPtPrimMCiso[i] = 0;
278     fhEPrimMC    [i] = 0;
279     fhPtPrimMC   [i] = 0;
280     fhEtaPrimMC  [i] = 0;
281     fhPhiPrimMC  [i] = 0;
282   }
283   
284   // Pile-Up
285   
286   for(Int_t i = 0 ; i < 7 ; i++)
287   {
288     fhPtInConePileUp[i] = 0 ;
289     fhEIsoPileUp    [i] = 0 ;
290     fhPtIsoPileUp   [i] = 0 ;
291     fhENoIsoPileUp  [i] = 0 ;
292     fhPtNoIsoPileUp [i] = 0 ;
293   }
294 }
295
296 //_______________________________________________________________________________________________
297 void AliAnaParticleIsolation::CalculateCaloUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
298                                                   Float_t & etaBandPtSum, Float_t & phiBandPtSum)
299 {
300   // Get the clusters pT or sum of pT in phi/eta bands or at 45 degrees from trigger
301   
302   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
303   
304   Float_t conesize   = GetIsolationCut()->GetConeSize();
305   
306   //Select the Calorimeter
307   TObjArray * pl = 0x0;
308   if      (GetCalorimeter() == kPHOS )
309     pl    = GetPHOSClusters();
310   else if (GetCalorimeter() == kEMCAL)
311     pl    = GetEMCALClusters();
312   
313   if(!pl) return ;
314   
315   //Get vertex for cluster momentum calculation
316   Double_t vertex[] = {0,0,0} ; //vertex ;
317   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
318     GetReader()->GetVertex(vertex);
319   
320   Float_t ptTrig    = pCandidate->Pt() ;
321   Float_t phiTrig   = pCandidate->Phi();
322   Float_t etaTrig   = pCandidate->Eta();
323   
324   for(Int_t icluster=0; icluster < pl->GetEntriesFast(); icluster++)
325   {
326     AliVCluster* cluster = (AliVCluster *) pl->At(icluster);
327     
328     if(!cluster)
329     {
330       printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Cluster not available?");
331       continue;
332     }
333     
334     //Do not count the candidate (photon or pi0) or the daughters of the candidate
335     if(cluster->GetID() == pCandidate->GetCaloLabel(0) ||
336        cluster->GetID() == pCandidate->GetCaloLabel(1)   ) continue ;
337     
338     //Remove matched clusters to tracks if Neutral and Track info is used
339     if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged &&
340        IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ;
341     
342     cluster->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
343     
344     //exclude particles in cone
345     Float_t rad = GetIsolationCut()->Radius(etaTrig, phiTrig, fMomentum.Eta(), fMomentum.Phi());
346     
347     // histo of eta and phi for all clusters
348     fhEtaPhiCluster->Fill(fMomentum.Eta(), fMomentum.Phi());
349     if(rad < conesize) {
350         // histos for all clusters in cone
351       fhEtaPhiInConeCluster->Fill(fMomentum.Eta(), fMomentum.Phi());
352       continue ;
353     }
354     //fill histogram for UE in phi band in EMCal acceptance
355     if(fMomentum.Eta() > (etaTrig-conesize) && fMomentum.Eta()  < (etaTrig+conesize))
356     {
357       phiBandPtSum+=fMomentum.Pt();
358       fhPhiBandCluster->Fill(fMomentum.Eta(),fMomentum.Phi());
359       
360     }
361     
362     //fill histogram for UE in eta band in EMCal acceptance
363     if(fMomentum.Phi() > (phiTrig-conesize) && fMomentum.Phi() < (phiTrig+conesize))
364     {
365       etaBandPtSum+=fMomentum.Pt();
366       fhEtaBandCluster->Fill(fMomentum.Eta(),fMomentum.Phi());
367     }
368   }
369   
370   fhConeSumPtEtaBandUECluster          ->Fill(ptTrig  ,       etaBandPtSum);
371   fhConeSumPtPhiBandUECluster          ->Fill(ptTrig  ,       phiBandPtSum);
372   fhConeSumPtEtaBandUEClusterTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSum);
373   fhConeSumPtPhiBandUEClusterTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSum);
374   
375 }
376
377 //________________________________________________________________________________________________
378 void AliAnaParticleIsolation::CalculateCaloCellUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
379                                                       Float_t & etaBandPtSumCells, Float_t & phiBandPtSumCells)
380 {
381   // Get the cells amplitude or sum of amplitude in phi/eta bands or at 45 degrees from trigger
382   
383   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
384   
385   Float_t conesize = GetIsolationCut()->GetConeSize();
386   
387   Float_t phiTrig = pCandidate->Phi();
388   if(phiTrig<0) phiTrig += TMath::TwoPi();
389   Float_t etaTrig = pCandidate->Eta();
390   
391   if(pCandidate->GetDetectorTag()==kEMCAL)
392   {
393     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
394     Int_t absId = -999;
395     
396     if (eGeom->GetAbsCellIdFromEtaPhi(etaTrig,phiTrig,absId))
397     {
398       if(!eGeom->CheckAbsCellId(absId)) return ;
399       
400       // Get absolute (col,row) of trigger particle
401       Int_t nSupMod = eGeom->GetSuperModuleNumber(absId);
402       Int_t nModule = -1;
403       Int_t imEta=-1, imPhi=-1;
404       Int_t ieta =-1, iphi =-1;
405       
406       if (eGeom->GetCellIndex(absId,nSupMod,nModule,imPhi,imEta))
407       {
408         eGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,imPhi,imEta,iphi,ieta);
409         
410         Int_t colTrig = ieta;
411         if (nSupMod % 2) colTrig = AliEMCALGeoParams::fgkEMCALCols + ieta ;
412         Int_t rowTrig = iphi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
413         
414         Int_t sqrSize = int(conesize/0.0143);
415         
416         AliVCaloCells * cells = GetEMCALCells();
417         
418         Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 24*(16/3) 5 full-size Sectors (2 SM) + 1 one-third Sector (2 SM)
419         Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols;
420         //  printf("nTotalRows %i, nTotalCols %i\n",nTotalRows,nTotalCols);
421         // Loop on cells in eta band
422         
423                   Int_t irowmin = rowTrig-sqrSize;
424         if(irowmin<0) irowmin=0;
425         Int_t irowmax = rowTrig+sqrSize;
426         if(irowmax>AliEMCALGeoParams::fgkEMCALRows) irowmax=AliEMCALGeoParams::fgkEMCALRows;
427         
428         
429         for(Int_t irow = irowmin; irow <irowmax; irow++)
430         {
431           for(Int_t icol = 0; icol < nTotalCols; icol++)
432           {
433             Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
434             if(inSector==5) continue;
435             Int_t inSupMod = -1;
436             Int_t icolLoc  = -1;
437             if(icol < AliEMCALGeoParams::fgkEMCALCols)
438             {
439               inSupMod = 2*inSector + 1;
440               icolLoc  = icol;
441             }
442             else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
443             {
444               inSupMod = 2*inSector;
445               icolLoc  = icol-AliEMCALGeoParams::fgkEMCALCols;
446             }
447             
448             Int_t irowLoc  = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;
449             
450             // Exclude cells in cone
451             if(TMath::Abs(icol-colTrig) < sqrSize || TMath::Abs(irow-rowTrig) < sqrSize){
452               continue ;
453             }
454             Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
455             if(!eGeom->CheckAbsCellId(iabsId)) continue;
456             etaBandPtSumCells += cells->GetCellAmplitude(iabsId);
457             fhEtaBandCell->Fill(colTrig,rowTrig);
458             
459             //          printf("ETA inSupMod %i,irowLoc %i,icolLoc %i, iabsId %i, etaBandPtSumCells %f\n",nSupMod,irowLoc,icolLoc,iabsId,etaBandPtSumCells);
460           }
461         }
462         Int_t icolmin = colTrig-sqrSize;
463         if(icolmin<0) icolmin=0;
464         Int_t icolmax = colTrig+sqrSize;
465         if(icolmax>AliEMCALGeoParams::fgkEMCALCols) icolmax=AliEMCALGeoParams::fgkEMCALCols;
466               
467         // Loop on cells in phi band
468         for(Int_t icol = icolmin; icol < icolmax; icol++)
469         {
470           for(Int_t irow = 0; irow < nTotalRows; irow++)
471           {
472             Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
473             if(inSector==5) continue;
474             Int_t inSupMod = -1;
475             Int_t icolLoc  = -1;
476             //    printf("icol %i, irow %i, inSector %i\n",icol,irow ,inSector);
477             if(icol < AliEMCALGeoParams::fgkEMCALCols)
478             {
479               //        printf("icol < AliEMCALGeoParams::fgkEMCALCols %i\n",AliEMCALGeoParams::fgkEMCALCols );
480               inSupMod = 2*inSector + 1;
481               icolLoc  = icol;
482             }
483             else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
484             {
485               //      printf("icol > AliEMCALGeoParams::fgkEMCALCols -1 %i\n",AliEMCALGeoParams::fgkEMCALCols -1 );
486               inSupMod = 2*inSector;
487               icolLoc  = icol-AliEMCALGeoParams::fgkEMCALCols;
488             }
489             
490             Int_t irowLoc  = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;   // Stesso problema di sopra //
491             
492             // Exclude cells in cone
493             if(TMath::Abs(icol-colTrig) < sqrSize) {
494               //printf("TMath::Abs(icol-colTrig) %i < sqrSize %i\n",TMath::Abs(icol-colTrig) ,sqrSize);continue ;
495             }
496             if(TMath::Abs(irow-rowTrig) < sqrSize) {
497               //printf("TMath::Abs(irow-rowTrig) %i < sqrSize %i\n",TMath::Abs(irow-rowTrig) ,sqrSize);continue ;
498             }
499             
500             Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
501             if(!eGeom->CheckAbsCellId(iabsId)) {printf("!eGeom->CheckAbsCellId(iabsId=%i) inSupMod %i irowLoc %i icolLoc %i \n",iabsId,inSupMod, irowLoc, icolLoc);continue;}
502             phiBandPtSumCells += cells->GetCellAmplitude(iabsId);
503             fhPhiBandCell->Fill(colTrig,rowTrig);
504             //printf("inSupMod %i,irowLoc %i,icolLoc %i, iabsId %i, phiBandPtSumCells %f\n",nSupMod,irowLoc,icolLoc,iabsId,phiBandPtSumCells);
505           }
506         }
507       }
508     }
509   }
510   
511   Float_t ptTrig = pCandidate->Pt();
512   
513   fhConeSumPtEtaBandUECell          ->Fill(ptTrig ,        etaBandPtSumCells);
514   fhConeSumPtPhiBandUECell          ->Fill(ptTrig ,        phiBandPtSumCells);
515   fhConeSumPtEtaBandUECellTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSumCells);
516   fhConeSumPtPhiBandUECellTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSumCells);
517   
518 }
519
520 //________________________________________________________________________________________________
521 void AliAnaParticleIsolation::CalculateTrackUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
522                                                    Float_t & etaBandPtSum, Float_t & phiBandPtSum)
523 {
524   // Get the track pT or sum of pT in phi/eta bands or at 45 degrees from trigger
525   
526   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyNeutral ) return ;
527   
528   Float_t conesize   = GetIsolationCut()->GetConeSize();
529   
530   Double_t sumptPerp= 0. ;
531   Float_t ptTrig    = pCandidate->Pt() ;
532   Float_t phiTrig   = pCandidate->Phi();
533   Float_t etaTrig   = pCandidate->Eta();
534   
535   TObjArray * trackList   = GetCTSTracks() ;
536   for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
537   {
538     AliVTrack* track = (AliVTrack *) trackList->At(itrack);
539     
540     if(!track)
541     {
542       printf("AliAnaParticleIsolation::CalculateTrackUEBand() - Track not available?");
543       continue;
544     }
545     
546     //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
547     if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
548        track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)   ) continue ;
549     
550     // histo of eta:phi for all tracks
551     fhEtaPhiTrack->Fill(track->Eta(),track->Phi());
552     
553     //exclude particles in cone
554     Float_t rad = GetIsolationCut()->Radius(etaTrig, phiTrig, track->Eta(), track->Phi());
555     if(rad < conesize) {
556         // histo of eta:phi for all tracks in cone
557       fhEtaPhiInConeTrack->Fill(track->Eta(),track->Phi());
558       continue ;
559     }
560     
561     //fill histogram for UE in phi band
562     if(track->Eta() > (etaTrig-conesize) && track->Eta()  < (etaTrig+conesize))
563     {
564       phiBandPtSum+=track->Pt();
565       fhPhiBandTrack->Fill(track->Eta(),track->Phi());
566     }
567     
568     //fill histogram for UE in eta band in EMCal acceptance
569     if(track->Phi() > (phiTrig-conesize) && track->Phi() < (phiTrig+conesize))
570     {
571       etaBandPtSum+=track->Pt();
572       fhEtaBandTrack->Fill(track->Eta(),track->Phi());
573     }
574     
575     //fill the histograms at +-45 degrees in phi from trigger particle, perpedicular to trigger axis in phi
576     Double_t dPhi = phiTrig - track->Phi() + TMath::PiOver2();
577     Double_t dEta = etaTrig - track->Eta();
578     Double_t arg  = dPhi*dPhi + dEta*dEta;
579     if(TMath::Sqrt(arg) < conesize)
580     {
581       fhPtInPerpCone->Fill(ptTrig,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
582       sumptPerp+=track->Pt();
583     }
584     
585     dPhi = phiTrig - track->Phi() - TMath::PiOver2();
586     arg  = dPhi*dPhi + dEta*dEta;
587     if(TMath::Sqrt(arg) < conesize)
588     {
589       fhPtInPerpCone->Fill(ptTrig,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
590       sumptPerp+=track->Pt();
591     }
592   }
593   
594   fhPerpConeSumPt                    ->Fill(ptTrig ,        sumptPerp   );
595   fhConeSumPtEtaBandUETrack          ->Fill(ptTrig ,        etaBandPtSum);
596   fhConeSumPtPhiBandUETrack          ->Fill(ptTrig ,        phiBandPtSum);
597   fhConeSumPtEtaBandUETrackTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSum);
598   fhConeSumPtPhiBandUETrackTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSum);
599   
600 }
601
602
603
604 //_____________________________________________________________________________________________________________________________________
605 void AliAnaParticleIsolation::CalculateNormalizeUEBandPerUnitArea(AliAODPWG4ParticleCorrelation * pCandidate, Float_t coneptsumCluster,
606                                                                   Float_t coneptsumCell,          Float_t coneptsumTrack,
607                                                                   Float_t &etaBandptsumTrackNorm, Float_t &etaBandptsumClusterNorm)
608 {
609   //normalize phi/eta band per area unit
610   
611   Float_t etaUEptsumTrack   = 0 ;
612   Float_t phiUEptsumTrack   = 0 ;
613   Float_t etaUEptsumCluster = 0 ;
614   Float_t phiUEptsumCluster = 0 ;
615   Float_t etaUEptsumCell    = 0 ;
616   Float_t phiUEptsumCell    = 0 ;
617   
618   Int_t   partTypeInCone    = GetIsolationCut()->GetParticleTypeInCone();
619   
620   // Do the normalization
621   
622   Float_t conesize  = GetIsolationCut()->GetConeSize();
623   Float_t coneA     = conesize*conesize*TMath::Pi(); // A = pi R^2, isolation cone area
624   Float_t ptTrig    = pCandidate->Pt() ;
625   Float_t phiTrig   = pCandidate->Phi();
626   Float_t etaTrig   = pCandidate->Eta();
627   
628   
629   // ------ //
630   // Tracks //
631   // ------ //
632   Float_t phiUEptsumTrackNorm  = 0 ;
633   Float_t etaUEptsumTrackNorm  = 0 ;
634   Float_t coneptsumTrackSubPhi = 0 ;
635   Float_t coneptsumTrackSubEta = 0 ;
636   Float_t coneptsumTrackSubPhiNorm = 0 ;
637   Float_t coneptsumTrackSubEtaNorm = 0 ;
638   etaBandptsumTrackNorm = 0 ;
639   
640   if( partTypeInCone!=AliIsolationCut::kOnlyNeutral )
641   {
642     // Sum the pT in the phi or eta band for clusters or tracks
643     CalculateTrackUEBand   (pCandidate,etaUEptsumTrack  ,phiUEptsumTrack  );// rajouter ici l'histo eta phi
644     
645     //Fill histos
646     fhConeSumPtVSUETracksEtaBand->Fill(coneptsumTrack,etaUEptsumTrack);
647     fhConeSumPtVSUETracksPhiBand->Fill(coneptsumTrack,phiUEptsumTrack);
648     
649     
650     Float_t correctConeSumTrack    = 1;
651     Float_t correctConeSumTrackPhi = 1;
652     
653     GetIsolationCut()->CalculateUEBandTrackNormalization(GetReader(),etaTrig, phiTrig,
654                                                          phiUEptsumTrack,etaUEptsumTrack,
655                                                          phiUEptsumTrackNorm,etaUEptsumTrackNorm,
656                                                          correctConeSumTrack,correctConeSumTrackPhi);
657     
658     coneptsumTrackSubPhi = coneptsumTrack - phiUEptsumTrackNorm;
659     coneptsumTrackSubEta = coneptsumTrack - etaUEptsumTrackNorm;
660     
661     etaBandptsumTrackNorm = etaUEptsumTrackNorm;
662     
663     fhConeSumPtPhiUESubTrack           ->Fill(ptTrig ,          coneptsumTrackSubPhi);
664     fhConeSumPtPhiUESubTrackTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumTrackSubPhi);
665     fhConeSumPtEtaUESubTrack           ->Fill(ptTrig ,          coneptsumTrackSubEta);
666     fhConeSumPtEtaUESubTrackTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumTrackSubEta);
667     
668     fhFractionTrackOutConeEta          ->Fill(ptTrig ,         correctConeSumTrack-1);
669     fhFractionTrackOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig,correctConeSumTrack-1);
670     
671     if(coneptsumTrack > 0)
672     {
673       coneptsumTrackSubPhiNorm = coneptsumTrackSubPhi/coneptsumTrack;
674         coneptsumTrackSubEtaNorm = coneptsumTrackSubEta/coneptsumTrack;
675     }
676     
677     fhConeSumPtSubvsConeSumPtTotPhiTrack    ->Fill(coneptsumTrack,coneptsumTrackSubPhi);
678     fhConeSumPtSubNormvsConeSumPtTotPhiTrack->Fill(coneptsumTrack,coneptsumTrackSubPhiNorm);
679     fhConeSumPtSubvsConeSumPtTotEtaTrack    ->Fill(coneptsumTrack,coneptsumTrackSubEta);
680     fhConeSumPtSubNormvsConeSumPtTotEtaTrack->Fill(coneptsumTrack,coneptsumTrackSubEtaNorm);
681     
682   }
683   
684   // ------------------------ //
685   // EMCal Clusters and cells //
686   // ------------------------ //
687   Float_t phiUEptsumClusterNorm  = 0 ;
688   Float_t etaUEptsumClusterNorm  = 0 ;
689   Float_t coneptsumClusterSubPhi = 0 ;
690   Float_t coneptsumClusterSubEta = 0 ;
691   Float_t coneptsumClusterSubPhiNorm = 0 ;
692   Float_t coneptsumClusterSubEtaNorm = 0 ;
693   Float_t phiUEptsumCellNorm     = 0 ;
694   Float_t etaUEptsumCellNorm     = 0 ;
695   Float_t coneptsumCellSubPhi    = 0 ;
696   Float_t coneptsumCellSubEta    = 0 ;
697   Float_t coneptsumCellSubPhiNorm = 0 ;
698   Float_t coneptsumCellSubEtaNorm = 0 ;
699   etaBandptsumClusterNorm = 0;
700   
701   if( partTypeInCone!=AliIsolationCut::kOnlyCharged )
702   {
703     
704     // -------------- //
705     // EMCal clusters //
706     // -------------- //
707     
708     // Sum the pT in the phi or eta band for clusters or tracks
709     CalculateCaloUEBand    (pCandidate,etaUEptsumCluster,phiUEptsumCluster);// rajouter ici l'histo eta phi
710     
711     //Fill histos
712     fhConeSumPtVSUEClusterEtaBand->Fill(coneptsumCluster,etaUEptsumCluster);
713     fhConeSumPtVSUEClusterPhiBand->Fill(coneptsumCluster,phiUEptsumCluster);
714     
715     
716     Float_t correctConeSumClusterEta = 1;
717     Float_t correctConeSumClusterPhi = 1;
718     
719     GetIsolationCut()->CalculateUEBandClusterNormalization(GetReader(),etaTrig, phiTrig,
720                                                            phiUEptsumCluster,etaUEptsumCluster,
721                                                            phiUEptsumClusterNorm,etaUEptsumClusterNorm,
722                                                            correctConeSumClusterEta,correctConeSumClusterPhi);
723     
724     // In case that cone is out of eta and phi side, we are over correcting, not too often with the current cuts ...
725     // Comment if not used
726     //  Float_t coneBadCellsCoeff   =1;
727     //  Float_t etaBandBadCellsCoeff=1;
728     //  Float_t phiBandBadCellsCoeff=1;
729     //  GetIsolationCut()->GetCoeffNormBadCell(pCandidate,   GetReader(),coneBadCellsCoeff,etaBandBadCellsCoeff,phiBandBadCellsCoeff) ;
730     
731     //coneptsumCluster=coneptsumCluster*coneBadCellsCoeff*correctConeSumClusterEta*correctConeSumClusterPhi;
732     
733     coneptsumClusterSubPhi = coneptsumCluster - phiUEptsumClusterNorm;
734     coneptsumClusterSubEta = coneptsumCluster - etaUEptsumClusterNorm;
735     
736     etaBandptsumClusterNorm = etaUEptsumClusterNorm;
737     
738     fhConeSumPtPhiUESubCluster           ->Fill(ptTrig ,          coneptsumClusterSubPhi);
739     fhConeSumPtPhiUESubClusterTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumClusterSubPhi);
740     fhConeSumPtEtaUESubCluster           ->Fill(ptTrig ,          coneptsumClusterSubEta);
741     fhConeSumPtEtaUESubClusterTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumClusterSubEta);
742     
743     fhFractionClusterOutConeEta          ->Fill(ptTrig ,          correctConeSumClusterEta-1);
744     fhFractionClusterOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumClusterEta-1);
745     fhFractionClusterOutConePhi          ->Fill(ptTrig ,          correctConeSumClusterPhi-1);
746     fhFractionClusterOutConePhiTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumClusterPhi-1);
747     
748     if(coneptsumCluster!=0)
749     {
750             coneptsumClusterSubPhiNorm = coneptsumClusterSubPhi/coneptsumCluster;
751         coneptsumClusterSubEtaNorm = coneptsumClusterSubEta/coneptsumCluster;
752     }
753     
754     fhConeSumPtSubvsConeSumPtTotPhiCluster    ->Fill(coneptsumCluster,coneptsumClusterSubPhi);
755     fhConeSumPtSubNormvsConeSumPtTotPhiCluster->Fill(coneptsumCluster,coneptsumClusterSubPhiNorm);
756     fhConeSumPtSubvsConeSumPtTotEtaCluster    ->Fill(coneptsumCluster,coneptsumClusterSubEta);
757     fhConeSumPtSubNormvsConeSumPtTotEtaCluster->Fill(coneptsumCluster,coneptsumClusterSubEtaNorm);
758     
759     // ----------- //
760     // EMCal Cells //
761     // ----------- //
762     
763     if(fFillCellHistograms)
764     {
765       // Sum the pT in the phi or eta band for clusters or tracks
766       CalculateCaloCellUEBand(pCandidate,etaUEptsumCell   ,phiUEptsumCell   );
767       
768       // Move to AliIsolationCut the calculation not the histograms??
769       
770       //Careful here if EMCal limits changed .. 2010 (4 SM) to 2011-12 (10 SM), for the moment consider 100 deg in phi
771       Float_t emcEtaSize = 0.7*2; // TO FIX
772       Float_t emcPhiSize = TMath::DegToRad()*100.; // TO FIX
773       
774       if(((2*conesize*emcPhiSize)-coneA)!=0)phiUEptsumCellNorm = phiUEptsumCell*(coneA / ((2*conesize*emcPhiSize)-coneA));
775       if(((2*conesize*emcEtaSize)-coneA)!=0)etaUEptsumCellNorm = etaUEptsumCell*(coneA / ((2*conesize*emcEtaSize)-coneA));
776       
777       // Need to correct coneptsumCluster by the fraction of the cone out of the calorimeter cut acceptance!
778       
779       Float_t correctConeSumCellEta = 1;
780       if(TMath::Abs(etaTrig)+conesize > emcEtaSize/2.)
781       {
782         Float_t excess = TMath::Abs(etaTrig) + conesize - emcEtaSize/2.;
783         correctConeSumCellEta = GetIsolationCut()->CalculateExcessAreaFraction(excess);
784         //printf("Excess EMC-Eta %2.3f, coneA %2.2f,  excessA %2.2f, angle %2.2f,factor %2.2f\n",excess,coneA, excessA, angle*TMath::RadToDeg(), correctConeSumClusterEta);
785         // Need to correct phi band surface if part of the cone falls out of track cut acceptance!
786         if(((2*(conesize-excess)*emcPhiSize)-(coneA-correctConeSumCellEta))!=0)phiUEptsumCellNorm = phiUEptsumCell*(coneA / ((2*(conesize-excess)*emcPhiSize)-(coneA-correctConeSumCellEta)));
787       }
788       
789       Float_t correctConeSumCellPhi = 1;
790       //printf("EMCPhiTrig %2.2f, conesize %2.2f, sum %2.2f, rest %2.2f \n",phiTrig*TMath::RadToDeg(),conesize*TMath::RadToDeg(),(phiTrig+conesize)*TMath::RadToDeg(),(phiTrig-conesize)*TMath::RadToDeg() );
791       if((phiTrig+conesize > 180*TMath::DegToRad()) ||
792          (phiTrig-conesize <  80*TMath::DegToRad()))
793       {
794         Float_t excess = 0;
795         if( phiTrig+conesize > 180*TMath::DegToRad() ) excess = conesize + phiTrig - 180*TMath::DegToRad() ;
796         else                                           excess = conesize - phiTrig +  80*TMath::DegToRad() ;
797         
798         correctConeSumCellPhi = GetIsolationCut()->CalculateExcessAreaFraction(excess);
799         //printf("Excess EMC-Phi %2.3f, coneA %2.2f,  excessA %2.2f, angle %2.2f,factor %2.2f\n",excess,coneA, excessA, angle*TMath::RadToDeg(), correctConeSumClusterPhi);
800         
801         // Need to correct eta band surface if part of the cone falls out of track cut acceptance!
802         if(((2*(conesize-excess)*emcEtaSize)-(coneA-correctConeSumCellPhi))!=0)etaUEptsumCellNorm = etaUEptsumCell*(coneA / ((2*(conesize-excess)*emcEtaSize)-(coneA-correctConeSumCellPhi)));
803         
804       }
805       
806       // In case that cone is out of eta and phi side, we are over correcting, not too often with the current cuts ...
807       coneptsumCellSubPhi = coneptsumCell*correctConeSumCellEta*correctConeSumCellPhi - phiUEptsumCellNorm;
808       coneptsumCellSubEta = coneptsumCell*correctConeSumCellEta*correctConeSumCellPhi - etaUEptsumCellNorm;
809       
810       fhConeSumPtPhiUESubCell           ->Fill(ptTrig ,          coneptsumCellSubPhi);
811       fhConeSumPtPhiUESubCellTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumCellSubPhi);
812       fhConeSumPtEtaUESubCell           ->Fill(ptTrig ,          coneptsumCellSubEta);
813       fhConeSumPtEtaUESubCellTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumCellSubEta);
814       
815       fhFractionCellOutConeEta          ->Fill(ptTrig ,          correctConeSumCellEta-1);
816       fhFractionCellOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumCellEta-1);
817       fhFractionCellOutConePhi          ->Fill(ptTrig ,          correctConeSumCellPhi-1);
818       fhFractionCellOutConePhiTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumCellPhi-1);
819       if(coneptsumCell!=0)
820       {
821         coneptsumCellSubPhiNorm = coneptsumCellSubPhi/coneptsumCell;
822         coneptsumCellSubEtaNorm = coneptsumCellSubEta/coneptsumCell;
823       }
824       
825       fhConeSumPtSubvsConeSumPtTotPhiCell    ->Fill(coneptsumCell,coneptsumCellSubPhi);
826       fhConeSumPtSubNormvsConeSumPtTotPhiCell->Fill(coneptsumCell,coneptsumCellSubPhiNorm);
827       fhConeSumPtSubvsConeSumPtTotEtaCell    ->Fill(coneptsumCell,coneptsumCellSubEta);
828       fhConeSumPtSubNormvsConeSumPtTotEtaCell->Fill(coneptsumCell,coneptsumCellSubEtaNorm);
829     }
830   }
831   
832   if( partTypeInCone==AliIsolationCut::kNeutralAndCharged )
833   {
834     // --------------------------- //
835     // Tracks and clusters in cone //
836     // --------------------------- //
837     
838     Double_t sumPhiUESub = coneptsumClusterSubPhi + coneptsumTrackSubPhi;
839     Double_t sumEtaUESub = coneptsumClusterSubEta + coneptsumTrackSubEta;
840     
841     fhConeSumPtPhiUESub          ->Fill(ptTrig ,          sumPhiUESub);
842     fhConeSumPtPhiUESubTrigEtaPhi->Fill(etaTrig, phiTrig, sumPhiUESub);
843     fhConeSumPtEtaUESub          ->Fill(ptTrig ,          sumEtaUESub);
844     fhConeSumPtEtaUESubTrigEtaPhi->Fill(etaTrig, phiTrig, sumEtaUESub);
845     
846     fhEtaBandClustervsTrack    ->Fill(etaUEptsumCluster    ,etaUEptsumTrack    );
847     fhPhiBandClustervsTrack    ->Fill(phiUEptsumCluster    ,phiUEptsumTrack    );
848     fhEtaBandNormClustervsTrack->Fill(etaUEptsumClusterNorm,etaUEptsumTrackNorm);
849     fhPhiBandNormClustervsTrack->Fill(phiUEptsumClusterNorm,phiUEptsumTrackNorm);
850     
851     fhConeSumPtEtaUESubClustervsTrack->Fill(coneptsumClusterSubEta,coneptsumTrackSubEta);
852     fhConeSumPtPhiUESubClustervsTrack->Fill(coneptsumClusterSubPhi,coneptsumTrackSubPhi);
853     
854     // ------------------------ //
855     // Tracks and cells in cone //
856     // ------------------------ //
857     
858     if(fFillCellHistograms)
859     {
860       Double_t sumPhiUESubTrackCell = coneptsumCellSubPhi + coneptsumTrackSubPhi;
861       Double_t sumEtaUESubTrackCell = coneptsumCellSubEta + coneptsumTrackSubEta;
862       
863       fhConeSumPtPhiUESubTrackCell          ->Fill(ptTrig ,          sumPhiUESubTrackCell);
864       fhConeSumPtPhiUESubTrackCellTrigEtaPhi->Fill(etaTrig, phiTrig, sumPhiUESubTrackCell);
865       fhConeSumPtEtaUESubTrackCell          ->Fill(ptTrig ,          sumEtaUESubTrackCell);
866       fhConeSumPtEtaUESubTrackCellTrigEtaPhi->Fill(etaTrig, phiTrig, sumEtaUESubTrackCell);
867       
868       fhEtaBandCellvsTrack    ->Fill(etaUEptsumCell    ,etaUEptsumTrack    );
869       fhPhiBandCellvsTrack    ->Fill(phiUEptsumCell    ,phiUEptsumTrack    );
870       fhEtaBandNormCellvsTrack->Fill(etaUEptsumCellNorm,etaUEptsumTrackNorm);
871       fhPhiBandNormCellvsTrack->Fill(phiUEptsumCellNorm,phiUEptsumTrackNorm);
872       
873       fhConeSumPtEtaUESubCellvsTrack->Fill(coneptsumCellSubEta,coneptsumTrackSubEta);
874       fhConeSumPtPhiUESubCellvsTrack->Fill(coneptsumCellSubPhi,coneptsumTrackSubPhi);
875     }
876     
877   }
878 }
879
880
881 //______________________________________________________________________________________________________________
882 void AliAnaParticleIsolation::CalculateCaloSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
883                                                         Float_t & coneptsumCluster, Float_t & coneptLeadCluster)
884 {
885   // Get the cluster pT or sum of pT in isolation cone
886   coneptLeadCluster = 0;
887   coneptsumCluster  = 0;
888   
889   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
890   
891   //Recover reference arrays with clusters and tracks
892   TObjArray * refclusters = aodParticle->GetObjArray(GetAODObjArrayName()+"Clusters");
893   if(!refclusters) return ;
894   
895   Float_t ptTrig = aodParticle->Pt();
896   
897   //Get vertex for cluster momentum calculation
898   Double_t vertex[] = {0,0,0} ; //vertex ;
899   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
900     GetReader()->GetVertex(vertex);
901   
902   Float_t ptcone = 0;
903   
904   for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
905   {
906     AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
907     calo->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
908     ptcone = fMomentum.Pt();
909     
910     fhPtInCone       ->Fill(ptTrig, ptcone);
911     fhPtClusterInCone->Fill(ptTrig, ptcone);
912     
913     if(IsPileUpAnalysisOn())
914     {
915       if(GetReader()->IsPileUpFromSPD())               fhPtInConePileUp[0]->Fill(ptTrig,ptcone);
916       if(GetReader()->IsPileUpFromEMCal())             fhPtInConePileUp[1]->Fill(ptTrig,ptcone);
917       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtInConePileUp[2]->Fill(ptTrig,ptcone);
918       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtInConePileUp[3]->Fill(ptTrig,ptcone);
919       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtInConePileUp[4]->Fill(ptTrig,ptcone);
920       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtInConePileUp[5]->Fill(ptTrig,ptcone);
921       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(ptTrig,ptcone);
922     }
923     
924     if(IsHighMultiplicityAnalysisOn()) fhPtInConeCent->Fill(GetEventCentrality(),ptcone);
925     
926     coneptsumCluster+=ptcone;
927     if(ptcone > coneptLeadCluster) coneptLeadCluster = ptcone;
928   }
929   
930   fhConeSumPtCluster ->Fill(ptTrig, coneptsumCluster );
931   fhConePtLeadCluster->Fill(ptTrig, coneptLeadCluster);
932 }
933
934 //______________________________________________________________________________________________________
935 void AliAnaParticleIsolation::CalculateCaloCellSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
936                                                             Float_t & coneptsumCell)
937 {
938   // Get the cell amplityde or sum of amplitudes in isolation cone
939   // Mising: Remove signal cells in cone in case the trigger is a cluster!
940   
941   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
942   
943   Float_t conesize = GetIsolationCut()->GetConeSize();
944   
945   Float_t  ptTrig  = aodParticle->Pt();
946   Float_t  phiTrig = aodParticle->Phi();
947   if(phiTrig<0) phiTrig += TMath::TwoPi();
948   Float_t  etaTrig = aodParticle->Eta();
949   
950   if(aodParticle->GetDetectorTag()==kEMCAL)
951   {
952     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
953     Int_t absId = -999;
954     
955     if (eGeom->GetAbsCellIdFromEtaPhi(etaTrig,phiTrig,absId))
956     {
957       if(!eGeom->CheckAbsCellId(absId)) return ;
958       
959       // Get absolute (col,row) of trigger particle
960       Int_t nSupMod = eGeom->GetSuperModuleNumber(absId);
961       Int_t nModule = -1;
962       Int_t imEta=-1, imPhi=-1;
963       Int_t ieta =-1, iphi =-1;
964       
965       if (eGeom->GetCellIndex(absId,nSupMod,nModule,imPhi,imEta))
966       {
967         Int_t iEta=-1, iPhi=-1;
968         eGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,imPhi,imEta,iphi,ieta);
969         
970         Int_t colTrig = iEta;
971         if (nSupMod % 2) colTrig = AliEMCALGeoParams::fgkEMCALCols + iEta ;
972         Int_t rowTrig = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
973         
974         Int_t sqrSize = int(conesize/0.0143);
975         
976         AliVCaloCells * cells = GetEMCALCells();
977         
978         // Loop on cells in cone
979         for(Int_t irow = rowTrig-sqrSize; irow < rowTrig+sqrSize; irow++)
980         {
981           for(Int_t icol = colTrig-sqrSize; icol < colTrig+sqrSize; icol++)
982           {
983             Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
984             if(inSector==5) continue;
985             
986             Int_t inSupMod = -1;
987             Int_t icolLoc  = -1;
988             if(icol < AliEMCALGeoParams::fgkEMCALCols)
989             {
990               inSupMod = 2*inSector + 1;
991               icolLoc  = icol;
992             }
993             else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
994             {
995               inSupMod = 2*inSector;
996               icolLoc  = icol-AliEMCALGeoParams::fgkEMCALCols;
997             }
998             
999             Int_t irowLoc  = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;
1000             
1001             Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
1002             if(!eGeom->CheckAbsCellId(iabsId)) continue;
1003             
1004             fhPtCellInCone->Fill(ptTrig, cells->GetCellAmplitude(iabsId));
1005             coneptsumCell += cells->GetCellAmplitude(iabsId);
1006           }
1007         }
1008       }
1009     }
1010   }
1011   
1012   fhConeSumPtCell->Fill(ptTrig,coneptsumCell);
1013   
1014 }
1015
1016 //___________________________________________________________________________________________________________
1017 void AliAnaParticleIsolation::CalculateTrackSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
1018                                                          Float_t & coneptsumTrack, Float_t & coneptLeadTrack)
1019 {
1020   // Get the track pT or sum of pT in isolation cone
1021   
1022   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyNeutral ) return ;
1023   
1024   //Recover reference arrays with clusters and tracks
1025   TObjArray * reftracks   = aodParticle->GetObjArray(GetAODObjArrayName()+"Tracks");
1026   if(!reftracks) return ;
1027   
1028   Float_t  ptTrig = aodParticle->Pt();
1029   Double_t bz     = GetReader()->GetInputEvent()->GetMagneticField();
1030   
1031   for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
1032   {
1033     AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
1034     Float_t pTtrack = track->Pt();
1035     
1036     fhPtInCone     ->Fill(ptTrig,pTtrack);
1037     fhPtTrackInCone->Fill(ptTrig,pTtrack);
1038     
1039     if(IsPileUpAnalysisOn())
1040     {
1041       ULong_t status = track->GetStatus();
1042       Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1043       //Double32_t tof = track->GetTOFsignal()*1e-3;
1044       Int_t trackBC = track->GetTOFBunchCrossing(bz);
1045       
1046       if     ( okTOF && trackBC!=0 ) fhPtTrackInConeOtherBC->Fill(ptTrig,pTtrack);
1047       else if( okTOF && trackBC==0 ) fhPtTrackInConeBC0    ->Fill(ptTrig,pTtrack);
1048       
1049       Int_t vtxBC = GetReader()->GetVertexBC();
1050       if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA) fhPtTrackInConeVtxBC0->Fill(ptTrig,pTtrack);
1051       
1052       if(GetReader()->IsPileUpFromSPD())             { fhPtInConePileUp[0]->Fill(ptTrig,pTtrack);
1053         if(okTOF && trackBC!=0 )                         fhPtTrackInConeOtherBCPileUpSPD->Fill(ptTrig,pTtrack);
1054         if(okTOF && trackBC==0 )                         fhPtTrackInConeBC0PileUpSPD    ->Fill(ptTrig,pTtrack); }
1055       if(GetReader()->IsPileUpFromEMCal())             fhPtInConePileUp[1]->Fill(ptTrig,pTtrack);
1056       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtInConePileUp[2]->Fill(ptTrig,pTtrack);
1057       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtInConePileUp[3]->Fill(ptTrig,pTtrack);
1058       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtInConePileUp[4]->Fill(ptTrig,pTtrack);
1059       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtInConePileUp[5]->Fill(ptTrig,pTtrack);
1060       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(ptTrig,pTtrack);
1061     }
1062     
1063     if(IsHighMultiplicityAnalysisOn()) fhPtInConeCent->Fill(GetEventCentrality(),pTtrack);
1064     
1065     coneptsumTrack+=pTtrack;
1066     if(pTtrack > coneptLeadTrack) coneptLeadTrack = pTtrack;
1067   }
1068
1069   fhConeSumPtTrack ->Fill(ptTrig, coneptsumTrack );
1070   fhConePtLeadTrack->Fill(ptTrig, coneptLeadTrack);
1071
1072 }
1073
1074 //_________________________________________________________________
1075 void AliAnaParticleIsolation::FillPileUpHistograms(Int_t clusterID)
1076 {
1077   // Fill some histograms to understand pile-up
1078   
1079   if(clusterID < 0 )
1080   {
1081     printf("AliAnaParticleIsolation::FillPileUpHistograms(), ID of cluster = %d, not possible! ", clusterID);
1082     return;
1083   }
1084   
1085   Int_t iclus = -1;
1086   TObjArray* clusters = 0x0;
1087   if     (GetCalorimeter() == kEMCAL) clusters = GetEMCALClusters();
1088   else if(GetCalorimeter() == kPHOS ) clusters = GetPHOSClusters();
1089   
1090   Float_t energy = 0;
1091   Float_t time   = -1000;
1092   
1093   if(clusters)
1094   {
1095     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
1096     energy = cluster->E();
1097     time   = cluster->GetTOF()*1e9;
1098   }
1099   
1100   //printf("E %f, time %f\n",energy,time);
1101   AliVEvent * event = GetReader()->GetInputEvent();
1102   
1103   fhTimeENoCut->Fill(energy,time);
1104   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
1105   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
1106   
1107   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
1108   
1109   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
1110   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
1111   
1112   // N pile up vertices
1113   Int_t nVerticesSPD    = -1;
1114   Int_t nVerticesTracks = -1;
1115   
1116   if      (esdEv)
1117   {
1118     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
1119     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
1120     
1121   }//ESD
1122   else if (aodEv)
1123   {
1124     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
1125     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
1126   }//AOD
1127   
1128   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
1129   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
1130   
1131   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
1132   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
1133   
1134   Int_t ncont = -1;
1135   Float_t z1 = -1, z2 = -1;
1136   Float_t diamZ = -1;
1137   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
1138   {
1139     if      (esdEv)
1140     {
1141       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
1142       ncont=pv->GetNContributors();
1143       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
1144       z2 = pv->GetZ();
1145       diamZ = esdEv->GetDiamondZ();
1146     }//ESD
1147     else if (aodEv)
1148     {
1149       AliAODVertex *pv=aodEv->GetVertex(iVert);
1150       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
1151       ncont=pv->GetNContributors();
1152       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
1153       z2=pv->GetZ();
1154       diamZ = aodEv->GetDiamondZ();
1155     }// AOD
1156     
1157     Double_t distZ  = TMath::Abs(z2-z1);
1158     diamZ  = TMath::Abs(z2-diamZ);
1159     
1160     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
1161     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
1162     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
1163     
1164   }// loop
1165 }
1166
1167 //_____________________________________________________________________________________________________________________
1168 void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(AliAODPWG4ParticleCorrelation  *pCandidate,
1169                                                                             Float_t coneptsum, Float_t coneleadpt,
1170                                                                             Int_t mcIndex)
1171 {
1172   // Fill Track matching and Shower Shape control histograms
1173   if(!fFillTMHisto && !fFillSSHisto && !fFillBackgroundBinHistograms && !fFillTaggedDecayHistograms) return;
1174   
1175   Int_t  clusterID = pCandidate->GetCaloLabel(0) ;
1176   Int_t  nMaxima   = pCandidate->GetNLM();
1177   Int_t  mcTag     = pCandidate->GetTag() ;
1178   Bool_t isolated  = pCandidate->IsIsolated();
1179   
1180   if(clusterID < 0 )
1181   {
1182     printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(), ID of cluster = %d, not possible! \n", clusterID);
1183     return;
1184   }
1185   
1186
1187   Float_t m02    = pCandidate->GetM02() ;
1188   Float_t energy = pCandidate->E();
1189   Float_t pt     = pCandidate->Pt();
1190   Float_t eta    = pCandidate->Eta();
1191   Float_t phi    = pCandidate->Phi();
1192   if(phi<0) phi+= TMath::TwoPi();
1193   
1194   // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
1195   if(fFillTaggedDecayHistograms)
1196   {
1197     Int_t decayTag = pCandidate->DecayTag();
1198     if(decayTag < 0) decayTag = 0;
1199
1200     for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1201     {
1202       if(!GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit])) continue;
1203       
1204       if(fFillSSHisto) fhPtLambda0Decay[isolated][ibit]->Fill(pt,m02);
1205       
1206       // In case it was not done on the trigger selection task
1207       // apply here a shower shape cut, not too strong, to select photons
1208       if( m02 < 0.3 ) continue;
1209       
1210       fhPtDecay    [isolated][ibit]->Fill(pt);
1211       fhEtaPhiDecay[isolated][ibit]->Fill(eta,phi);
1212      
1213       if(IsDataMC())
1214       {
1215         fhPtDecayMC[isolated][ibit][mcIndex]->Fill(pt);
1216
1217         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1218           fhPtDecayMC[isolated][ibit][kmcPhoton]->Fill(pt);
1219         
1220         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
1221           fhPtDecayMC[isolated][ibit][kmcPi0DecayLostPair]->Fill(pt);
1222       }
1223     } // bit loop
1224   } // decay histograms
1225
1226   
1227   // Get the max pt leading in cone or the sum of pt in cone
1228   // assign a bin to the candidate, depending on both quantities
1229   // see the shower shape in those bins.
1230   if(fFillBackgroundBinHistograms)
1231   {
1232     // Get the background bin for this cone and trigger
1233     Int_t ptsumBin  = -1;
1234     Int_t leadptBin = -1;
1235     
1236     if( GetDebug() > 1 )
1237       printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms() - pT cand: %2.2f, In cone pT: Sum %2.2f, Lead %2.2f, n bins %d\n",
1238              pt,coneptsum,coneleadpt,fNBkgBin);
1239     
1240     for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1241     {
1242       if( coneptsum  >= fBkgBinLimit[ibin] && coneptsum  < fBkgBinLimit[ibin+1]) ptsumBin  = ibin;
1243       if( coneleadpt >= fBkgBinLimit[ibin] && coneleadpt < fBkgBinLimit[ibin+1]) leadptBin = ibin;
1244     }
1245     
1246     // Fill the histograms per bin of pt lead or pt sum
1247     if( GetDebug() > 1 && ptsumBin  >=0 ) printf("\t Sum bin %d [%2.2f,%2.2f]\n" , ptsumBin ,fBkgBinLimit[ptsumBin] ,fBkgBinLimit[ptsumBin +1]);
1248     if( GetDebug() > 1 && leadptBin >=0 ) printf("\t Lead bin %d [%2.2f,%2.2f]\n", leadptBin,fBkgBinLimit[leadptBin],fBkgBinLimit[leadptBin+1]);
1249     
1250     if( leadptBin >=0 )
1251     {
1252       fhPtLeadConeBin[leadptBin]->Fill(pt);
1253       if(fFillSSHisto) fhPtLeadConeBinLambda0[leadptBin]->Fill(pt,m02);
1254     }
1255     
1256     if( ptsumBin  >=0 )
1257     {
1258       fhSumPtConeBin[ptsumBin]->Fill(pt);
1259       if(fFillSSHisto) fhSumPtConeBinLambda0[ptsumBin]->Fill(pt,m02);
1260     }
1261     
1262     // Check if it was a decay
1263     if(fFillTaggedDecayHistograms && m02 < 0.3)
1264     {
1265       Int_t decayTag = pCandidate->DecayTag();
1266       if(decayTag < 0) decayTag = 0;
1267       for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1268       {
1269         if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
1270         {
1271           Int_t leadptBinDecay = leadptBin+ibit*fNBkgBin;
1272           Int_t  ptsumBinDecay =  ptsumBin+ibit*fNBkgBin;
1273           if( leadptBin >=0 ) fhPtLeadConeBinDecay[leadptBinDecay]->Fill(pt);
1274           if( ptsumBin  >=0 ) fhSumPtConeBinDecay [ ptsumBinDecay]->Fill(pt);
1275         }
1276       }
1277     }
1278     
1279     if( GetDebug() > 1 && leadptBin == 0 )
1280       printf("No track/clusters in isolation cone: cand pt %2.2f GeV/c, track multiplicity %d, N clusters %d\n",
1281              pt, GetTrackMultiplicity(),GetEMCALClusters()->GetEntriesFast());
1282     
1283     if(IsDataMC())
1284     {
1285       Int_t leadptBinMC = leadptBin+mcIndex*fNBkgBin;
1286       Int_t  ptsumBinMC =  ptsumBin+mcIndex*fNBkgBin;
1287      
1288       if( leadptBin >=0 )
1289       {
1290         fhPtLeadConeBinMC[leadptBinMC]->Fill(pt);
1291         if(fFillSSHisto) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,m02);
1292       }
1293       
1294       if( ptsumBin  >=0 )
1295       {
1296          fhSumPtConeBinMC [ ptsumBinMC]->Fill(pt);
1297         if(fFillSSHisto)  fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,m02);
1298       }
1299
1300       
1301       if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1302       {
1303         leadptBinMC = leadptBin+kmcPhoton*fNBkgBin;
1304         ptsumBinMC  =  ptsumBin+kmcPhoton*fNBkgBin;
1305         if( leadptBin >=0 )
1306         {
1307           fhPtLeadConeBinMC[leadptBinMC]->Fill(pt);
1308           if(fFillSSHisto) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,m02);
1309         }
1310         
1311         if( ptsumBin  >=0 )
1312         {
1313           fhSumPtConeBinMC [ ptsumBinMC]->Fill(pt);
1314           if(fFillSSHisto)  fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,m02);
1315         }
1316       }
1317     }
1318   }
1319   
1320   if(fFillPtTrigBinHistograms)
1321   {
1322     // Get the background bin for this cone and trigger
1323     Int_t ptTrigBin  = -1;
1324     
1325     for(Int_t ibin = 0; ibin < fNPtTrigBin; ibin++)
1326     {
1327       if( pt  >= fPtTrigBinLimit[ibin] && coneptsum  < fPtTrigBinLimit[ibin+1]) ptTrigBin  = ibin;
1328     }
1329     
1330     // Fill the histograms per pT candidate bin of pt lead or pt sum
1331     if( GetDebug() > 1 && ptTrigBin  >=0 ) printf("Trigger pT %f, bin %d [%2.2f,%2.2f]\n" , pt , ptTrigBin, fPtTrigBinLimit[ptTrigBin] ,fPtTrigBinLimit[ptTrigBin +1]);
1332     
1333     if( ptTrigBin >=0 )
1334     {
1335       fhPtTrigBinPtLeadCone[ptTrigBin]->Fill(coneleadpt);
1336       fhPtTrigBinSumPtCone [ptTrigBin]->Fill(coneptsum );
1337       if(fFillSSHisto)
1338       {
1339         fhPtTrigBinLambda0vsPtLeadCone[ptTrigBin]->Fill(coneleadpt,m02);
1340         fhPtTrigBinLambda0vsSumPtCone [ptTrigBin]->Fill(coneptsum ,m02);
1341       }
1342     }
1343     
1344     // Check if it was a decay
1345     if( fFillTaggedDecayHistograms && m02 < 0.3 )
1346     {
1347       Int_t decayTag = pCandidate->DecayTag();
1348       if(decayTag < 0) decayTag = 0;
1349       for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1350       {
1351         if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
1352         {
1353           Int_t binDecay = ptTrigBin+ibit*fNPtTrigBin;
1354           if( binDecay > 0 )
1355           {
1356             fhPtTrigBinPtLeadConeDecay[binDecay]->Fill(coneleadpt);
1357             fhPtTrigBinSumPtConeDecay [binDecay]->Fill(coneptsum );
1358           }
1359         }
1360       }
1361     }
1362     
1363     if( IsDataMC() )
1364     {
1365       Int_t ptTrigBinMC = ptTrigBin+mcIndex*fNPtTrigBin;
1366       
1367       if( ptTrigBin >=0 )
1368       {
1369         fhPtTrigBinPtLeadConeMC[ptTrigBinMC]->Fill(coneleadpt);
1370         fhPtTrigBinSumPtConeMC [ptTrigBinMC]->Fill(coneptsum );
1371         if(fFillSSHisto)
1372         {
1373           fhPtTrigBinLambda0vsPtLeadConeMC[ptTrigBinMC]->Fill(coneleadpt,m02);
1374           fhPtTrigBinLambda0vsSumPtConeMC [ptTrigBinMC]->Fill(coneptsum ,m02);
1375         }
1376       }
1377       
1378       if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1379       {
1380         ptTrigBinMC = ptTrigBin+kmcPhoton*fNPtTrigBin;
1381         if( ptTrigBin >=0 )
1382         {
1383           fhPtTrigBinPtLeadConeMC[ptTrigBinMC]->Fill(coneleadpt);
1384           fhPtTrigBinSumPtConeMC [ptTrigBinMC]->Fill(coneptsum );
1385           if(fFillSSHisto)
1386           {
1387             fhPtTrigBinLambda0vsPtLeadConeMC[ptTrigBinMC]->Fill(coneleadpt,m02);
1388             fhPtTrigBinLambda0vsSumPtConeMC [ptTrigBinMC]->Fill(coneptsum ,m02);
1389           }
1390         }
1391       } // photon MC
1392     } // MC
1393   } // pT trigger bins
1394   
1395   // Shower shape dependent histograms
1396   if(fFillSSHisto)
1397   {
1398     fhELambda0 [isolated]->Fill(energy, m02);
1399     fhPtLambda0[isolated]->Fill(pt,     m02);
1400     //fhELambda1 [isolated]->Fill(energy, m20);
1401     
1402     if(IsDataMC())
1403     {
1404       if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1405         fhPtLambda0MC[kmcPhoton][isolated]->Fill(pt,m02);
1406       
1407       if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
1408         fhPtLambda0MC[kmcPi0DecayLostPair][isolated]->Fill(pt,m02);
1409       
1410       fhPtLambda0MC[mcIndex][isolated]->Fill(pt,m02);
1411     }
1412     
1413     if(GetCalorimeter() == kEMCAL &&  GetFirstSMCoveredByTRD() >= 0 &&
1414        GetModuleNumber(pCandidate) >= GetFirstSMCoveredByTRD()  )
1415     {
1416       fhELambda0TRD [isolated]->Fill(energy, m02 );
1417       fhPtLambda0TRD[isolated]->Fill(pt    , m02 );
1418       //fhELambda1TRD [isolated]->Fill(energy, m20 );
1419     }
1420     
1421     if(fFillNLMHistograms)
1422     {
1423       fhNLocMax[isolated]->Fill(energy,nMaxima);
1424       if     (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,m02); fhELambda1LocMax1[isolated]->Fill(energy,m02); }
1425       else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,m02); fhELambda1LocMax2[isolated]->Fill(energy,m02); }
1426       else                { fhELambda0LocMaxN[isolated]->Fill(energy,m02); fhELambda1LocMaxN[isolated]->Fill(energy,m02); }
1427     }
1428   } // SS histo fill
1429   
1430   // Track matching dependent histograms
1431   if(fFillTMHisto)
1432   {
1433     Int_t iclus = -1;
1434     TObjArray* clusters = 0x0;
1435     if     (GetCalorimeter() == kEMCAL) clusters = GetEMCALClusters();
1436     else if(GetCalorimeter() == kPHOS ) clusters = GetPHOSClusters();
1437     
1438     if(!clusters) return;
1439     
1440     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
1441     
1442     Float_t dZ  = cluster->GetTrackDz();
1443     Float_t dR  = cluster->GetTrackDx();
1444     
1445     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1446     {
1447       dR = 2000., dZ = 2000.;
1448       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1449     }
1450     
1451     //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
1452     if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
1453     {
1454       fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
1455       fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
1456       if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
1457     }
1458     
1459     // Check dEdx and E/p of matched clusters
1460     
1461     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1462     {
1463       
1464       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1465       
1466       if(track)
1467       {
1468         Float_t dEdx = track->GetTPCsignal();
1469         fhdEdx[isolated]->Fill(cluster->E(), dEdx);
1470         
1471         Float_t eOverp = cluster->E()/track->P();
1472         fhEOverP[isolated]->Fill(cluster->E(),  eOverp);
1473       }
1474       //else
1475       //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1476       
1477       
1478       if(IsDataMC())
1479       {
1480         if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
1481         {
1482           if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
1483                      GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
1484           else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
1485           else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
1486           else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
1487           
1488         }
1489         else
1490         {
1491           if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
1492                      GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
1493           else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
1494           else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
1495           else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
1496         }
1497       }  // MC
1498       
1499     } // match window
1500     
1501   }// TM histos fill
1502   
1503 }
1504
1505 //______________________________________________________
1506 TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
1507 {
1508   //Save parameters used for analysis
1509   TString parList ; //this will be list of parameters used for this analysis.
1510   const Int_t buffersize = 255;
1511   char onePar[buffersize] ;
1512   
1513   snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
1514   parList+=onePar ;
1515   snprintf(onePar, buffersize,"Calorimeter: %s\n",GetCalorimeterString().Data()) ;
1516   parList+=onePar ;
1517   snprintf(onePar, buffersize,"Isolation Cand Detector: %s\n",fIsoDetectorString.Data()) ;
1518   parList+=onePar ;
1519   snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
1520   parList+=onePar ;
1521   snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
1522   parList+=onePar ;
1523   snprintf(onePar, buffersize,"fFillTMHisto=%d (Flag for track matching histograms) \n",fFillTMHisto) ;
1524   parList+=onePar ;
1525   snprintf(onePar, buffersize,"fFillSSHisto=%d (Flag for shower shape histograms) \n",fFillSSHisto) ;
1526   parList+=onePar ;
1527   
1528   if(fMakeSeveralIC)
1529   {
1530     snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
1531     parList+=onePar ;
1532     snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
1533     parList+=onePar ;
1534     
1535     for(Int_t icone = 0; icone < fNCones ; icone++)
1536     {
1537       snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
1538       parList+=onePar ;
1539     }
1540     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1541     {
1542       snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
1543       parList+=onePar ;
1544     }
1545     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1546     {
1547       snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
1548       parList+=onePar ;
1549     }
1550     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1551     {
1552       snprintf(onePar, buffersize,"fSumPtThresholds[%d]=%1.2f (isolation sum pt threshold) \n",ipt, fSumPtThresholds[ipt]) ;
1553       parList+=onePar ;
1554     }
1555   }
1556   
1557   //Get parameters set in base class.
1558   parList += GetBaseParametersList() ;
1559   
1560   //Get parameters set in IC class.
1561   if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
1562   
1563   return new TObjString(parList) ;
1564   
1565 }
1566
1567 //________________________________________________________
1568 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
1569 {
1570   // Create histograms to be saved in output file and
1571   // store them in outputContainer
1572   TList * outputContainer = new TList() ;
1573   outputContainer->SetName("IsolatedParticleHistos") ;
1574   
1575   Int_t   nptbins  = GetHistogramRanges()->GetHistoPtBins();
1576   Int_t   nphibins = GetHistogramRanges()->GetHistoPhiBins();
1577   Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
1578   Float_t ptmax    = GetHistogramRanges()->GetHistoPtMax();
1579   Float_t phimax   = GetHistogramRanges()->GetHistoPhiMax();
1580   Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
1581   Float_t ptmin    = GetHistogramRanges()->GetHistoPtMin();
1582   Float_t phimin   = GetHistogramRanges()->GetHistoPhiMin();
1583   Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
1584   Int_t   ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();
1585   Float_t ssmax    = GetHistogramRanges()->GetHistoShowerShapeMax();
1586   Float_t ssmin    = GetHistogramRanges()->GetHistoShowerShapeMin();
1587   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
1588   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
1589   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
1590   
1591   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1592   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1593   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1594   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1595   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1596   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1597   
1598   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
1599   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
1600   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1601   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1602   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
1603   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1604   
1605   Int_t   nptsumbins    = GetHistogramRanges()->GetHistoNPtSumBins();
1606   Float_t ptsummax      = GetHistogramRanges()->GetHistoPtSumMax();
1607   Float_t ptsummin      = GetHistogramRanges()->GetHistoPtSumMin();
1608   Int_t   nptinconebins = GetHistogramRanges()->GetHistoNPtInConeBins();
1609   Float_t ptinconemax   = GetHistogramRanges()->GetHistoPtInConeMax();
1610   Float_t ptinconemin   = GetHistogramRanges()->GetHistoPtInConeMin();
1611   
1612   //Float_t ptthre    = GetIsolationCut()->GetPtThreshold();
1613   //Float_t ptsumthre = GetIsolationCut()->GetSumPtThreshold();
1614   //Float_t ptfrac    = GetIsolationCut()->GetPtFraction();
1615   Float_t r         = GetIsolationCut()->GetConeSize();
1616   Int_t   method    = GetIsolationCut()->GetICMethod() ;
1617   Int_t   particle  = GetIsolationCut()->GetParticleTypeInCone() ;
1618   
1619   TString sThreshold = "";
1620   if      ( method == AliIsolationCut::kSumPtIC )
1621   {
1622     sThreshold = Form(", %2.2f < #Sigma #it{p}_{T}^{in cone} < %2.2f GeV/#it{c}",
1623                       GetIsolationCut()->GetSumPtThreshold(), GetIsolationCut()->GetSumPtThresholdMax());
1624     if(GetIsolationCut()->GetSumPtThresholdMax() > 200)
1625       sThreshold = Form(", #Sigma #it{p}_{T}^{in cone} = %2.2f GeV/#it{c}",
1626                         GetIsolationCut()->GetSumPtThreshold());
1627   }
1628   else if ( method == AliIsolationCut::kPtThresIC)
1629   {
1630     sThreshold = Form(", %2.2f < #it{p}_{T}^{th} < %2.2f GeV/#it{c}",
1631                       GetIsolationCut()->GetPtThreshold(),GetIsolationCut()->GetPtThresholdMax());
1632     if(GetIsolationCut()->GetSumPtThreshold() > 200)
1633       sThreshold = Form(", #it{p}_{T}^{th} = %2.2f GeV/#it{c}",
1634                         GetIsolationCut()->GetPtThreshold());
1635   }
1636   else if ( method == AliIsolationCut::kPtFracIC)
1637     sThreshold = Form(", #Sigma #it{p}_{T}^{in cone}/#it{p}_{T}^{trig} = %2.2f" ,
1638                       GetIsolationCut()->GetPtFraction());
1639   
1640   TString sParticle = ", x^{0,#pm}";
1641   if      ( particle == AliIsolationCut::kOnlyNeutral )  sParticle = ", x^{0}";
1642   else if ( particle == AliIsolationCut::kOnlyCharged )  sParticle = ", x^{#pm}";
1643   
1644   TString parTitle = Form("#it{R} = %2.2f%s%s",GetIsolationCut()->GetConeSize(), sThreshold.Data(),sParticle.Data());
1645   
1646   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1647   
1648   // MC histograms title and name
1649   TString mcPartType[] = { "#gamma", "#gamma_{prompt}", "#gamma_{fragmentation}",
1650     "#pi^{0} (merged #gamma)","#gamma_{#pi decay}","#gamma_{#pi decay} lost companion",
1651     "#gamma_{#eta decay}","#gamma_{other decay}",
1652     "e^{#pm}","hadrons?"} ;
1653   
1654   TString mcPartName[] = { "Photon","PhotonPrompt","PhotonFrag",
1655     "Pi0","Pi0Decay","Pi0DecayLostPair","EtaDecay","OtherDecay",
1656     "Electron","Hadron"} ;
1657   
1658   // Primary MC histograms title and name
1659   TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}",
1660     "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","#pi^{0}"} ;
1661   
1662   TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay",
1663     "PhotonPrompt","PhotonFrag","PhotonISR","Pi0"} ;
1664   
1665   // Not Isolated histograms, reference histograms
1666   
1667   fhENoIso  = new TH1F("hENoIso",
1668                        Form("Number of not isolated leading particles vs #it{p}_{T}, %s",parTitle.Data()),
1669                        nptbins,ptmin,ptmax);
1670   fhENoIso->SetYTitle("#it{counts}");
1671   fhENoIso->SetXTitle("E (GeV/#it{c})");
1672   outputContainer->Add(fhENoIso) ;
1673   
1674   fhPtNoIso  = new TH1F("hPtNoIso",
1675                         Form("Number of not isolated leading particles vs #it{p}_{T}, %s",parTitle.Data()),
1676                         nptbins,ptmin,ptmax);
1677   fhPtNoIso->SetYTitle("#it{counts}");
1678   fhPtNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1679   outputContainer->Add(fhPtNoIso) ;
1680   
1681   fhEtaPhiNoIso  = new TH2F("hEtaPhiNoIso",
1682                             Form("Number of not isolated leading particles #eta vs #phi, %s",parTitle.Data()),
1683                             netabins,etamin,etamax,nphibins,phimin,phimax);
1684   fhEtaPhiNoIso->SetXTitle("#eta");
1685   fhEtaPhiNoIso->SetYTitle("#phi");
1686   outputContainer->Add(fhEtaPhiNoIso) ;
1687   
1688   if(IsDataMC())
1689   {
1690     // For histograms in arrays, index in the array, corresponding to any particle origin
1691     
1692     for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1693     {
1694       fhPtNoIsoMC[imc]  = new TH1F(Form("hPtNoIsoMC%s",mcPartName[imc].Data()),
1695                                    Form("#it{p}_{T} of NOT isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1696                                    nptbins,ptmin,ptmax);
1697       fhPtNoIsoMC[imc]->SetYTitle("#it{counts}");
1698       fhPtNoIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1699       outputContainer->Add(fhPtNoIsoMC[imc]) ;
1700       
1701       fhPtIsoMC[imc]  = new TH1F(Form("hPtMC%s",mcPartName[imc].Data()),
1702                                  Form("#it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1703                                  nptbins,ptmin,ptmax);
1704       fhPtIsoMC[imc]->SetYTitle("#it{counts}");
1705       fhPtIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1706       outputContainer->Add(fhPtIsoMC[imc]) ;
1707       
1708       fhPhiIsoMC[imc]  = new TH2F(Form("hPhiMC%s",mcPartName[imc].Data()),
1709                                   Form("#phi vs #it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1710                                   nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1711       fhPhiIsoMC[imc]->SetYTitle("#phi");
1712       fhPhiIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1713       outputContainer->Add(fhPhiIsoMC[imc]) ;
1714       
1715       fhEtaIsoMC[imc]  = new TH2F(Form("hEtaMC%s",mcPartName[imc].Data()),
1716                                   Form("#phi vs #it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1717                                   nptbins,ptmin,ptmax,netabins,etamin,etamax);
1718       fhEtaIsoMC[imc]->SetYTitle("#eta");
1719       fhEtaIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1720       outputContainer->Add(fhEtaIsoMC[imc]) ;
1721     }
1722   }
1723   
1724   // Histograms for tagged candidates as decay
1725   if(fFillTaggedDecayHistograms)
1726   {
1727     TString isoName [] = {"NoIso","Iso"};
1728     TString isoTitle[] = {"Not isolated"  ,"isolated"};
1729     
1730     for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1731     {
1732       for(Int_t iso = 0; iso < 2; iso++)
1733       {
1734         if(fMakeSeveralIC && iso) continue;
1735         fhPtDecay[iso][ibit]  =
1736         new TH1F(Form("hPtDecay%s_bit%d",isoName[iso].Data(),fDecayBits[ibit]),
1737                  Form("Number of %s leading pi0 decay particles vs #it{p}_{T}, bit %d, %s",isoTitle[iso].Data(),fDecayBits[ibit],parTitle.Data()),
1738                  nptbins,ptmin,ptmax);
1739         fhPtDecay[iso][ibit]->SetYTitle("#it{counts}");
1740         fhPtDecay[iso][ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1741         outputContainer->Add(fhPtDecay[iso][ibit]) ;
1742         
1743         fhEtaPhiDecay[iso][ibit]  =
1744         new TH2F(Form("hEtaPhiDecay%s_bit%d",isoName[iso].Data(),fDecayBits[ibit]),
1745                  Form("Number of %s leading Pi0 decay particles #eta vs #phi, bit %d, %s",isoTitle[iso].Data(),fDecayBits[ibit],parTitle.Data()),
1746                  netabins,etamin,etamax,nphibins,phimin,phimax);
1747         fhEtaPhiDecay[iso][ibit]->SetXTitle("#eta");
1748         fhEtaPhiDecay[iso][ibit]->SetYTitle("#phi");
1749         outputContainer->Add(fhEtaPhiDecay[iso][ibit]) ;
1750         
1751         if(fFillSSHisto)
1752         {
1753           fhPtLambda0Decay[iso][ibit]  = new TH2F
1754           (Form("hPtLambda0Decay%s_bit%d",isoName[iso].Data(),fDecayBits[ibit]),
1755            Form("%s cluster : #it{p}_{T} vs #lambda_{0}, decay bit %d, %s",isoTitle[iso].Data(), fDecayBits[ibit], parTitle.Data()),
1756            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1757           fhPtLambda0Decay[iso][ibit]->SetYTitle("#lambda_{0}^{2}");
1758           fhPtLambda0Decay[iso][ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1759           outputContainer->Add(fhPtLambda0Decay[iso][ibit]) ;
1760         }
1761         
1762         if(IsDataMC())
1763         {
1764           for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1765           {
1766             fhPtDecayMC[iso][ibit][imc]  =
1767             new TH1F(Form("hPtDecay%s_bit%d_MC%s",isoName[iso].Data(),fDecayBits[ibit],mcPartName[imc].Data()),
1768                      Form("#it{p}_{T} of %s, decay bit %d,  %s, %s",isoTitle[iso].Data(),fDecayBits[ibit],mcPartType[imc].Data(),parTitle.Data()),
1769                      nptbins,ptmin,ptmax);
1770             fhPtDecayMC[iso][ibit][imc]->SetYTitle("#it{counts}");
1771             fhPtDecayMC[iso][ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1772             outputContainer->Add(fhPtDecayMC[iso][ibit][imc]) ;
1773           }// MC particle loop
1774         }// MC
1775       } // bit loop
1776     } //iso loop
1777   }// decay
1778   
1779   if(!fMakeSeveralIC)
1780   {
1781     TString isoName [] = {"NoIso","Iso"};
1782     TString isoTitle[] = {"Not isolated"  ,"isolated"};
1783     
1784     fhEIso   = new TH1F("hE",
1785                         Form("Number of isolated particles vs E, %s",parTitle.Data()),
1786                         nptbins,ptmin,ptmax);
1787     fhEIso->SetYTitle("d#it{N} / d#it{E}");
1788     fhEIso->SetXTitle("#it{E} (GeV/#it{c})");
1789     outputContainer->Add(fhEIso) ;
1790     
1791     fhPtIso  = new TH1F("hPt",
1792                         Form("Number of isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1793                         nptbins,ptmin,ptmax);
1794     fhPtIso->SetYTitle("d#it{N} / #it{p}_{T}");
1795     fhPtIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1796     outputContainer->Add(fhPtIso) ;
1797     
1798     fhPhiIso  = new TH2F("hPhi",
1799                          Form("Number of isolated particles vs #phi, %s",parTitle.Data()),
1800                          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1801     fhPhiIso->SetYTitle("#phi");
1802     fhPhiIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1803     outputContainer->Add(fhPhiIso) ;
1804     
1805     fhEtaIso  = new TH2F("hEta",
1806                          Form("Number of isolated particles vs #eta, %s",parTitle.Data()),
1807                          nptbins,ptmin,ptmax,netabins,etamin,etamax);
1808     fhEtaIso->SetYTitle("#eta");
1809     fhEtaIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1810     outputContainer->Add(fhEtaIso) ;
1811     
1812     fhEtaPhiIso  = new TH2F("hEtaPhiIso",
1813                             Form("Number of isolated particles #eta vs #phi, %s",parTitle.Data()),
1814                             netabins,etamin,etamax,nphibins,phimin,phimax);
1815     fhEtaPhiIso->SetXTitle("#eta");
1816     fhEtaPhiIso->SetYTitle("#phi");
1817     outputContainer->Add(fhEtaPhiIso) ;
1818     
1819     if(IsHighMultiplicityAnalysisOn())
1820     {
1821       fhPtCentralityIso  = new TH2F("hPtCentrality",
1822                                     Form("centrality vs #it{p}_{T} for isolated particles, %s",parTitle.Data()),
1823                                     nptbins,ptmin,ptmax, 100,0,100);
1824       fhPtCentralityIso->SetYTitle("centrality");
1825       fhPtCentralityIso->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1826       outputContainer->Add(fhPtCentralityIso) ;
1827       
1828       fhPtEventPlaneIso  = new TH2F("hPtEventPlane",
1829                                     Form("event plane angle vs #it{p}_{T} for isolated particles, %s",parTitle.Data()),
1830                                     nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1831       fhPtEventPlaneIso->SetYTitle("Event plane angle (rad)");
1832       fhPtEventPlaneIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1833       outputContainer->Add(fhPtEventPlaneIso) ;
1834     }
1835     
1836     if(fFillNLMHistograms)
1837     {
1838       fhPtNLocMaxIso  = new TH2F("hPtNLocMax",
1839                                  Form("Number of isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1840                                  nptbins,ptmin,ptmax,10,0,10);
1841       fhPtNLocMaxIso->SetYTitle("#it{NLM}");
1842       fhPtNLocMaxIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1843       
1844       fhPtNLocMaxNoIso  = new TH2F("hPtNLocMaxNoIso",
1845                                    Form("Number of not isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1846                                    nptbins,ptmin,ptmax,10,0,10);
1847       fhPtNLocMaxNoIso->SetYTitle("#it{NLM}");
1848       fhPtNLocMaxNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1849       outputContainer->Add(fhPtNLocMaxNoIso) ;
1850     }
1851
1852     fhConePtLead  = new TH2F("hConePtLead",
1853                             Form("Track or Cluster  leading #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1854                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1855     fhConePtLead->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
1856     fhConePtLead->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1857     outputContainer->Add(fhConePtLead) ;
1858     
1859     fhConeSumPt  = new TH2F("hConePtSum",
1860                             Form("Track and Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1861                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1862     fhConeSumPt->SetYTitle("#Sigma #it{p}_{T}");
1863     fhConeSumPt->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1864     outputContainer->Add(fhConeSumPt) ;
1865     
1866     fhConeSumPtTrigEtaPhi  = new TH2F("hConePtSumTrigEtaPhi",
1867                                       Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1868                                       netabins,etamin,etamax,nphibins,phimin,phimax);
1869     fhConeSumPtTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1870     fhConeSumPtTrigEtaPhi->SetXTitle("#eta_{trigger}");
1871     fhConeSumPtTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1872     outputContainer->Add(fhConeSumPtTrigEtaPhi) ;
1873     
1874     fhPtInCone  = new TH2F("hPtInCone",
1875                            Form("#it{p}_{T} of clusters and tracks in isolation cone for #it{R} =  %2.2f",r),
1876                            nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1877     fhPtInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1878     fhPtInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1879     outputContainer->Add(fhPtInCone) ;
1880     
1881     if(fFillBackgroundBinHistograms)
1882     {
1883       fhPtLeadConeBin              = new TH1F*[fNBkgBin];
1884       fhSumPtConeBin               = new TH1F*[fNBkgBin];
1885       if(fFillSSHisto)
1886       {
1887         fhPtLeadConeBinLambda0     = new TH2F*[fNBkgBin];
1888         fhSumPtConeBinLambda0      = new TH2F*[fNBkgBin];
1889       }
1890       
1891       if(fFillTaggedDecayHistograms)
1892       {
1893         fhPtLeadConeBinDecay       = new TH1F*[fNBkgBin*fNDecayBits];
1894         fhSumPtConeBinDecay        = new TH1F*[fNBkgBin*fNDecayBits];
1895       }
1896       
1897       if(IsDataMC())
1898       {
1899         fhPtLeadConeBinMC          = new TH1F*[fNBkgBin*fgkNmcTypes];
1900         fhSumPtConeBinMC           = new TH1F*[fNBkgBin*fgkNmcTypes];
1901         
1902         if(fFillSSHisto)
1903         {
1904           fhPtLeadConeBinLambda0MC = new TH2F*[fNBkgBin*fgkNmcTypes];
1905           fhSumPtConeBinLambda0MC  = new TH2F*[fNBkgBin*fgkNmcTypes];
1906         }
1907       }
1908       
1909       for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1910       {
1911         fhPtLeadConeBin[ibin]  = new TH1F
1912         (Form("hPtLeadCone_Bin%d",ibin),
1913          Form("cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1914               fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1915         fhPtLeadConeBin[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
1916         fhPtLeadConeBin[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1917         outputContainer->Add(fhPtLeadConeBin[ibin]) ;
1918         
1919         fhSumPtConeBin[ibin]  = new TH1F
1920         (Form("hSumPtCone_Bin%d",ibin),
1921          Form("in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, %s",
1922               fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1923         fhSumPtConeBin[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
1924         fhSumPtConeBin[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1925         outputContainer->Add(fhSumPtConeBin[ibin]) ;
1926         
1927         if(fFillTaggedDecayHistograms)
1928         {
1929           for(Int_t idecay = 0; idecay < fNDecayBits; idecay++)
1930           {
1931             Int_t bindecay = ibin+idecay*fNBkgBin;
1932
1933             fhPtLeadConeBinDecay[bindecay]  = new TH1F
1934             (Form("hPtLeadCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
1935              Form("Decay bit %d, cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1936                   fDecayBits[idecay],fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1937             fhPtLeadConeBinDecay[bindecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
1938             fhPtLeadConeBinDecay[bindecay]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1939             outputContainer->Add(fhPtLeadConeBinDecay[bindecay]) ;
1940             
1941             fhSumPtConeBinDecay[bindecay]  = new TH1F
1942             (Form("hSumPtCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
1943              Form("Decay bit %d, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c},  %s",
1944                   fDecayBits[idecay],fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
1945             fhSumPtConeBinDecay[bindecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
1946             fhSumPtConeBinDecay[bindecay]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1947             outputContainer->Add(fhSumPtConeBinDecay[bindecay]) ;
1948           }
1949         }
1950         
1951         if(IsDataMC())
1952         {
1953           for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1954           {
1955             Int_t binmc = ibin+imc*fNBkgBin;
1956             fhPtLeadConeBinMC[binmc]  = new TH1F
1957             (Form("hPtLeadCone_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
1958              Form("in cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, MC %s, %s",
1959                   fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
1960             fhPtLeadConeBinMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
1961             fhPtLeadConeBinMC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1962             outputContainer->Add(fhPtLeadConeBinMC[binmc]) ;
1963             
1964             fhSumPtConeBinMC[binmc]  = new TH1F
1965             (Form("hSumPtCone_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
1966              Form("in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, MC %s, %s",
1967                   fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
1968             fhSumPtConeBinMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
1969             fhSumPtConeBinMC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1970             outputContainer->Add(fhSumPtConeBinMC[binmc]) ;
1971           } // MC particle loop
1972         }
1973         
1974         if(fFillSSHisto)
1975         {
1976           fhPtLeadConeBinLambda0[ibin]  = new TH2F
1977           (Form("hPtLeadConeLambda0_Bin%d",ibin),
1978            Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, %s",
1979                 fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1980           fhPtLeadConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
1981           fhPtLeadConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1982           outputContainer->Add(fhPtLeadConeBinLambda0[ibin]) ;
1983           
1984           fhSumPtConeBinLambda0[ibin]  = new TH2F
1985           (Form("hSumPtConeLambda0_Bin%d",ibin),
1986            Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, %s",
1987                 fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1988           fhSumPtConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
1989           fhSumPtConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1990           outputContainer->Add(fhSumPtConeBinLambda0[ibin]) ;
1991           
1992           if(IsDataMC())
1993           {
1994             for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1995             {
1996               Int_t binmc = ibin+imc*fNBkgBin;
1997               fhPtLeadConeBinLambda0MC[binmc]  = new TH2F
1998               (Form("hPtLeadConeLambda0_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
1999                Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f GeV/#it{c}, MC %s, %s",
2000                     fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2001               fhPtLeadConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
2002               fhPtLeadConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2003               outputContainer->Add(fhPtLeadConeBinLambda0MC[binmc]) ;
2004               
2005               fhSumPtConeBinLambda0MC[binmc]  = new TH2F
2006               (Form("hSumPtConeLambda0_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
2007                Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f GeV/#it{c}, MC %s, %s",
2008                     fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2009               fhSumPtConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
2010               fhSumPtConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2011               outputContainer->Add(fhSumPtConeBinLambda0MC[binmc]) ;
2012             } // MC particle loop
2013           }
2014         } // shower shape on
2015       } // pt in cone bin loop
2016     } // bkg cone pt bin histograms
2017
2018     if(fFillPtTrigBinHistograms)
2019     {
2020       fhPtTrigBinPtLeadCone = new TH1F*[fNPtTrigBin];
2021       fhPtTrigBinSumPtCone  = new TH1F*[fNPtTrigBin];
2022       
2023       fhPtTrigBinPtLeadConeDecay = new TH1F*[fNPtTrigBin];
2024       fhPtTrigBinSumPtConeDecay  = new TH1F*[fNPtTrigBin];
2025       
2026       if(IsDataMC())
2027       {
2028         fhPtTrigBinPtLeadConeMC = new TH1F*[fNPtTrigBin*fgkNmcTypes];
2029         fhPtTrigBinSumPtConeMC  = new TH1F*[fNPtTrigBin*fgkNmcTypes];
2030       }
2031
2032       if(fFillSSHisto)
2033       {
2034         fhPtTrigBinLambda0vsPtLeadCone = new TH2F*[fNPtTrigBin];
2035         fhPtTrigBinLambda0vsSumPtCone  = new TH2F*[fNPtTrigBin];
2036         
2037         if(IsDataMC())
2038         {
2039           fhPtTrigBinLambda0vsPtLeadConeMC = new TH2F*[fNPtTrigBin*fgkNmcTypes];
2040           fhPtTrigBinLambda0vsSumPtConeMC  = new TH2F*[fNPtTrigBin*fgkNmcTypes];
2041         }
2042       }
2043       
2044       for(Int_t ibin = 0; ibin < fNPtTrigBin; ibin++)
2045       {
2046         fhPtTrigBinPtLeadCone[ibin]  = new TH1F
2047         (Form("hPtTrigBin_PtLeadCone_Bin%d",ibin),
2048          Form("#it{p}_{T}^{lead. in cone}, %2.2f<#it{p}_{T}^{cand}<%2.2f GeV/#it{c}, %s",
2049               fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
2050         fhPtTrigBinPtLeadCone[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
2051         fhPtTrigBinPtLeadCone[ibin]->SetXTitle("#it{p}_{T}^{in cone} (GeV/#it{c})");
2052         outputContainer->Add(fhPtTrigBinPtLeadCone[ibin]) ;
2053         
2054         fhPtTrigBinSumPtCone[ibin]  = new TH1F
2055         (Form("hPtTrigBin_SumPtCone_Bin%d",ibin),
2056          Form("#Sigma #it{p}_{T}^{in cone} %2.2f <#it{p}_{T}^{cand}< %2.2f GeV/#it{c}, %s",
2057               fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], parTitle.Data()),nptsumbins,ptsummin,ptsummax);
2058         fhPtTrigBinSumPtCone[ibin]->SetYTitle("d #it{N} / d #it{p}_{T}");
2059         fhPtTrigBinSumPtCone[ibin]->SetXTitle("#Sigma #it{p}_{T}^{in cone} (GeV/#it{c})");
2060         outputContainer->Add(fhPtTrigBinSumPtCone[ibin]) ;
2061
2062         if(fFillTaggedDecayHistograms)
2063         {
2064           for(Int_t idecay = 0; idecay < fNDecayBits; idecay++)
2065           {
2066             Int_t binDecay = ibin+idecay*fNPtTrigBin;
2067             
2068             fhPtTrigBinPtLeadConeDecay[binDecay]  = new TH1F
2069             (Form("hPtTrigBin_PtLeadCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
2070              Form("Decay bit %d, #it{p}_{T}^{lead. in cone}, %2.2f<#it{p}_{T}^{cand}<%2.2f GeV/#it{c}, %s",
2071                   fDecayBits[idecay],fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax);
2072             fhPtTrigBinPtLeadConeDecay[binDecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
2073             fhPtTrigBinPtLeadConeDecay[binDecay]->SetXTitle("#it{p}_{T}^{lead in cone} (GeV/#it{c})");
2074             outputContainer->Add(fhPtTrigBinPtLeadConeDecay[binDecay]) ;
2075             
2076             fhPtTrigBinSumPtConeDecay[binDecay]  = new TH1F
2077             (Form("hPtTrigBin_SumPtCone_Bin%d_DecayBit%d",ibin,fDecayBits[idecay]),
2078              Form("Decay bit %d, #Sigma #it{p}_{T}^{in cone} %2.2f <#it{p}_{T}^{cand}< %2.2f GeV/#it{c}, %s",
2079                   fDecayBits[idecay],fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], parTitle.Data()),nptsumbins,ptsummin,ptsummax);
2080             fhPtTrigBinSumPtConeDecay[binDecay]->SetYTitle("d #it{N} / d #it{p}_{T}");
2081             fhPtTrigBinSumPtConeDecay[binDecay]->SetXTitle("#Sigma #it{p}_{T}^{in cone} (GeV/#it{c})");
2082             outputContainer->Add(fhPtTrigBinSumPtConeDecay[binDecay]) ;
2083           }
2084         }
2085         
2086         if(IsDataMC())
2087         {
2088           for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
2089           {
2090             Int_t binmc = ibin+imc*fNPtTrigBin;
2091             fhPtTrigBinPtLeadConeMC[binmc]  = new TH1F
2092             (Form("hPtTrigBin_PtLeadCone_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
2093              Form("#it{p}_{T}^{lead. in cone}, %2.2f<#it{p}_{T}^{cand}<%2.2f GeV/#it{c}, MC %s, %s",
2094                   fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax);
2095             fhPtTrigBinPtLeadConeMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
2096             fhPtTrigBinPtLeadConeMC[binmc]->SetXTitle("#it{p}_{T}^{lead in cone} (GeV/#it{c})");
2097             outputContainer->Add(fhPtTrigBinPtLeadConeMC[binmc]) ;
2098             
2099             fhPtTrigBinSumPtConeMC[binmc]  = new TH1F
2100             (Form("hPtTrigBin_SumPtCone_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
2101              Form("#Sigma #it{p}_{T}^{in cone}, %2.2f <#it{p}_{T}^{cand}< %2.2f GeV/#it{c}, MC %s, %s",
2102                   fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptsumbins,ptsummin,ptsummax);
2103             fhPtTrigBinSumPtConeMC[binmc]->SetYTitle("d #it{N} / d #it{p}_{T}");
2104             fhPtTrigBinSumPtConeMC[binmc]->SetXTitle("#Sigma #it{p}_{T}^{in cone} (GeV/#it{c})");
2105             outputContainer->Add(fhPtTrigBinSumPtConeMC[binmc]) ;
2106           } // MC particle loop
2107         } // MC
2108
2109         if(fFillSSHisto)
2110         {
2111           fhPtTrigBinLambda0vsPtLeadCone[ibin]  = new TH2F
2112           (Form("hPtTrigBin_PtLeadConeVSLambda0_Bin%d",ibin),
2113            Form("#lambda_{0} vs #it{p}_{T}^{lead. in cone}, %2.2f<#it{p}_{T}^{cand}<%2.2f GeV/#it{c}, %s",
2114                 fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2115           fhPtTrigBinLambda0vsPtLeadCone[ibin]->SetYTitle("#lambda_{0}^{2}");
2116           fhPtTrigBinLambda0vsPtLeadCone[ibin]->SetXTitle("#it{p}_{T}^{lead in cone} (GeV/#it{c})");
2117           outputContainer->Add(fhPtTrigBinLambda0vsPtLeadCone[ibin]) ;
2118           
2119           fhPtTrigBinLambda0vsSumPtCone[ibin]  = new TH2F
2120           (Form("hPtTrigBin_SumPtConeVSLambda0_Bin%d",ibin),
2121            Form("#lambda_{0} vs #Sigma #it{p}_{T}^{in cone} %2.2f <#it{p}_{T}^{cand}< %2.2f GeV/#it{c}, %s",
2122                 fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], parTitle.Data()),nptsumbins,ptsummin,ptsummax,ssbins,ssmin,ssmax);
2123           fhPtTrigBinLambda0vsSumPtCone[ibin]->SetYTitle("#lambda_{0}^{2}");
2124           fhPtTrigBinLambda0vsSumPtCone[ibin]->SetXTitle("#Sigma #it{p}_{T}^{in cone} (GeV/#it{c})");
2125           outputContainer->Add(fhPtTrigBinLambda0vsSumPtCone[ibin]) ;
2126           
2127           if(IsDataMC())
2128           {
2129             for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
2130             {
2131               Int_t binmc = ibin+imc*fNPtTrigBin;
2132               fhPtTrigBinLambda0vsPtLeadConeMC[binmc]  = new TH2F
2133               (Form("hPtTrigBin_PtLeadConeVSLambda0_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
2134                Form("#lambda_{0} vs #it{p}_{T}^{lead. in cone}, %2.2f<#it{p}_{T}^{cand}<%2.2f GeV/#it{c}, MC %s, %s",
2135                     fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2136               fhPtTrigBinLambda0vsPtLeadConeMC[binmc]->SetYTitle("#lambda_{0}^{2}");
2137               fhPtTrigBinLambda0vsPtLeadConeMC[binmc]->SetXTitle("#it{p}_{T}^{lead in cone} (GeV/#it{c})");
2138               outputContainer->Add(fhPtTrigBinLambda0vsPtLeadConeMC[binmc]) ;
2139               
2140               fhPtTrigBinLambda0vsSumPtConeMC[binmc]  = new TH2F
2141               (Form("hPtTrigBin_SumPtConeVSLambda0_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
2142                Form("#lambda_{0} vs #Sigma #it{p}_{T}^{in cone}, %2.2f <#it{p}_{T}^{cand}< %2.2f GeV/#it{c}, MC %s, %s",
2143                     fPtTrigBinLimit[ibin],fPtTrigBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptsumbins,ptsummin,ptsummax,ssbins,ssmin,ssmax);
2144               fhPtTrigBinLambda0vsSumPtConeMC[binmc]->SetYTitle("#lambda_{0}^{2}");
2145               fhPtTrigBinLambda0vsSumPtConeMC[binmc]->SetXTitle("#Sigma #it{p}_{T}^{in cone} (GeV/#it{c})");
2146               outputContainer->Add(fhPtTrigBinLambda0vsSumPtConeMC[binmc]) ;
2147             } // MC particle loop
2148           } // MC
2149         } // SS histo
2150       } // pt trig bin loop
2151     } // pt trig bin histograms
2152     
2153     if(IsHighMultiplicityAnalysisOn())
2154     {
2155       fhPtInConeCent  = new TH2F("hPtInConeCent",
2156                                  Form("#it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2157                                  100,0,100,nptinconebins,ptinconemin,ptinconemax);
2158       fhPtInConeCent->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2159       fhPtInConeCent->SetXTitle("centrality");
2160       outputContainer->Add(fhPtInConeCent) ;
2161     }
2162     
2163     // Cluster only histograms
2164     if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyCharged)
2165     {
2166       fhConeSumPtCluster  = new TH2F("hConePtSumCluster",
2167                                      Form("Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2168                                      nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2169       fhConeSumPtCluster->SetYTitle("#Sigma #it{p}_{T}");
2170       fhConeSumPtCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2171       outputContainer->Add(fhConeSumPtCluster) ;
2172       
2173       fhConePtLeadCluster  = new TH2F("hConeLeadPtCluster",
2174                                     Form("Cluster leading in isolation cone for #it{R} =  %2.2f",r),
2175                                     nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2176       fhConePtLeadCluster->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
2177       fhConePtLeadCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2178       outputContainer->Add(fhConePtLeadCluster) ;
2179
2180       
2181       if(fFillCellHistograms)
2182       {
2183         fhConeSumPtCell  = new TH2F("hConePtSumCell",
2184                                     Form("Cell #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2185                                     nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2186         fhConeSumPtCell->SetYTitle("#Sigma #it{p}_{T}");
2187         fhConeSumPtCell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2188         outputContainer->Add(fhConeSumPtCell) ;
2189       }
2190       
2191       if(fFillUEBandSubtractHistograms)
2192       {
2193         fhConeSumPtEtaBandUECluster  = new TH2F("hConePtSumEtaBandUECluster",
2194                                                 "#Sigma cluster #it{p}_{T} in UE Eta Band",
2195                                                 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2196         fhConeSumPtEtaBandUECluster->SetYTitle("#Sigma #it{p}_{T}");
2197         fhConeSumPtEtaBandUECluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2198         outputContainer->Add(fhConeSumPtEtaBandUECluster) ;
2199         
2200         fhConeSumPtPhiBandUECluster  = new TH2F("hConePtSumPhiBandUECluster",
2201                                                 "#Sigma cluster #it{p}_{T} UE Phi Band",
2202                                                 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2203         fhConeSumPtPhiBandUECluster->SetYTitle("#Sigma #it{p}_{T}");
2204         fhConeSumPtPhiBandUECluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2205         outputContainer->Add(fhConeSumPtPhiBandUECluster) ;
2206         
2207         fhConeSumPtEtaBandUEClusterTrigEtaPhi  = new TH2F("hConePtSumEtaBandUEClusterTrigEtaPhi",
2208                                                           "Trigger #eta vs #phi, #Sigma cluster #it{p}_{T} in UE Eta Band",
2209                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
2210         fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2211         fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
2212         fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2213         outputContainer->Add(fhConeSumPtEtaBandUEClusterTrigEtaPhi) ;
2214         
2215         fhConeSumPtPhiBandUEClusterTrigEtaPhi  = new TH2F("hConePtSumPhiBandUEClusterTrigEtaPhi",
2216                                                           "Trigger #eta vs #phi, #Sigma cluster #it{p}_{T} UE Phi Band",
2217                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
2218         fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2219         fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
2220         fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2221         outputContainer->Add(fhConeSumPtPhiBandUEClusterTrigEtaPhi) ;
2222         if(fFillCellHistograms)
2223         {
2224           
2225           fhConeSumPtEtaBandUECell  = new TH2F("hConePtSumEtaBandUECell",
2226                                                "#Sigma cell #it{p}_{T} in UE Eta Band",
2227                                                nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2228           fhConeSumPtEtaBandUECell->SetYTitle("#Sigma #it{p}_{T}");
2229           fhConeSumPtEtaBandUECell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2230           outputContainer->Add(fhConeSumPtEtaBandUECell) ;
2231           
2232           fhConeSumPtPhiBandUECell  = new TH2F("hConePtSumPhiBandUECell",
2233                                                "#Sigma cell #it{p}_{T} UE Phi Band",
2234                                                nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2235           fhConeSumPtPhiBandUECell->SetYTitle("#Sigma #it{p}_{T}");
2236           fhConeSumPtPhiBandUECell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2237           outputContainer->Add(fhConeSumPtPhiBandUECell) ;
2238           
2239           fhConeSumPtEtaBandUECellTrigEtaPhi  = new TH2F("hConePtSumEtaBandUECellTrigEtaPhi",
2240                                                          "Trigger #eta vs #phi, #Sigma cell #it{p}_{T} in UE Eta Band",
2241                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2242           fhConeSumPtEtaBandUECellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2243           fhConeSumPtEtaBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2244           fhConeSumPtEtaBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2245           outputContainer->Add(fhConeSumPtEtaBandUECellTrigEtaPhi) ;
2246           
2247           fhConeSumPtPhiBandUECellTrigEtaPhi  = new TH2F("hConePtSumPhiBandUECellTrigEtaPhi",
2248                                                          "Trigger #eta vs #phi, #Sigma cell #it{p}_{T} UE Phi Band",
2249                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2250           fhConeSumPtPhiBandUECellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2251           fhConeSumPtPhiBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2252           fhConeSumPtPhiBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2253           outputContainer->Add(fhConeSumPtPhiBandUECellTrigEtaPhi) ;
2254         }
2255         
2256         fhEtaBandCluster  = new TH2F("hEtaBandCluster",
2257                                      Form("#eta vs #phi of clusters in #eta band isolation cone for #it{R} =  %2.2f",r),
2258                                      netabins,-1,1,nphibins,0,TMath::TwoPi());
2259         fhEtaBandCluster->SetXTitle("#eta");
2260         fhEtaBandCluster->SetYTitle("#phi");
2261         outputContainer->Add(fhEtaBandCluster) ;
2262         
2263         fhPhiBandCluster  = new TH2F("hPhiBandCluster",
2264                                      Form("#eta vs #phi of clusters in #phi band isolation cone for #it{R} =  %2.2f",r),
2265                                      netabins,-1,1,nphibins,0,TMath::TwoPi());
2266         fhPhiBandCluster->SetXTitle("#eta");
2267         fhPhiBandCluster->SetYTitle("#phi");
2268         outputContainer->Add(fhPhiBandCluster) ;
2269         
2270         fhEtaPhiInConeCluster= new TH2F("hEtaPhiInConeCluster",
2271                                         Form("#eta vs #phi of clusters in cone for #it{R} =  %2.2f",r),
2272                                         netabins,-1,1,nphibins,0,TMath::TwoPi());
2273         fhEtaPhiInConeCluster->SetXTitle("#eta");
2274         fhEtaPhiInConeCluster->SetYTitle("#phi");
2275         outputContainer->Add(fhEtaPhiInConeCluster) ;
2276         
2277         fhEtaPhiCluster= new TH2F("hEtaPhiCluster",
2278                                   Form("#eta vs #phi of all clusters"),
2279                                   netabins,-1,1,nphibins,0,TMath::TwoPi());
2280         fhEtaPhiCluster->SetXTitle("#eta");
2281         fhEtaPhiCluster->SetYTitle("#phi");
2282         outputContainer->Add(fhEtaPhiCluster) ;
2283         
2284       }
2285       
2286       fhPtClusterInCone  = new TH2F("hPtClusterInCone",
2287                                     Form("#it{p}_{T} of clusters in isolation cone for #it{R} =  %2.2f",r),
2288                                     nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2289       fhPtClusterInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2290       fhPtClusterInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2291       outputContainer->Add(fhPtClusterInCone) ;
2292       
2293       if(fFillCellHistograms)
2294       {
2295         fhPtCellInCone  = new TH2F("hPtCellInCone",
2296                                    Form("#it{p}_{T} of cells in isolation cone for #it{R} =  %2.2f",r),
2297                                    nptbins,ptmin,ptmax,1000,0,50);
2298         fhPtCellInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2299         fhPtCellInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2300         outputContainer->Add(fhPtCellInCone) ;
2301         
2302         fhEtaBandCell  = new TH2F("hEtaBandCell",
2303                                   Form("#col vs #row of cells in #eta band isolation cone for #it{R} =  %2.2f",r),
2304                                   96,0,95,128,0,127);
2305         fhEtaBandCell->SetXTitle("#col");
2306         fhEtaBandCell->SetYTitle("#row");
2307         outputContainer->Add(fhEtaBandCell) ;
2308         
2309         fhPhiBandCell  = new TH2F("hPhiBandCell",
2310                                   Form("#col vs #row of cells in #phi band isolation cone for #it{R} =  %2.2f",r),
2311                                   96,0,95,128,0,127);
2312         fhPhiBandCell->SetXTitle("#col");
2313         fhPhiBandCell->SetYTitle("#row");
2314         outputContainer->Add(fhPhiBandCell) ;
2315       }
2316       
2317       if(fFillUEBandSubtractHistograms)
2318       {
2319         fhConeSumPtEtaUESubCluster  = new TH2F("hConeSumPtEtaUESubCluster",
2320                                                Form("Clusters #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2321                                                nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2322         fhConeSumPtEtaUESubCluster->SetYTitle("#Sigma #it{p}_{T}");
2323         fhConeSumPtEtaUESubCluster->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2324         outputContainer->Add(fhConeSumPtEtaUESubCluster) ;
2325         
2326         fhConeSumPtPhiUESubCluster  = new TH2F("hConeSumPtPhiUESubCluster",
2327                                                Form("Clusters #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2328                                                nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2329         fhConeSumPtPhiUESubCluster->SetYTitle("#Sigma #it{p}_{T}");
2330         fhConeSumPtPhiUESubCluster->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2331         outputContainer->Add(fhConeSumPtPhiUESubCluster) ;
2332         
2333         fhConeSumPtEtaUESubClusterTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubClusterTrigEtaPhi",
2334                                                          Form("Trigger #eta vs #phi, Clusters #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2335                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2336         fhConeSumPtEtaUESubClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2337         fhConeSumPtEtaUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
2338         fhConeSumPtEtaUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2339         outputContainer->Add(fhConeSumPtEtaUESubClusterTrigEtaPhi) ;
2340         
2341         fhConeSumPtPhiUESubClusterTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubClusterTrigEtaPhi",
2342                                                          Form("Trigger #eta vs #phi, Clusters #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2343                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2344         fhConeSumPtPhiUESubClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2345         fhConeSumPtPhiUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
2346         fhConeSumPtPhiUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2347         outputContainer->Add(fhConeSumPtPhiUESubClusterTrigEtaPhi) ;
2348         
2349         if(fFillCellHistograms)
2350         {
2351           fhConeSumPtEtaUESubCell  = new TH2F("hConeSumPtEtaUESubCell",
2352                                               Form("Cells #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2353                                               nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2354           fhConeSumPtEtaUESubCell->SetYTitle("#Sigma #it{p}_{T}");
2355           fhConeSumPtEtaUESubCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2356           outputContainer->Add(fhConeSumPtEtaUESubCell) ;
2357           
2358           fhConeSumPtPhiUESubCell  = new TH2F("hConeSumPtPhiUESubCell",
2359                                               Form("Cells #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2360                                               nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2361           fhConeSumPtPhiUESubCell->SetYTitle("#Sigma #it{p}_{T}");
2362           fhConeSumPtPhiUESubCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2363           outputContainer->Add(fhConeSumPtPhiUESubCell) ;
2364           
2365           fhConeSumPtEtaUESubCellTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubCellTrigEtaPhi",
2366                                                         Form("Trigger #eta vs #phi, Cells #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2367                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2368           fhConeSumPtEtaUESubCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2369           fhConeSumPtEtaUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2370           fhConeSumPtEtaUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2371           outputContainer->Add(fhConeSumPtEtaUESubCellTrigEtaPhi) ;
2372           
2373           fhConeSumPtPhiUESubCellTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubCellTrigEtaPhi",
2374                                                         Form("Trigger #eta vs #phi, Cells #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2375                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2376           fhConeSumPtPhiUESubCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2377           fhConeSumPtPhiUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2378           fhConeSumPtPhiUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2379           outputContainer->Add(fhConeSumPtPhiUESubCellTrigEtaPhi) ;
2380         }
2381         
2382         fhFractionClusterOutConeEta  = new TH2F("hFractionClusterOutConeEta",
2383                                                 Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #eta acceptance",r),
2384                                                 nptbins,ptmin,ptmax,100,0,1);
2385         fhFractionClusterOutConeEta->SetYTitle("#it{fraction}");
2386         fhFractionClusterOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2387         outputContainer->Add(fhFractionClusterOutConeEta) ;
2388         
2389         fhFractionClusterOutConeEtaTrigEtaPhi  = new TH2F("hFractionClusterOutConeEtaTrigEtaPhi",
2390                                                           Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #eta acceptance, in trigger #eta-#phi ",r),
2391                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
2392         fhFractionClusterOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2393         fhFractionClusterOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2394         fhFractionClusterOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2395         outputContainer->Add(fhFractionClusterOutConeEtaTrigEtaPhi) ;
2396         
2397         fhFractionClusterOutConePhi  = new TH2F("hFractionClusterOutConePhi",
2398                                                 Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #phi acceptance",r),
2399                                                 nptbins,ptmin,ptmax,100,0,1);
2400         fhFractionClusterOutConePhi->SetYTitle("#it{fraction}");
2401         fhFractionClusterOutConePhi->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2402         outputContainer->Add(fhFractionClusterOutConePhi) ;
2403         
2404         fhFractionClusterOutConePhiTrigEtaPhi  = new TH2F("hFractionClusterOutConePhiTrigEtaPhi",
2405                                                           Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #phi acceptance, in trigger #eta-#phi ",r),
2406                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
2407         fhFractionClusterOutConePhiTrigEtaPhi->SetZTitle("#it{fraction}");
2408         fhFractionClusterOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
2409         fhFractionClusterOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2410         outputContainer->Add(fhFractionClusterOutConePhiTrigEtaPhi) ;
2411         
2412         fhConeSumPtSubvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCluster",
2413                                                           Form("#Sigma #it{p}_{T} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2414                                                           nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2415         fhConeSumPtSubvsConeSumPtTotPhiCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2416         fhConeSumPtSubvsConeSumPtTotPhiCluster->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2417         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCluster);
2418         
2419         fhConeSumPtSubNormvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCluster",
2420                                                               Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2421                                                               nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2422         fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2423         fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2424         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCluster);
2425         
2426         fhConeSumPtSubvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCluster",
2427                                                           Form("#Sigma #it{p}_{T} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2428                                                           nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2429         fhConeSumPtSubvsConeSumPtTotEtaCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2430         fhConeSumPtSubvsConeSumPtTotEtaCluster->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2431         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCluster);
2432         
2433         fhConeSumPtSubNormvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCluster",
2434                                                               Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2435                                                               nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2436         fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2437         fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2438         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCluster);
2439         
2440         fhConeSumPtVSUEClusterEtaBand  = new TH2F("hConeSumPtVSUEClusterEtaBand",
2441                                                   Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in eta band for cluster (before normalization), R=%2.2f",r),
2442                                                   nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
2443         fhConeSumPtVSUEClusterEtaBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2444         fhConeSumPtVSUEClusterEtaBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2445         outputContainer->Add(fhConeSumPtVSUEClusterEtaBand);
2446         
2447         fhConeSumPtVSUEClusterPhiBand  = new TH2F("hConeSumPtVSUEClusterPhiBand",
2448                                                   Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in phi band for cluster (before normalization), R=%2.2f",r),
2449                                                   nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
2450         fhConeSumPtVSUEClusterPhiBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2451         fhConeSumPtVSUEClusterPhiBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2452         outputContainer->Add(fhConeSumPtVSUEClusterPhiBand);
2453         
2454         if(fFillCellHistograms)
2455         {
2456           fhFractionCellOutConeEta  = new TH2F("hFractionCellOutConeEta",
2457                                                Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #eta acceptance",r),
2458                                                nptbins,ptmin,ptmax,100,0,1);
2459           fhFractionCellOutConeEta->SetYTitle("#it{fraction}");
2460           fhFractionCellOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2461           outputContainer->Add(fhFractionCellOutConeEta) ;
2462           
2463           fhFractionCellOutConeEtaTrigEtaPhi  = new TH2F("hFractionCellOutConeEtaTrigEtaPhi",
2464                                                          Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #eta acceptance, in trigger #eta-#phi ",r),
2465                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2466           fhFractionCellOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2467           fhFractionCellOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2468           fhFractionCellOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2469           outputContainer->Add(fhFractionCellOutConeEtaTrigEtaPhi) ;
2470           
2471           fhFractionCellOutConePhi  = new TH2F("hFractionCellOutConePhi",
2472                                                Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #phi acceptance",r),
2473                                                nptbins,ptmin,ptmax,100,0,1);
2474           fhFractionCellOutConePhi->SetYTitle("#it{fraction}");
2475           fhFractionCellOutConePhi->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2476           outputContainer->Add(fhFractionCellOutConePhi) ;
2477           
2478           fhFractionCellOutConePhiTrigEtaPhi  = new TH2F("hFractionCellOutConePhiTrigEtaPhi",
2479                                                          Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #phi acceptance, in trigger #eta-#phi ",r),
2480                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2481           fhFractionCellOutConePhiTrigEtaPhi->SetZTitle("#it{fraction}");
2482           fhFractionCellOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
2483           fhFractionCellOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2484           outputContainer->Add(fhFractionCellOutConePhiTrigEtaPhi) ;
2485           
2486           
2487           fhConeSumPtSubvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCell",
2488                                                          Form("#Sigma #it{p}_{T} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2489                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2490           fhConeSumPtSubvsConeSumPtTotPhiCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2491           fhConeSumPtSubvsConeSumPtTotPhiCell->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2492           outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCell);
2493           
2494           fhConeSumPtSubNormvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCell",
2495                                                              Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2496                                                              nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2497           fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2498           fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2499           outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCell);
2500           
2501           fhConeSumPtSubvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCell",
2502                                                          Form("#Sigma #it{p}_{T} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2503                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2504           fhConeSumPtSubvsConeSumPtTotEtaCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2505           fhConeSumPtSubvsConeSumPtTotEtaCell->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2506           outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCell);
2507           
2508           fhConeSumPtSubNormvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCell",
2509                                                              Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2510                                                              nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2511           fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2512           fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2513           outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCell);
2514         }
2515       }
2516     }
2517     
2518     // Track only histograms
2519     if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyNeutral)
2520     {
2521       fhConeSumPtTrack  = new TH2F("hConePtSumTrack",
2522                                    Form("Track #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2523                                    nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2524       fhConeSumPtTrack->SetYTitle("#Sigma #it{p}_{T}");
2525       fhConeSumPtTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2526       outputContainer->Add(fhConeSumPtTrack) ;
2527
2528       fhConePtLeadTrack  = new TH2F("hConeLeadPtTrack",
2529                                    Form("Track leading in isolation cone for #it{R} =  %2.2f",r),
2530                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2531       fhConePtLeadTrack->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
2532       fhConePtLeadTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2533       outputContainer->Add(fhConePtLeadTrack) ;
2534       
2535       fhPtTrackInCone  = new TH2F("hPtTrackInCone",
2536                                   Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f",r),
2537                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2538       fhPtTrackInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2539       fhPtTrackInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2540       outputContainer->Add(fhPtTrackInCone) ;
2541       
2542       
2543       if(fFillUEBandSubtractHistograms)
2544       {
2545         fhConeSumPtEtaBandUETrack  = new TH2F("hConePtSumEtaBandUETrack",
2546                                               "#Sigma track #it{p}_{T} in UE Eta Band",
2547                                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2548         fhConeSumPtEtaBandUETrack->SetYTitle("#Sigma #it{p}_{T}");
2549         fhConeSumPtEtaBandUETrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2550         outputContainer->Add(fhConeSumPtEtaBandUETrack) ;
2551         
2552         fhConeSumPtPhiBandUETrack  = new TH2F("hConePtSumPhiBandUETrack",
2553                                               "#Sigma track #it{p}_{T} in UE Phi Band",
2554                                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax*8);
2555         fhConeSumPtPhiBandUETrack->SetYTitle("#Sigma #it{p}_{T}");
2556         fhConeSumPtPhiBandUETrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2557         outputContainer->Add(fhConeSumPtPhiBandUETrack) ;
2558         
2559         
2560         fhConeSumPtEtaBandUETrackTrigEtaPhi  = new TH2F("hConePtSumEtaBandUETrackTrigEtaPhi",
2561                                                         "Trigger #eta vs #phi, #Sigma track #it{p}_{T} in UE Eta Band",
2562                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2563         fhConeSumPtEtaBandUETrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2564         fhConeSumPtEtaBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2565         fhConeSumPtEtaBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2566         outputContainer->Add(fhConeSumPtEtaBandUETrackTrigEtaPhi) ;
2567         
2568         fhConeSumPtPhiBandUETrackTrigEtaPhi  = new TH2F("hConePtSumPhiBandUETrackTrigEtaPhi",
2569                                                         "Trigger #eta vs #phi, #Sigma track #it{p}_{T} in UE Phi Band",
2570                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2571         fhConeSumPtPhiBandUETrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2572         fhConeSumPtPhiBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2573         fhConeSumPtPhiBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2574         outputContainer->Add(fhConeSumPtPhiBandUETrackTrigEtaPhi) ;
2575         
2576         fhEtaBandTrack  = new TH2F("hEtaBandTrack",
2577                                    Form("#eta vs #phi of tracks in #eta band isolation cone for #it{R} =  %2.2f",r),
2578                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
2579         fhEtaBandTrack->SetXTitle("#eta");
2580         fhEtaBandTrack->SetYTitle("#phi");
2581         outputContainer->Add(fhEtaBandTrack) ;
2582         
2583         fhPhiBandTrack  = new TH2F("hPhiBandTrack",
2584                                    Form("#eta vs #phi of tracks in #phi band isolation cone for #it{R} =  %2.2f",r),
2585                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
2586         fhPhiBandTrack->SetXTitle("#eta");
2587         fhPhiBandTrack->SetYTitle("#phi");
2588         outputContainer->Add(fhPhiBandTrack) ;
2589         
2590         fhConeSumPtEtaUESubTrack  = new TH2F("hConeSumPtEtaUESubTrack",
2591                                              Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2592                                              nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2593         fhConeSumPtEtaUESubTrack->SetYTitle("#Sigma #it{p}_{T}");
2594         fhConeSumPtEtaUESubTrack->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2595         outputContainer->Add(fhConeSumPtEtaUESubTrack) ;
2596         
2597         fhConeSumPtPhiUESubTrack  = new TH2F("hConeSumPtPhiUESubTrack",
2598                                              Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2599                                              nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2600         fhConeSumPtPhiUESubTrack->SetYTitle("#Sigma #it{p}_{T}");
2601         fhConeSumPtPhiUESubTrack->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2602         outputContainer->Add(fhConeSumPtPhiUESubTrack) ;
2603         
2604         fhConeSumPtEtaUESubTrackTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrackTrigEtaPhi",
2605                                                        Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2606                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
2607         fhConeSumPtEtaUESubTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2608         fhConeSumPtEtaUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2609         fhConeSumPtEtaUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2610         outputContainer->Add(fhConeSumPtEtaUESubTrackTrigEtaPhi) ;
2611         
2612         fhConeSumPtPhiUESubTrackTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrackTrigEtaPhi",
2613                                                        Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2614                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
2615         fhConeSumPtPhiUESubTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2616         fhConeSumPtPhiUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2617         fhConeSumPtPhiUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2618         outputContainer->Add(fhConeSumPtPhiUESubTrackTrigEtaPhi) ;
2619         
2620         fhFractionTrackOutConeEta  = new TH2F("hFractionTrackOutConeEta",
2621                                               Form("Fraction of the isolation cone #it{R} =  %2.2f, out of tracks #eta acceptance",r),
2622                                               nptbins,ptmin,ptmax,100,0,1);
2623         fhFractionTrackOutConeEta->SetYTitle("#it{fraction}");
2624         fhFractionTrackOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2625         outputContainer->Add(fhFractionTrackOutConeEta) ;
2626         
2627         fhFractionTrackOutConeEtaTrigEtaPhi  = new TH2F("hFractionTrackOutConeEtaTrigEtaPhi",
2628                                                         Form("Fraction of the isolation cone #it{R} =  %2.2f, out of tracks #eta acceptance, in trigger #eta-#phi ",r),
2629                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2630         fhFractionTrackOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2631         fhFractionTrackOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2632         fhFractionTrackOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2633         outputContainer->Add(fhFractionTrackOutConeEtaTrigEtaPhi) ;
2634         
2635         fhConeSumPtSubvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubvsConeSumPtTotPhiTrack",
2636                                                         Form("#Sigma #it{p}_{T} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2637                                                         nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2638         fhConeSumPtSubvsConeSumPtTotPhiTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2639         fhConeSumPtSubvsConeSumPtTotPhiTrack->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2640         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiTrack);
2641         
2642         fhConeSumPtSubNormvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiTrack",
2643                                                             Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2644                                                             nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2645         fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2646         fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2647         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiTrack);
2648         
2649         fhConeSumPtSubvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubvsConeSumPtTotEtaTrack",
2650                                                         Form("#Sigma #it{p}_{T} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2651                                                         nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2652         fhConeSumPtSubvsConeSumPtTotEtaTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2653         fhConeSumPtSubvsConeSumPtTotEtaTrack->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2654         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaTrack);
2655         
2656         fhConeSumPtSubNormvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaTrack",
2657                                                             Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2658                                                             nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2659         fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2660         fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2661         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaTrack);
2662         
2663         
2664         // UE in perpendicular cone
2665         fhPerpConeSumPt  = new TH2F("hPerpConePtSum",
2666                                     Form("#Sigma #it{p}_{T} in isolation cone at #pm 45 degree phi from trigger particle, #it{R} =  %2.2f",r),
2667                                     nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2668         fhPerpConeSumPt->SetYTitle("#Sigma #it{p}_{T}");
2669         fhPerpConeSumPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2670         outputContainer->Add(fhPerpConeSumPt) ;
2671         
2672         fhPtInPerpCone  = new TH2F("hPtInPerpCone",
2673                                    Form("#it{p}_{T} in isolation cone at #pm 45 degree phi from trigger particle, #it{R} =  %2.2f",r),
2674                                    nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2675         fhPtInPerpCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2676         fhPtInPerpCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2677         outputContainer->Add(fhPtInPerpCone) ;
2678         
2679         fhEtaPhiTrack= new TH2F("hEtaPhiTrack",
2680                                 Form("#eta vs #phi of all Tracks"),
2681                                 netabins,-1,1,nphibins,0,TMath::TwoPi());
2682         fhEtaPhiTrack->SetXTitle("#eta");
2683         fhEtaPhiTrack->SetYTitle("#phi");
2684         outputContainer->Add(fhEtaPhiTrack) ;
2685         
2686         fhEtaPhiInConeTrack= new TH2F("hEtaPhiInConeTrack",
2687                                       Form("#eta vs #phi of Tracks in cone for #it{R} =  %2.2f",r),
2688                                       netabins,-1,1,nphibins,0,TMath::TwoPi());
2689         fhEtaPhiInConeTrack->SetXTitle("#eta");
2690         fhEtaPhiInConeTrack->SetYTitle("#phi");
2691         outputContainer->Add(fhEtaPhiInConeTrack) ;
2692         
2693         fhConeSumPtVSUETracksEtaBand  = new TH2F("hConeSumPtVSUETracksEtaBand",
2694                                                  Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in eta band for tracks (before normalization), R=%2.2f",r),
2695                                                  nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
2696         fhConeSumPtVSUETracksEtaBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2697         fhConeSumPtVSUETracksEtaBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2698         outputContainer->Add(fhConeSumPtVSUETracksEtaBand);
2699         
2700         fhConeSumPtVSUETracksPhiBand  = new TH2F("hConeSumPtVSUETracksPhiBand",
2701                                                  Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in phi band for tracks (before normalization), R=%2.2f",r),
2702                                                  nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
2703         fhConeSumPtVSUETracksPhiBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2704         fhConeSumPtVSUETracksPhiBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2705         outputContainer->Add(fhConeSumPtVSUETracksPhiBand);
2706       }
2707     }
2708     
2709     if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged )
2710     {
2711       fhConeSumPtClustervsTrack   = new TH2F("hConePtSumClustervsTrack",
2712                                              Form("Track vs Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2713                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2714       fhConeSumPtClustervsTrack->SetXTitle("#Sigma #it{p}_{T}^{cluster} (GeV/#it{c})");
2715       fhConeSumPtClustervsTrack->SetYTitle("#Sigma #it{p}_{T}^{track} (GeV/#it{c})");
2716       outputContainer->Add(fhConeSumPtClustervsTrack) ;
2717
2718       fhConeSumPtClusterTrackFrac   = new TH2F("hConePtSumClusterTrackFraction",
2719                                              Form("#Sigma #it{p}_{T}^{cluster}/#Sigma #it{p}_{T}^{track} in isolation cone for #it{R} =  %2.2f",r),
2720                                              nptbins,ptmin,ptmax,200,0,5);
2721       fhConeSumPtClusterTrackFrac->SetYTitle("#Sigma #it{p}^{cluster}_{T} /#Sigma #it{p}_{T}^{track}");
2722       fhConeSumPtClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
2723       outputContainer->Add(fhConeSumPtClusterTrackFrac) ;
2724
2725       
2726       fhConePtLeadClustervsTrack   = new TH2F("hConePtLeadClustervsTrack",
2727                                              Form("Track vs Cluster lead #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2728                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2729       fhConePtLeadClustervsTrack->SetXTitle("#it{p}^{leading cluster}_{T} (GeV/#it{c})");
2730       fhConePtLeadClustervsTrack->SetYTitle("#it{p}^{leading track}_{T} (GeV/#it{c})");
2731       outputContainer->Add(fhConePtLeadClustervsTrack) ;
2732       
2733       fhConePtLeadClusterTrackFrac   = new TH2F("hConePtLeadClusterTrackFraction",
2734                                                Form(" #it{p}^{leading cluster}_{T}/#it{p}^{leading track}_{T} in isolation cone for #it{R} =  %2.2f",r),
2735                                                nptbins,ptmin,ptmax,200,0,5);
2736       fhConePtLeadClusterTrackFrac->SetYTitle("#it{p}^{leading cluster}_{T}/ #it{p}^{leading track}_{T}");
2737       fhConePtLeadClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
2738       outputContainer->Add(fhConePtLeadClusterTrackFrac) ;
2739
2740       
2741       if(fFillCellHistograms)
2742       {
2743         fhConeSumPtCellvsTrack   = new TH2F("hConePtSumCellvsTrack",
2744                                             Form("Track vs cell #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2745                                             nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2746         fhConeSumPtCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2747         fhConeSumPtCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2748         outputContainer->Add(fhConeSumPtCellvsTrack) ;
2749         
2750         fhConeSumPtCellTrack = new TH2F("hConePtSumCellTrack",
2751                                         Form("Track and Cell #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2752                                         nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2753         fhConeSumPtCellTrack->SetYTitle("#Sigma #it{p}_{T}");
2754         fhConeSumPtCellTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2755         outputContainer->Add(fhConeSumPtCellTrack) ;
2756         
2757         fhConeSumPtCellTrackTrigEtaPhi  = new TH2F("hConePtSumCellTrackTrigEtaPhi",
2758                                                    Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2759                                                    netabins,etamin,etamax,nphibins,phimin,phimax);
2760         fhConeSumPtCellTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2761         fhConeSumPtCellTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2762         fhConeSumPtCellTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2763         outputContainer->Add(fhConeSumPtCellTrackTrigEtaPhi) ;
2764         
2765       }
2766       
2767       if(fFillUEBandSubtractHistograms)
2768       {
2769         fhConeSumPtEtaUESub  = new TH2F("hConeSumPtEtaUESub",
2770                                         Form("#Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2771                                         nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2772         fhConeSumPtEtaUESub->SetYTitle("#Sigma #it{p}_{T}");
2773         fhConeSumPtEtaUESub->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2774         outputContainer->Add(fhConeSumPtEtaUESub) ;
2775         
2776         fhConeSumPtPhiUESub  = new TH2F("hConeSumPtPhiUESub",
2777                                         Form("#Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2778                                         nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2779         fhConeSumPtPhiUESub->SetYTitle("#Sigma #it{p}_{T}");
2780         fhConeSumPtPhiUESub->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2781         outputContainer->Add(fhConeSumPtPhiUESub) ;
2782         
2783         fhConeSumPtEtaUESubTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrigEtaPhi",
2784                                                   Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2785                                                   netabins,etamin,etamax,nphibins,phimin,phimax);
2786         fhConeSumPtEtaUESubTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2787         fhConeSumPtEtaUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2788         fhConeSumPtEtaUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2789         outputContainer->Add(fhConeSumPtEtaUESubTrigEtaPhi) ;
2790         
2791         fhConeSumPtPhiUESubTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrigEtaPhi",
2792                                                   Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2793                                                   netabins,etamin,etamax,nphibins,phimin,phimax);
2794         fhConeSumPtPhiUESubTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2795         fhConeSumPtPhiUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2796         fhConeSumPtPhiUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2797         outputContainer->Add(fhConeSumPtPhiUESubTrigEtaPhi) ;
2798         
2799         fhConeSumPtEtaUESubClustervsTrack   = new TH2F("hConePtSumEtaUESubClustervsTrack",
2800                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} =  %2.2f",r),
2801                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2802         fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2803         fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2804         outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2805         
2806         fhConeSumPtPhiUESubClustervsTrack   = new TH2F("hConePhiUESubPtSumClustervsTrack",
2807                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} =  %2.2f",r),
2808                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2809         fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2810         fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2811         outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2812         
2813         fhEtaBandClustervsTrack   = new TH2F("hEtaBandClustervsTrack",
2814                                              Form("Track vs Cluster #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2815                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2816         fhEtaBandClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2817         fhEtaBandClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2818         outputContainer->Add(fhEtaBandClustervsTrack) ;
2819         
2820         fhPhiBandClustervsTrack   = new TH2F("hPhiBandClustervsTrack",
2821                                              Form("Track vs Cluster #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2822                                              nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2823         fhPhiBandClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2824         fhPhiBandClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2825         outputContainer->Add(fhPhiBandClustervsTrack) ;
2826         
2827         fhEtaBandNormClustervsTrack   = new TH2F("hEtaBandNormClustervsTrack",
2828                                                  Form("Track vs Cluster #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2829                                                  nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2830         fhEtaBandNormClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2831         fhEtaBandNormClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2832         outputContainer->Add(fhEtaBandNormClustervsTrack) ;
2833         
2834         fhPhiBandNormClustervsTrack   = new TH2F("hPhiBandNormClustervsTrack",
2835                                                  Form("Track vs Cluster #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2836                                                  nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2837         fhPhiBandNormClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2838         fhPhiBandNormClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2839         outputContainer->Add(fhPhiBandNormClustervsTrack) ;
2840         
2841         fhConeSumPtEtaUESubClustervsTrack   = new TH2F("hConePtSumEtaUESubClustervsTrack",
2842                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} =  %2.2f",r),
2843                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2844         fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2845         fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2846         outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2847         
2848         fhConeSumPtPhiUESubClustervsTrack   = new TH2F("hConePhiUESubPtSumClustervsTrack",
2849                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} =  %2.2f",r),
2850                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2851         fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2852         fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2853         outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2854         
2855         if(fFillCellHistograms)
2856         {
2857           
2858           fhConeSumPtEtaUESubCellvsTrack   = new TH2F("hConePtSumEtaUESubCellvsTrack",
2859                                                       Form("Track vs Cell #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} =  %2.2f",r),
2860                                                       2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2861           fhConeSumPtEtaUESubCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2862           fhConeSumPtEtaUESubCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2863           outputContainer->Add(fhConeSumPtEtaUESubCellvsTrack) ;
2864           
2865           fhConeSumPtPhiUESubCellvsTrack   = new TH2F("hConePhiUESubPtSumCellvsTrack",
2866                                                       Form("Track vs Cell #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} =  %2.2f",r),
2867                                                       2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2868           fhConeSumPtPhiUESubCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2869           fhConeSumPtPhiUESubCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2870           outputContainer->Add(fhConeSumPtPhiUESubCellvsTrack) ;
2871           
2872           fhEtaBandCellvsTrack   = new TH2F("hEtaBandCellvsTrack",
2873                                             Form("Track vs Cell #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2874                                             nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2875           fhEtaBandCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2876           fhEtaBandCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2877           outputContainer->Add(fhEtaBandCellvsTrack) ;
2878           
2879           fhPhiBandCellvsTrack   = new TH2F("hPhiBandCellvsTrack",
2880                                             Form("Track vs Cell #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2881                                             nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2882           fhPhiBandCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2883           fhPhiBandCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2884           outputContainer->Add(fhPhiBandCellvsTrack) ;
2885           
2886           fhEtaBandNormCellvsTrack   = new TH2F("hEtaBandNormCellvsTrack",
2887                                                 Form("Track vs Cell #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2888                                                 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2889           fhEtaBandNormCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2890           fhEtaBandNormCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2891           outputContainer->Add(fhEtaBandNormCellvsTrack) ;
2892           
2893           fhPhiBandNormCellvsTrack   = new TH2F("hPhiBandNormCellvsTrack",
2894                                                 Form("Track vs Cell #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2895                                                 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2896           fhPhiBandNormCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2897           fhPhiBandNormCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2898           outputContainer->Add(fhPhiBandNormCellvsTrack) ;
2899           
2900           fhConeSumPtEtaUESubTrackCell  = new TH2F("hConeSumPtEtaUESubTrackCell",
2901                                                    Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2902                                                    nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2903           fhConeSumPtEtaUESubTrackCell->SetYTitle("#Sigma #it{p}_{T}");
2904           fhConeSumPtEtaUESubTrackCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2905           outputContainer->Add(fhConeSumPtEtaUESubTrackCell) ;
2906           
2907           fhConeSumPtPhiUESubTrackCell  = new TH2F("hConeSumPtPhiUESubTrackCell",
2908                                                    Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2909                                                    nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2910           fhConeSumPtPhiUESubTrackCell->SetYTitle("#Sigma #it{p}_{T}");
2911           fhConeSumPtPhiUESubTrackCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2912           outputContainer->Add(fhConeSumPtPhiUESubTrackCell) ;
2913           
2914           fhConeSumPtEtaUESubTrackCellTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrackCellTrigEtaPhi",
2915                                                              Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2916                                                              netabins,etamin,etamax,nphibins,phimin,phimax);
2917           fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2918           fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2919           fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2920           outputContainer->Add(fhConeSumPtEtaUESubTrackCellTrigEtaPhi) ;
2921           
2922           fhConeSumPtPhiUESubTrackCellTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrackCellTrigEtaPhi",
2923                                                              Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2924                                                              netabins,etamin,etamax,nphibins,phimin,phimax);
2925           fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2926           fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2927           fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2928           outputContainer->Add(fhConeSumPtPhiUESubTrackCellTrigEtaPhi) ;
2929         }
2930       }
2931     }
2932     
2933     for(Int_t iso = 0; iso < 2; iso++)
2934     {
2935       if(fFillTMHisto)
2936       {
2937         fhTrackMatchedDEta[iso]  = new TH2F
2938         (Form("hTrackMatchedDEta%s",isoName[iso].Data()),
2939          Form("%s - d#eta of cluster-track vs cluster energy, %s",isoTitle[iso].Data(),parTitle.Data()),
2940          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2941         fhTrackMatchedDEta[iso]->SetYTitle("d#eta");
2942         fhTrackMatchedDEta[iso]->SetXTitle("E_{cluster} (GeV)");
2943         
2944         fhTrackMatchedDPhi[iso]  = new TH2F
2945         (Form("hTrackMatchedDPhi%s",isoName[iso].Data()),
2946          Form("%s - d#phi of cluster-track vs cluster energy, %s",isoTitle[iso].Data(),parTitle.Data()),
2947          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2948         fhTrackMatchedDPhi[iso]->SetYTitle("d#phi (rad)");
2949         fhTrackMatchedDPhi[iso]->SetXTitle("E_{cluster} (GeV)");
2950         
2951         fhTrackMatchedDEtaDPhi[iso]  = new TH2F
2952         (Form("hTrackMatchedDEtaDPhi%s",isoName[iso].Data()),
2953          Form("%s - d#eta vs d#phi of cluster-track, %s",isoTitle[iso].Data(),parTitle.Data()),
2954          nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2955         fhTrackMatchedDEtaDPhi[iso]->SetYTitle("d#phi (rad)");
2956         fhTrackMatchedDEtaDPhi[iso]->SetXTitle("d#eta");
2957         
2958         outputContainer->Add(fhTrackMatchedDEta[iso]) ;
2959         outputContainer->Add(fhTrackMatchedDPhi[iso]) ;
2960         outputContainer->Add(fhTrackMatchedDEtaDPhi[iso]) ;
2961         
2962         fhdEdx[iso]  = new TH2F
2963         (Form("hdEdx%s",isoName[iso].Data()),
2964          Form("%s - Matched track <d#it{E}/d#it{x}> vs cluster #it{E}, %s",isoTitle[iso].Data(),parTitle.Data()),
2965          nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2966         fhdEdx[iso]->SetXTitle("#it{E} (GeV)");
2967         fhdEdx[iso]->SetYTitle("<d#it{E}/d#it{x}>");
2968         outputContainer->Add(fhdEdx[iso]);
2969         
2970         fhEOverP[iso]  = new TH2F
2971         (Form("hEOverP%s",isoName[iso].Data()),
2972          Form("%s - Matched track #it{E}/#it{p} vs cluster, %s",isoTitle[iso].Data(),parTitle.Data()),
2973          nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2974         fhEOverP[iso]->SetXTitle("#it{E} (GeV)");
2975         fhEOverP[iso]->SetYTitle("#it{E}/#it{p}");
2976         outputContainer->Add(fhEOverP[iso]);
2977         
2978         if(IsDataMC())
2979         {
2980           fhTrackMatchedMCParticle[iso]  = new TH2F
2981           (Form("hTrackMatchedMCParticle%s",isoName[iso].Data()),
2982            Form("%s - Origin of particle vs cluster #it{E}, %s",isoTitle[iso].Data(),parTitle.Data()),
2983            nptbins,ptmin,ptmax,8,0,8);
2984           fhTrackMatchedMCParticle[iso]->SetXTitle("#it{E} (GeV)");
2985           //fhTrackMatchedMCParticle[iso]->SetYTitle("Particle type");
2986           
2987           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(1 ,"Photon");
2988           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(2 ,"Electron");
2989           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
2990           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(4 ,"Rest");
2991           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
2992           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
2993           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
2994           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
2995           
2996           outputContainer->Add(fhTrackMatchedMCParticle[iso]);
2997         }
2998       }
2999       
3000       if(fFillSSHisto)
3001       {
3002         fhELambda0[iso]  = new TH2F
3003         (Form("hELambda0%s",isoName[iso].Data()),
3004          Form("%s cluster : #it{E} vs #lambda_{0}, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3005         fhELambda0[iso]->SetYTitle("#lambda_{0}^{2}");
3006         fhELambda0[iso]->SetXTitle("#it{E} (GeV)");
3007         outputContainer->Add(fhELambda0[iso]) ;
3008         
3009 //        fhELambda1[iso]  = new TH2F
3010 //        (Form("hELambda1%s",isoName[iso].Data()),
3011 //         Form("%s cluster: #it{E} vs #lambda_{1}, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3012 //        fhELambda1[iso]->SetYTitle("#lambda_{1}^{2}");
3013 //        fhELambda1[iso]->SetXTitle("#it{E} (GeV)");
3014 //        outputContainer->Add(fhELambda1[iso]) ;
3015         
3016         fhPtLambda0[iso]  = new TH2F
3017         (Form("hPtLambda0%s",isoName[iso].Data()),
3018          Form("%s cluster : #it{p}_{T} vs #lambda_{0}, %s",isoTitle[iso].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3019         fhPtLambda0[iso]->SetYTitle("#lambda_{0}^{2}");
3020         fhPtLambda0[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3021         outputContainer->Add(fhPtLambda0[iso]) ;
3022          
3023         if(IsDataMC())
3024         {
3025           for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
3026           {
3027             fhPtLambda0MC[imc][iso]  = new TH2F(Form("hPtLambda0%s_MC%s",isoName[iso].Data(),mcPartName[imc].Data()),
3028                                                 Form("%s cluster : #it{p}_{T} vs #lambda_{0}: %s %s",isoTitle[iso].Data(),mcPartType[imc].Data(),parTitle.Data()),
3029                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3030             fhPtLambda0MC[imc][iso]->SetYTitle("#lambda_{0}^{2}");
3031             fhPtLambda0MC[imc][iso]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3032             outputContainer->Add( fhPtLambda0MC[imc][iso]) ;
3033           }
3034         }
3035         
3036         if(fIsoDetector==kEMCAL &&  GetFirstSMCoveredByTRD() >= 0)
3037         {
3038           fhPtLambda0TRD[iso]  = new TH2F
3039           (Form("hPtLambda0TRD%s",isoName[iso].Data()),
3040            Form("%s cluster: #it{p}_{T} vs #lambda_{0}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3041           fhPtLambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
3042           fhPtLambda0TRD[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3043           outputContainer->Add(fhPtLambda0TRD[iso]) ;
3044           
3045           fhELambda0TRD[iso]  = new TH2F
3046           (Form("hELambda0TRD%s",isoName[iso].Data()),
3047            Form("%s cluster: #it{E} vs #lambda_{0}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3048           fhELambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
3049           fhELambda0TRD[iso]->SetXTitle("#it{E} (GeV)");
3050           outputContainer->Add(fhELambda0TRD[iso]) ;
3051           
3052 //          fhELambda1TRD[iso]  = new TH2F
3053 //          (Form("hELambda1TRD%s",isoName[iso].Data()),
3054 //           Form("%s cluster: #it{E} vs #lambda_{1}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3055 //          fhELambda1TRD[iso]->SetYTitle("#lambda_{1}^{2}");
3056 //          fhELambda1TRD[iso]->SetXTitle("#it{E} (GeV)");
3057 //          outputContainer->Add(fhELambda1TRD[iso]) ;
3058         }
3059         
3060         if(fFillNLMHistograms)
3061         {
3062           fhNLocMax[iso] = new TH2F
3063           (Form("hNLocMax%s",isoName[iso].Data()),
3064            Form("%s - Number of local maxima in cluster, %s",isoTitle[iso].Data(),parTitle.Data()),
3065            nptbins,ptmin,ptmax,10,0,10);
3066           fhNLocMax[iso]->SetYTitle("#it{NLM}");
3067           fhNLocMax[iso]->SetXTitle("#it{E} (GeV)");
3068           outputContainer->Add(fhNLocMax[iso]) ;
3069           
3070           fhELambda0LocMax1[iso]  = new TH2F
3071           (Form("hELambda0LocMax1%s",isoName[iso].Data()),
3072            Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{0}, #it{NLM}=1, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3073           fhELambda0LocMax1[iso]->SetYTitle("#lambda_{0}^{2}");
3074           fhELambda0LocMax1[iso]->SetXTitle("#it{E} (GeV)");
3075           outputContainer->Add(fhELambda0LocMax1[iso]) ;
3076           
3077           fhELambda1LocMax1[iso]  = new TH2F
3078           (Form("hELambda1LocMax1%s",isoName[iso].Data()),
3079            Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{1}, #it{NLM}=1, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3080           fhELambda1LocMax1[iso]->SetYTitle("#lambda_{1}^{2}");
3081           fhELambda1LocMax1[iso]->SetXTitle("#it{E} (GeV)");
3082           outputContainer->Add(fhELambda1LocMax1[iso]) ;
3083           
3084           fhELambda0LocMax2[iso]  = new TH2F
3085           (Form("hELambda0LocMax2%s",isoName[iso].Data()),
3086            Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{0}, #it{NLM}=2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3087           fhELambda0LocMax2[iso]->SetYTitle("#lambda_{0}^{2}");
3088           fhELambda0LocMax2[iso]->SetXTitle("#it{E} (GeV)");
3089           outputContainer->Add(fhELambda0LocMax2[iso]) ;
3090           
3091           fhELambda1LocMax2[iso]  = new TH2F
3092           (Form("hELambda1LocMax2%s",isoName[iso].Data()),
3093            Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{1}, #it{NLM}=2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3094           fhELambda1LocMax2[iso]->SetYTitle("#lambda_{1}^{2}");
3095           fhELambda1LocMax2[iso]->SetXTitle("#it{E} (GeV)");
3096           outputContainer->Add(fhELambda1LocMax2[iso]) ;
3097           
3098           fhELambda0LocMaxN[iso]  = new TH2F
3099           ( Form("hELambda0LocMaxN%s",isoName[iso].Data()),
3100            Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{0}, #it{NLM}>2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3101           fhELambda0LocMaxN[iso]->SetYTitle("#lambda_{0}^{2}");
3102           fhELambda0LocMaxN[iso]->SetXTitle("#it{E} (GeV)");
3103           outputContainer->Add(fhELambda0LocMaxN[iso]) ;
3104           
3105           fhELambda1LocMaxN[iso]  = new TH2F
3106           (Form("hELambda1LocMaxN%s",isoName[iso].Data()),
3107            Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{1}, #it{NLM}>2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3108           fhELambda1LocMaxN[iso]->SetYTitle("#lambda_{1}^{2}");
3109           fhELambda1LocMaxN[iso]->SetXTitle("#it{E} (GeV)");
3110           outputContainer->Add(fhELambda1LocMaxN[iso]) ;
3111         } // NLM
3112       } // SS histo
3113     } // control histograms for isolated and non isolated objects
3114     
3115     
3116     if(IsPileUpAnalysisOn())
3117     {
3118       fhPtTrackInConeOtherBC  = new TH2F("hPtTrackInConeOtherBC",
3119                                          Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC!=0",r),
3120                                          nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
3121       fhPtTrackInConeOtherBC->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
3122       fhPtTrackInConeOtherBC->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3123       outputContainer->Add(fhPtTrackInConeOtherBC) ;
3124       
3125       fhPtTrackInConeOtherBCPileUpSPD  = new TH2F("hPtTrackInConeOtherBCPileUpSPD",
3126                                                   Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC!=0, pile-up from SPD",r),
3127                                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
3128       fhPtTrackInConeOtherBCPileUpSPD->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
3129       fhPtTrackInConeOtherBCPileUpSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3130       outputContainer->Add(fhPtTrackInConeOtherBCPileUpSPD) ;
3131       
3132       fhPtTrackInConeBC0  = new TH2F("hPtTrackInConeBC0",
3133                                      Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC==0",r),
3134                                      nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
3135       fhPtTrackInConeBC0->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
3136       fhPtTrackInConeBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3137       outputContainer->Add(fhPtTrackInConeBC0) ;
3138       
3139       fhPtTrackInConeVtxBC0  = new TH2F("hPtTrackInConeVtxBC0",
3140                                         Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC==0",r),
3141                                         nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
3142     &n