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