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