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