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