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
  36 // --- Analysis system ---
  37 #include "AliAnaParticleIsolation.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliIsolationCut.h"
  40 #include "AliNeutralMesonSelection.h"
  41 #include "AliAODPWG4ParticleCorrelation.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "AliVTrack.h"
  44 #include "AliVCluster.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliEMCALGeometry.h"
  48
  49 ClassImp(AliAnaParticleIsolation)
  50
  51 //______________________________________________________________________________
  52 AliAnaParticleIsolation::AliAnaParticleIsolation() :
  53 AliAnaCaloTrackCorrBaseClass(),   fCalorimeter(""),
  54 fReMakeIC(0),                     fMakeSeveralIC(0),
  55 fFillPileUpHistograms(0),
  56 fFillTMHisto(0),                  fFillSSHisto(0),
  57 // Several IC
  58 fNCones(0),                       fNPtThresFrac(0),
  59 fConeSizes(),                     fPtThresholds(),
  60 fPtFractions(),                   fSumPtThresholds(),
  61 // Histograms
  62 fhEIso(0),                        fhPtIso(0),
  63 fhPtCentralityIso(0),             fhPtEventPlaneIso(0),
  64 fhPtNLocMaxIso(0),
  65 fhPhiIso(0),                      fhEtaIso(0),                              fhEtaPhiIso(0),
  66 fhEtaPhiNoIso(0),
  67 fhENoIso(0),                      fhPtNoIso(0),                             fhPtNLocMaxNoIso(0),
  68 fhPtDecayIso(0),                  fhPtDecayNoIso(0),
  69 fhEtaPhiDecayIso(0),              fhEtaPhiDecayNoIso(0),
  70 fhPtInCone(0),
  71 fhPtClusterInCone(0),             fhPtCellInCone(0),                        fhPtTrackInCone(0),
  72 fhPtTrackInConeOtherBC(0),        fhPtTrackInConeOtherBCPileUpSPD(0),
  73 fhPtTrackInConeBC0(0),            fhPtTrackInConeVtxBC0(0),
  74 fhPtTrackInConeBC0PileUpSPD(0),
  75 fhPtInConePileUp(),               fhPtInConeCent(0),
  76 fhPerpConeSumPt(0),               fhPtInPerpCone(0),
  77 fhEtaPhiInConeCluster(0),fhEtaPhiCluster(0),fhEtaPhiInConeTrack(0),fhEtaPhiTrack(0),
  78 fhEtaBandCluster(0),              fhPhiBandCluster(0),
  79 fhEtaBandTrack(0),                fhPhiBandTrack(0),
  80 fhEtaBandCell(0),                 fhPhiBandCell(0),
  81 fhConeSumPt(0),                   fhConeSumPtCellTrack(0),
  82 fhConeSumPtCell(0),               fhConeSumPtCluster(0),                    fhConeSumPtTrack(0),
  83 fhConeSumPtEtaBandUECluster(0),             fhConeSumPtPhiBandUECluster(0),
  84 fhConeSumPtEtaBandUETrack(0),               fhConeSumPtPhiBandUETrack(0),
  85 fhConeSumPtEtaBandUECell(0),                fhConeSumPtPhiBandUECell(0),
  86 fhConeSumPtTrigEtaPhi(0),
  87 fhConeSumPtCellTrackTrigEtaPhi(0),
  88 fhConeSumPtEtaBandUEClusterTrigEtaPhi(0),   fhConeSumPtPhiBandUEClusterTrigEtaPhi(0),
  89 fhConeSumPtEtaBandUETrackTrigEtaPhi(0),     fhConeSumPtPhiBandUETrackTrigEtaPhi(0),
  90 fhConeSumPtEtaBandUECellTrigEtaPhi(0),      fhConeSumPtPhiBandUECellTrigEtaPhi(0),
  91 fhConeSumPtEtaUESub(0),                     fhConeSumPtPhiUESub(0),
  92 fhConeSumPtEtaUESubTrigEtaPhi(0),           fhConeSumPtPhiUESubTrigEtaPhi(0),
  93 fhConeSumPtEtaUESubTrackCell(0),            fhConeSumPtPhiUESubTrackCell(0),
  94 fhConeSumPtEtaUESubTrackCellTrigEtaPhi(0),  fhConeSumPtPhiUESubTrackCellTrigEtaPhi(0),
  95 fhConeSumPtEtaUESubCluster(0),              fhConeSumPtPhiUESubCluster(0),
  96 fhConeSumPtEtaUESubClusterTrigEtaPhi(0),    fhConeSumPtPhiUESubClusterTrigEtaPhi(0),
  97 fhConeSumPtEtaUESubCell(0),                 fhConeSumPtPhiUESubCell(0),
  98 fhConeSumPtEtaUESubCellTrigEtaPhi(0),       fhConeSumPtPhiUESubCellTrigEtaPhi(0),
  99 fhConeSumPtEtaUESubTrack(0),                fhConeSumPtPhiUESubTrack(0),
 100 fhConeSumPtEtaUESubTrackTrigEtaPhi(0),      fhConeSumPtPhiUESubTrackTrigEtaPhi(0),
 101 fhFractionTrackOutConeEta(0),               fhFractionTrackOutConeEtaTrigEtaPhi(0),
 102 fhFractionClusterOutConeEta(0),             fhFractionClusterOutConeEtaTrigEtaPhi(0),
 103 fhFractionClusterOutConePhi(0),             fhFractionClusterOutConePhiTrigEtaPhi(0),
 104 fhFractionCellOutConeEta(0),                fhFractionCellOutConeEtaTrigEtaPhi(0),
 105 fhFractionCellOutConePhi(0),                fhFractionCellOutConePhiTrigEtaPhi(0),
 106 fhConeSumPtClustervsTrack(0),
 107 fhConeSumPtEtaUESubClustervsTrack(0),       fhConeSumPtPhiUESubClustervsTrack(0),
 108 fhConeSumPtCellvsTrack(0),
 109 fhConeSumPtEtaUESubCellvsTrack(0),          fhConeSumPtPhiUESubCellvsTrack(0),
 110 fhEtaBandClustervsTrack(0),                 fhPhiBandClustervsTrack(0),
 111 fhEtaBandNormClustervsTrack(0),             fhPhiBandNormClustervsTrack(0),
 112 fhEtaBandCellvsTrack(0),                    fhPhiBandCellvsTrack(0),
 113 fhEtaBandNormCellvsTrack(0),                fhPhiBandNormCellvsTrack(0),
 114 fhConeSumPtSubvsConeSumPtTotPhiTrack(0),    fhConeSumPtSubNormvsConeSumPtTotPhiTrack(0),
 115 fhConeSumPtSubvsConeSumPtTotEtaTrack(0),    fhConeSumPtSubNormvsConeSumPtTotEtaTrack(0),
 116 fhConeSumPtSubvsConeSumPtTotPhiCluster(0),  fhConeSumPtSubNormvsConeSumPtTotPhiCluster(0),
 117 fhConeSumPtSubvsConeSumPtTotEtaCluster(0),  fhConeSumPtSubNormvsConeSumPtTotEtaCluster(0),
 118 fhConeSumPtSubvsConeSumPtTotPhiCell(0),     fhConeSumPtSubNormvsConeSumPtTotPhiCell(0),
 119 fhConeSumPtSubvsConeSumPtTotEtaCell(0),     fhConeSumPtSubNormvsConeSumPtTotEtaCell(0),
 120 fhConeSumPtVSUETracksEtaBand(0),            fhConeSumPtVSUETracksPhiBand(0),
 121 fhConeSumPtVSUEClusterEtaBand(0),           fhConeSumPtVSUEClusterPhiBand(0),
 122 // MC histograms
 123 fhPtIsoPrompt(0),                 fhPhiIsoPrompt(0),               fhEtaIsoPrompt(0),
 124 fhPtThresIsolatedPrompt(),        fhPtFracIsolatedPrompt(),        fhPtSumIsolatedPrompt(),
 125 fhPtIsoFragmentation(0),          fhPhiIsoFragmentation(0),        fhEtaIsoFragmentation(0),
 126 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
 127 fhPtIsoPi0(0),                    fhPhiIsoPi0(0),                  fhEtaIsoPi0(0),
 128 fhPtThresIsolatedPi0(),           fhPtFracIsolatedPi0(),           fhPtSumIsolatedPi0(),
 129 fhPtIsoPi0Decay(0),               fhPhiIsoPi0Decay(0),             fhEtaIsoPi0Decay(0),
 130 fhPtThresIsolatedPi0Decay(),      fhPtFracIsolatedPi0Decay(),      fhPtSumIsolatedPi0Decay(),
 131 fhPtIsoEtaDecay(0),               fhPhiIsoEtaDecay(0),             fhEtaIsoEtaDecay(0),
 132 fhPtThresIsolatedEtaDecay(),      fhPtFracIsolatedEtaDecay(),      fhPtSumIsolatedEtaDecay(),
 133 fhPtIsoOtherDecay(0),             fhPhiIsoOtherDecay(0),           fhEtaIsoOtherDecay(0),
 134 fhPtThresIsolatedOtherDecay(),    fhPtFracIsolatedOtherDecay(),    fhPtSumIsolatedOtherDecay(),
 135 //fhPtIsoConversion(0),             fhPhiIsoConversion(0),           fhEtaIsoConversion(0),
 136 //fhPtThresIsolatedConversion(),    fhPtFracIsolatedConversion(),    fhPtSumIsolatedConversion(),
 137 fhPtIsoHadron(0),                 fhPhiIsoHadron(0),               fhEtaIsoHadron(0),
 138 fhPtThresIsolatedHadron(),        fhPtFracIsolatedHadron(),        fhPtSumIsolatedHadron(),
 139 fhPtNoIsoPi0(0),                  fhPtNoIsoPi0Decay(0),
 140 fhPtNoIsoEtaDecay(0),             fhPtNoIsoOtherDecay(0),
 141 fhPtNoIsoPrompt(0),               fhPtIsoMCPhoton(0),              fhPtNoIsoMCPhoton(0),
 142 //fhPtNoIsoConversion(0),
 143 fhPtNoIsoFragmentation(0),        fhPtNoIsoHadron(0),
 144 // Hist several IC
 145 fhSumPtLeadingPt(),               fhPtLeadingPt(),
 146 fhPerpSumPtLeadingPt(),           fhPerpPtLeadingPt(),
 147 fhPtThresIsolated(),              fhPtFracIsolated(),              fhPtSumIsolated(),
 148 fhEtaPhiPtThresIso(),             fhEtaPhiPtThresDecayIso(),       fhPtPtThresDecayIso(),
 149 fhEtaPhiPtFracIso(),              fhEtaPhiPtFracDecayIso(),        fhPtPtFracDecayIso(),
 150 fhPtPtSumDecayIso(),              fhEtaPhiSumDensityIso(),         fhEtaPhiSumDensityDecayIso(),
 151 fhPtSumDensityIso(),              fhPtSumDensityDecayIso(),
 152 fhPtFracPtSumIso(),               fhPtFracPtSumDecayIso(),
 153 fhEtaPhiFracPtSumIso(),           fhEtaPhiFracPtSumDecayIso(),
 154 // Cluster control histograms
 155 fhTrackMatchedDEta(),             fhTrackMatchedDPhi(),           fhTrackMatchedDEtaDPhi(),
 156 fhdEdx(),                         fhEOverP(),                     fhTrackMatchedMCParticle(),
 157 fhELambda0() ,                    fhELambda1(),                   fhELambda0SSBkg(),
 158 fhELambda0TRD(),                  fhELambda1TRD(),
 159 fhELambda0MCPhoton(),             fhELambda0MCPi0(),              fhELambda0MCPi0Decay(),
 160 fhELambda0MCEtaDecay(),           fhELambda0MCOtherDecay(),       fhELambda0MCHadron(),
 161 // Number of local maxima in cluster
 162 fhNLocMax(),
 163 fhELambda0LocMax1(),              fhELambda1LocMax1(),
 164 fhELambda0LocMax2(),              fhELambda1LocMax2(),
 165 fhELambda0LocMaxN(),              fhELambda1LocMaxN(),
 166 // PileUp
 167 fhEIsoPileUp(),                   fhPtIsoPileUp(),
 168 fhENoIsoPileUp(),                 fhPtNoIsoPileUp(),
 169 fhTimeENoCut(0),                  fhTimeESPD(0),                  fhTimeESPDMulti(0),
 170 fhTimeNPileUpVertSPD(0),          fhTimeNPileUpVertTrack(0),
 171 fhTimeNPileUpVertContributors(0),
 172 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
 173 // Histograms settings
 174 fHistoNPtSumBins(0),              fHistoPtSumMax(0.),              fHistoPtSumMin(0.),
 175 fHistoNPtInConeBins(0),           fHistoPtInConeMax(0.),           fHistoPtInConeMin(0.)
 176 {
 177   //default ctor
 178
 179   //Initialize parameters
 180   InitParameters();
 181
 182   for(Int_t i = 0; i < 5 ; i++)
 183   {
 184     fConeSizes[i]      = 0 ;
 185
 186     fhPtSumIsolatedPrompt       [i] = 0 ;
 187     fhPtSumIsolatedFragmentation[i] = 0 ;
 188     fhPtSumIsolatedPi0Decay     [i] = 0 ;
 189     fhPtSumIsolatedPi0          [i] = 0 ;
 190     fhPtSumIsolatedEtaDecay     [i] = 0 ;
 191     fhPtSumIsolatedOtherDecay   [i] = 0 ;
 192 //  fhPtSumIsolatedConversion   [i] = 0 ;
 193     fhPtSumIsolatedHadron       [i] = 0 ;
 194
 195     for(Int_t j = 0; j < 5 ; j++)
 196     {
 197       fhPtThresIsolated             [i][j] = 0 ;
 198       fhPtFracIsolated              [i][j] = 0 ;
 199       fhPtSumIsolated               [i][j] = 0 ;
 200
 201       fhEtaPhiPtThresIso            [i][j] = 0 ;
 202       fhEtaPhiPtThresDecayIso       [i][j] = 0 ;
 203       fhPtPtThresDecayIso           [i][j] = 0 ;
 204
 205       fhEtaPhiPtFracIso             [i][j] = 0 ;
 206       fhEtaPhiPtFracDecayIso        [i][j] = 0 ;
 207       fhPtPtFracDecayIso            [i][j] = 0 ;
 208       fhPtPtSumDecayIso             [i][j] = 0 ;
 209       fhPtSumDensityIso             [i][j] = 0 ;
 210       fhPtSumDensityDecayIso        [i][j] = 0 ;
 211       fhEtaPhiSumDensityIso         [i][j] = 0 ;
 212       fhEtaPhiSumDensityDecayIso    [i][j] = 0 ;
 213       fhPtFracPtSumIso              [i][j] = 0 ;
 214       fhPtFracPtSumDecayIso         [i][j] = 0 ;
 215       fhEtaPhiFracPtSumIso          [i][j] = 0 ;
 216       fhEtaPhiFracPtSumDecayIso     [i][j] = 0 ;
 217
 218       fhPtThresIsolatedPrompt       [i][j] = 0 ;
 219       fhPtThresIsolatedFragmentation[i][j] = 0 ;
 220       fhPtThresIsolatedPi0Decay     [i][j] = 0 ;
 221       fhPtThresIsolatedPi0          [i][j] = 0 ;
 222       fhPtThresIsolatedEtaDecay     [i][j] = 0 ;
 223       fhPtThresIsolatedOtherDecay   [i][j] = 0 ;
 224 //    fhPtThresIsolatedConversion   [i][j] = 0 ;
 225       fhPtThresIsolatedHadron      [ i][j] = 0 ;
 226
 227       fhPtFracIsolatedPrompt        [i][j] = 0 ;
 228       fhPtFracIsolatedFragmentation [i][j] = 0 ;
 229       fhPtFracIsolatedPi0           [i][j] = 0 ;
 230       fhPtFracIsolatedPi0Decay      [i][j] = 0 ;
 231       fhPtFracIsolatedEtaDecay      [i][j] = 0 ;
 232       fhPtFracIsolatedOtherDecay    [i][j] = 0 ;
 233 //    fhPtFracIsolatedConversion    [i][j] = 0 ;
 234       fhPtFracIsolatedHadron        [i][j] = 0 ;
 235
 236     }
 237   }
 238
 239   for(Int_t i = 0; i < 5 ; i++)
 240   {
 241     fPtFractions    [i] = 0 ;
 242     fPtThresholds   [i] = 0 ;
 243     fSumPtThresholds[i] = 0 ;
 244   }
 245
 246
 247   for(Int_t i = 0; i < 2 ; i++)
 248   {
 249     fhTrackMatchedDEta[i] = 0 ;             fhTrackMatchedDPhi[i] = 0 ;           fhTrackMatchedDEtaDPhi  [i] = 0 ;
 250     fhdEdx            [i] = 0 ;             fhEOverP          [i] = 0 ;           fhTrackMatchedMCParticle[i] = 0 ;
 251     fhELambda0        [i] = 0 ;             fhELambda1        [i] = 0 ;
 252     fhELambda0TRD     [i] = 0 ;             fhELambda1TRD     [i] = 0 ;
 253
 254     fhELambda0MCPhoton  [i] = 0 ;           fhELambda0MCPi0       [i] = 0 ;       fhELambda0MCPi0Decay[i] = 0 ;
 255     fhELambda0MCEtaDecay[i] = 0 ;           fhELambda0MCOtherDecay[i] = 0 ;       fhELambda0MCHadron  [i] = 0 ;
 256
 257
 258     // Number of local maxima in cluster
 259     fhNLocMax        [i] = 0 ;
 260     fhELambda0LocMax1[i] = 0 ;              fhELambda1LocMax1[i] = 0 ;
 261     fhELambda0LocMax2[i] = 0 ;              fhELambda1LocMax2[i] = 0 ;
 262     fhELambda0LocMaxN[i] = 0 ;              fhELambda1LocMaxN[i] = 0 ;
 263
 264   }
 265
 266   // Pile-Up
 267
 268   for(Int_t i = 0 ; i < 7 ; i++)
 269   {
 270     fhPtInConePileUp[i] = 0 ;
 271     fhEIsoPileUp    [i] = 0 ;
 272     fhPtIsoPileUp   [i] = 0 ;
 273     fhENoIsoPileUp  [i] = 0 ;
 274     fhPtNoIsoPileUp [i] = 0 ;
 275   }
 276
 277 }
 278
 279 //_______________________________________________________________________________________________
 280 void AliAnaParticleIsolation::CalculateCaloUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
 281                                                   Float_t & etaBandPtSum, Float_t & phiBandPtSum)
 282 {
 283   // Get the clusters pT or sum of pT in phi/eta bands or at 45 degrees from trigger
 284
 285   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
 286
 287   Float_t conesize   = GetIsolationCut()->GetConeSize();
 288   TLorentzVector mom ;
 289
 290   //Select the Calorimeter
 291   TObjArray * pl = 0x0;
 292   if      (fCalorimeter == "PHOS" )
 293     pl    = GetPHOSClusters();
 294   else if (fCalorimeter == "EMCAL")
 295     pl    = GetEMCALClusters();
 296
 297   if(!pl) return ;
 298
 299   //Get vertex for cluster momentum calculation
 300   Double_t vertex[] = {0,0,0} ; //vertex ;
 301   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 302     GetReader()->GetVertex(vertex);
 303
 304   Float_t ptTrig    = pCandidate->Pt() ;
 305   Float_t phiTrig   = pCandidate->Phi();
 306   Float_t etaTrig   = pCandidate->Eta();
 307
 308   for(Int_t icluster=0; icluster < pl->GetEntriesFast(); icluster++)
 309   {
 310     AliVCluster* cluster = (AliVCluster *) pl->At(icluster);
 311
 312     if(!cluster)
 313     {
 314       printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Cluster not available?");
 315       continue;
 316     }
 317
 318     //Do not count the candidate (photon or pi0) or the daughters of the candidate
 319     if(cluster->GetID() == pCandidate->GetCaloLabel(0) ||
 320        cluster->GetID() == pCandidate->GetCaloLabel(1)   ) continue ;
 321
 322     //Remove matched clusters to tracks if Neutral and Track info is used
 323     if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged &&
 324         IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ;
 325
 326     cluster->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
 327
 328     //exclude particles in cone
 329     Float_t rad = GetIsolationCut()->Radius(etaTrig, phiTrig, mom.Eta(), mom.Phi());
 330
 331     // histo of eta and phi for all clusters
 332     fhEtaPhiCluster->Fill(mom.Eta(), mom.Phi());
 333     if(rad < conesize) {
 334         // histos for all clusters in cone
 335      fhEtaPhiInConeCluster->Fill(mom.Eta(), mom.Phi());
 336     continue ;
 337     }
 338     //fill histogram for UE in phi band in EMCal acceptance
 339     if(mom.Eta() > (etaTrig-conesize) && mom.Eta()  < (etaTrig+conesize))
 340     {
 341       phiBandPtSum+=mom.Pt();
 342       fhPhiBandCluster->Fill(mom.Eta(),mom.Phi());
 343
 344 }
 345
 346     //fill histogram for UE in eta band in EMCal acceptance
 347     if(mom.Phi() > (phiTrig-conesize) && mom.Phi() < (phiTrig+conesize))
 348     {
 349       etaBandPtSum+=mom.Pt();
 350       fhEtaBandCluster->Fill(mom.Eta(),mom.Phi());
 351     }
 352   }
 353
 354   fhConeSumPtEtaBandUECluster          ->Fill(ptTrig  ,       etaBandPtSum);
 355   fhConeSumPtPhiBandUECluster          ->Fill(ptTrig  ,       phiBandPtSum);
 356   fhConeSumPtEtaBandUEClusterTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSum);
 357   fhConeSumPtPhiBandUEClusterTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSum);
 358
 359 }
 360
 361 //________________________________________________________________________________________________
 362 void AliAnaParticleIsolation::CalculateCaloCellUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
 363                                                       Float_t & etaBandPtSumCells, Float_t & phiBandPtSumCells)
 364 {
 365   // Get the cells amplitude or sum of amplitude in phi/eta bands or at 45 degrees from trigger
 366
 367   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
 368
 369   Float_t conesize = GetIsolationCut()->GetConeSize();
 370
 371   Float_t phiTrig = pCandidate->Phi();
 372   if(phiTrig<0) phiTrig += TMath::TwoPi();
 373   Float_t etaTrig = pCandidate->Eta();
 374
 375   if(pCandidate->GetDetector()=="EMCAL")
 376   {
 377     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
 378     Int_t absId = -999;
 379
 380     if (eGeom->GetAbsCellIdFromEtaPhi(etaTrig,phiTrig,absId))
 381     {
 382       if(!eGeom->CheckAbsCellId(absId)) return ;
 383
 384       // Get absolute (col,row) of trigger particle
 385       Int_t nSupMod = eGeom->GetSuperModuleNumber(absId);
 386       Int_t nModule = -1;
 387       Int_t imEta=-1, imPhi=-1;
 388       Int_t ieta =-1, iphi =-1;
 389
 390       if (eGeom->GetCellIndex(absId,nSupMod,nModule,imPhi,imEta))
 391       {
 392         eGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,imPhi,imEta,iphi,ieta);
 393
 394         Int_t colTrig = ieta;
 395         if (nSupMod % 2) colTrig = AliEMCALGeoParams::fgkEMCALCols + ieta ;
 396         Int_t rowTrig = iphi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
 397
 398         Int_t sqrSize = int(conesize/0.0143);
 399
 400         AliVCaloCells * cells = GetEMCALCells();
 401
 402         Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 24*(16/3) 5 full-size Sectors (2 SM) + 1 one-third Sector (2 SM)
 403         Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols;
 404         //  printf("nTotalRows %i, nTotalCols %i\n",nTotalRows,nTotalCols);
 405         // Loop on cells in eta band
 406
 407                   Int_t irowmin = rowTrig-sqrSize;
 408           if(irowmin<0) irowmin=0;
 409   Int_t irowmax = rowTrig+sqrSize;
 410           if(irowmax>AliEMCALGeoParams::fgkEMCALRows) irowmax=AliEMCALGeoParams::fgkEMCALRows;
 411
 412
 413         for(Int_t irow = irowmin; irow <irowmax; irow++)
 414         {
 415           for(Int_t icol = 0; icol < nTotalCols; icol++)
 416           {
 417             Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
 418         if(inSector==5) continue;
 419           Int_t inSupMod = -1;
 420             Int_t icolLoc  = -1;
 421             if(icol < AliEMCALGeoParams::fgkEMCALCols)
 422             {
 423               inSupMod = 2*inSector + 1;
 424               icolLoc  = icol;
 425             }
 426             else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
 427             {
 428               inSupMod = 2*inSector;
 429               icolLoc  = icol-AliEMCALGeoParams::fgkEMCALCols;
 430             }
 431
 432             Int_t irowLoc  = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;
 433
 434             // Exclude cells in cone
 435             if(TMath::Abs(icol-colTrig) < sqrSize || TMath::Abs(irow-rowTrig) < sqrSize){
 436                  continue ;
 437             }
 438             Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
 439             if(!eGeom->CheckAbsCellId(iabsId)) continue;
 440             etaBandPtSumCells += cells->GetCellAmplitude(iabsId);
 441                 fhEtaBandCell->Fill(colTrig,rowTrig);
 442
 443     //          printf("ETA inSupMod %i,irowLoc %i,icolLoc %i, iabsId %i, etaBandPtSumCells %f\n",nSupMod,irowLoc,icolLoc,iabsId,etaBandPtSumCells);
 444   }
 445         }
 446     Int_t icolmin = colTrig-sqrSize;
 447           if(icolmin<0) icolmin=0;
 448           Int_t icolmax = colTrig+sqrSize;
 449           if(icolmax>AliEMCALGeoParams::fgkEMCALCols) icolmax=AliEMCALGeoParams::fgkEMCALCols;
 450
 451         // Loop on cells in phi band
 452         for(Int_t icol = icolmin; icol < icolmax; icol++)
 453         {
 454           for(Int_t irow = 0; irow < nTotalRows; irow++)
 455           {
 456             Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
 457                 if(inSector==5) continue;
 458             Int_t inSupMod = -1;
 459             Int_t icolLoc  = -1;
 460         //    printf("icol %i, irow %i, inSector %i\n",icol,irow ,inSector);
 461             if(icol < AliEMCALGeoParams::fgkEMCALCols)
 462             {
 463         //      printf("icol < AliEMCALGeoParams::fgkEMCALCols %i\n",AliEMCALGeoParams::fgkEMCALCols );
 464               inSupMod = 2*inSector + 1;
 465               icolLoc  = icol;
 466             }
 467             else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
 468             {
 469      //      printf("icol > AliEMCALGeoParams::fgkEMCALCols -1 %i\n",AliEMCALGeoParams::fgkEMCALCols -1 );
 470              inSupMod = 2*inSector;
 471               icolLoc  = icol-AliEMCALGeoParams::fgkEMCALCols;
 472             }
 473
 474             Int_t irowLoc  = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;   // Stesso problema di sopra //
 475
 476             // Exclude cells in cone
 477             if(TMath::Abs(icol-colTrig) < sqrSize) {
 478                 //printf("TMath::Abs(icol-colTrig) %i < sqrSize %i\n",TMath::Abs(icol-colTrig) ,sqrSize);continue ;
 479                 }
 480             if(TMath::Abs(irow-rowTrig) < sqrSize) {
 481                 //printf("TMath::Abs(irow-rowTrig) %i < sqrSize %i\n",TMath::Abs(irow-rowTrig) ,sqrSize);continue ;
 482                 }
 483
 484             Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
 485             if(!eGeom->CheckAbsCellId(iabsId)) {printf("!eGeom->CheckAbsCellId(iabsId=%i) inSupMod %i irowLoc %i icolLoc %i \n",iabsId,inSupMod, irowLoc, icolLoc);continue;}
 486             phiBandPtSumCells += cells->GetCellAmplitude(iabsId);
 487                 fhPhiBandCell->Fill(colTrig,rowTrig);
 488          //printf("inSupMod %i,irowLoc %i,icolLoc %i, iabsId %i, phiBandPtSumCells %f\n",nSupMod,irowLoc,icolLoc,iabsId,phiBandPtSumCells);
 489             }
 490         }
 491       }
 492     }
 493   }
 494
 495   Float_t ptTrig = pCandidate->Pt();
 496
 497   fhConeSumPtEtaBandUECell          ->Fill(ptTrig ,        etaBandPtSumCells);
 498   fhConeSumPtPhiBandUECell          ->Fill(ptTrig ,        phiBandPtSumCells);
 499   fhConeSumPtEtaBandUECellTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSumCells);
 500   fhConeSumPtPhiBandUECellTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSumCells);
 501
 502 }
 503
 504 //________________________________________________________________________________________________
 505 void AliAnaParticleIsolation::CalculateTrackUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
 506                                                    Float_t & etaBandPtSum, Float_t & phiBandPtSum)
 507 {
 508   // Get the track pT or sum of pT in phi/eta bands or at 45 degrees from trigger
 509
 510   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyNeutral ) return ;
 511
 512   Float_t conesize   = GetIsolationCut()->GetConeSize();
 513
 514   Double_t sumptPerp= 0. ;
 515   Float_t ptTrig    = pCandidate->Pt() ;
 516   Float_t phiTrig   = pCandidate->Phi();
 517   Float_t etaTrig   = pCandidate->Eta();
 518
 519   TObjArray * trackList   = GetCTSTracks() ;
 520   for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
 521   {
 522     AliVTrack* track = (AliVTrack *) trackList->At(itrack);
 523
 524     if(!track)
 525     {
 526       printf("AliAnaParticleIsolation::CalculateTrackUEBand() - Track not available?");
 527       continue;
 528     }
 529
 530     //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
 531     if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
 532        track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)   ) continue ;
 533
 534    // histo of eta:phi for all tracks
 535     fhEtaPhiTrack->Fill(track->Eta(),track->Phi());
 536
 537     //exclude particles in cone
 538     Float_t rad = GetIsolationCut()->Radius(etaTrig, phiTrig, track->Eta(), track->Phi());
 539     if(rad < conesize) {
 540         // histo of eta:phi for all tracks in cone
 541         fhEtaPhiInConeTrack->Fill(track->Eta(),track->Phi());
 542       continue ;
 543     }
 544
 545     //fill histogram for UE in phi band
 546     if(track->Eta() > (etaTrig-conesize) && track->Eta()  < (etaTrig+conesize))
 547     {
 548       phiBandPtSum+=track->Pt();
 549       fhPhiBandTrack->Fill(track->Eta(),track->Phi());
 550     }
 551
 552     //fill histogram for UE in eta band in EMCal acceptance
 553     if(track->Phi() > (phiTrig-conesize) && track->Phi() < (phiTrig+conesize))
 554     {
 555       etaBandPtSum+=track->Pt();
 556       fhEtaBandTrack->Fill(track->Eta(),track->Phi());
 557     }
 558
 559     //fill the histograms at +-45 degrees in phi from trigger particle, perpedicular to trigger axis in phi
 560     Double_t dPhi = phiTrig - track->Phi() + TMath::PiOver2();
 561     Double_t dEta = etaTrig - track->Eta();
 562     Double_t arg  = dPhi*dPhi + dEta*dEta;
 563     if(TMath::Sqrt(arg) < conesize)
 564     {
 565       fhPtInPerpCone->Fill(ptTrig,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
 566       sumptPerp+=track->Pt();
 567     }
 568
 569     dPhi = phiTrig - track->Phi() - TMath::PiOver2();
 570     arg  = dPhi*dPhi + dEta*dEta;
 571     if(TMath::Sqrt(arg) < conesize)
 572     {
 573       fhPtInPerpCone->Fill(ptTrig,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
 574       sumptPerp+=track->Pt();
 575     }
 576   }
 577
 578   fhPerpConeSumPt                    ->Fill(ptTrig ,        sumptPerp   );
 579   fhConeSumPtEtaBandUETrack          ->Fill(ptTrig ,        etaBandPtSum);
 580   fhConeSumPtPhiBandUETrack          ->Fill(ptTrig ,        phiBandPtSum);
 581   fhConeSumPtEtaBandUETrackTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSum);
 582   fhConeSumPtPhiBandUETrackTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSum);
 583
 584 }
 585
 586
 587
 588 //___________________________________________________________________________________________________________________________________
 589 void AliAnaParticleIsolation::CalculateNormalizeUEBandPerUnitArea(AliAODPWG4ParticleCorrelation * pCandidate,
 590                                                                   Float_t coneptsumCluster,  Float_t coneptsumCell, Float_t coneptsumTrack, Float_t &etaBandptsumTrackNorm, Float_t &etaBandptsumClusterNorm)
 591 {
 592   //normalize phi/eta band per area unit
 593
 594   Float_t etaUEptsumTrack   = 0 ;
 595   Float_t phiUEptsumTrack   = 0 ;
 596   Float_t etaUEptsumCluster = 0 ;
 597   Float_t phiUEptsumCluster = 0 ;
 598   Float_t etaUEptsumCell    = 0 ;
 599   Float_t phiUEptsumCell    = 0 ;
 600
 601   Int_t   partTypeInCone    = GetIsolationCut()->GetParticleTypeInCone();
 602
 603   // Do the normalization
 604
 605   Float_t conesize  = GetIsolationCut()->GetConeSize();
 606   Float_t coneA     = conesize*conesize*TMath::Pi(); // A = pi R^2, isolation cone area
 607   Float_t ptTrig    = pCandidate->Pt() ;
 608   Float_t phiTrig   = pCandidate->Phi();
 609   Float_t etaTrig   = pCandidate->Eta();
 610
 611
 612   // ------ //
 613   // Tracks //
 614   // ------ //
 615   Float_t phiUEptsumTrackNorm  = 0 ;
 616   Float_t etaUEptsumTrackNorm  = 0 ;
 617   Float_t coneptsumTrackSubPhi = 0 ;
 618   Float_t coneptsumTrackSubEta = 0 ;
 619   Float_t coneptsumTrackSubPhiNorm = 0 ;
 620   Float_t coneptsumTrackSubEtaNorm = 0 ;
 621   etaBandptsumTrackNorm = 0 ;
 622
 623   if( partTypeInCone!=AliIsolationCut::kOnlyNeutral )
 624   {
 625     // Sum the pT in the phi or eta band for clusters or tracks
 626     CalculateTrackUEBand   (pCandidate,etaUEptsumTrack  ,phiUEptsumTrack  );// rajouter ici l'histo eta phi
 627
 628   //Fill histos
 629   fhConeSumPtVSUETracksEtaBand->Fill(coneptsumTrack,etaUEptsumTrack);
 630   fhConeSumPtVSUETracksPhiBand->Fill(coneptsumTrack,phiUEptsumTrack);
 631
 632
 633     Float_t correctConeSumTrack    = 1;
 634     Float_t correctConeSumTrackPhi = 1;
 635
 636     GetIsolationCut()->CalculateUEBandTrackNormalization(GetReader(),etaTrig, phiTrig,
 637                                                            phiUEptsumTrack,etaUEptsumTrack,
 638                                                            phiUEptsumTrackNorm,etaUEptsumTrackNorm,
 639                                                            correctConeSumTrack,correctConeSumTrackPhi);
 640
 641     coneptsumTrackSubPhi = coneptsumTrack - phiUEptsumTrackNorm;
 642     coneptsumTrackSubEta = coneptsumTrack - etaUEptsumTrackNorm;
 643
 644     etaBandptsumTrackNorm = etaUEptsumTrackNorm;
 645
 646     fhConeSumPtPhiUESubTrack           ->Fill(ptTrig ,          coneptsumTrackSubPhi);
 647     fhConeSumPtPhiUESubTrackTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumTrackSubPhi);
 648     fhConeSumPtEtaUESubTrack           ->Fill(ptTrig ,          coneptsumTrackSubEta);
 649     fhConeSumPtEtaUESubTrackTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumTrackSubEta);
 650
 651     fhFractionTrackOutConeEta          ->Fill(ptTrig ,         correctConeSumTrack-1);
 652     fhFractionTrackOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig,correctConeSumTrack-1);
 653
 654     if(coneptsumTrack > 0)
 655     {
 656       coneptsumTrackSubPhiNorm = coneptsumTrackSubPhi/coneptsumTrack;
 657         coneptsumTrackSubEtaNorm = coneptsumTrackSubEta/coneptsumTrack;
 658     }
 659
 660     fhConeSumPtSubvsConeSumPtTotPhiTrack    ->Fill(coneptsumTrack,coneptsumTrackSubPhi);
 661     fhConeSumPtSubNormvsConeSumPtTotPhiTrack->Fill(coneptsumTrack,coneptsumTrackSubPhiNorm);
 662     fhConeSumPtSubvsConeSumPtTotEtaTrack    ->Fill(coneptsumTrack,coneptsumTrackSubEta);
 663     fhConeSumPtSubNormvsConeSumPtTotEtaTrack->Fill(coneptsumTrack,coneptsumTrackSubEtaNorm);
 664
 665   }
 666
 667   // ------------------------ //
 668   // EMCal Clusters and cells //
 669   // ------------------------ //
 670   Float_t phiUEptsumClusterNorm  = 0 ;
 671   Float_t etaUEptsumClusterNorm  = 0 ;
 672   Float_t coneptsumClusterSubPhi = 0 ;
 673   Float_t coneptsumClusterSubEta = 0 ;
 674   Float_t coneptsumClusterSubPhiNorm = 0 ;
 675   Float_t coneptsumClusterSubEtaNorm = 0 ;
 676   Float_t phiUEptsumCellNorm     = 0 ;
 677   Float_t etaUEptsumCellNorm     = 0 ;
 678   Float_t coneptsumCellSubPhi    = 0 ;
 679   Float_t coneptsumCellSubEta    = 0 ;
 680   Float_t coneptsumCellSubPhiNorm = 0 ;
 681   Float_t coneptsumCellSubEtaNorm = 0 ;
 682   etaBandptsumClusterNorm = 0;
 683
 684   if( partTypeInCone!=AliIsolationCut::kOnlyCharged )
 685   {
 686
 687     // -------------- //
 688     // EMCal clusters //
 689     // -------------- //
 690
 691     // Sum the pT in the phi or eta band for clusters or tracks
 692     CalculateCaloUEBand    (pCandidate,etaUEptsumCluster,phiUEptsumCluster);// rajouter ici l'histo eta phi
 693
 694   //Fill histos
 695   fhConeSumPtVSUEClusterEtaBand->Fill(coneptsumCluster,etaUEptsumCluster);
 696   fhConeSumPtVSUEClusterPhiBand->Fill(coneptsumCluster,phiUEptsumCluster);
 697
 698
 699     Float_t correctConeSumClusterEta = 1;
 700     Float_t correctConeSumClusterPhi = 1;
 701
 702     GetIsolationCut()->CalculateUEBandClusterNormalization(GetReader(),etaTrig, phiTrig,
 703                                                            phiUEptsumCluster,etaUEptsumCluster,
 704                                                            phiUEptsumClusterNorm,etaUEptsumClusterNorm,
 705                                                            correctConeSumClusterEta,correctConeSumClusterPhi);
 706
 707     // In case that cone is out of eta and phi side, we are over correcting, not too often with the current cuts ...
 708     // Comment if not used
 709     //  Float_t coneBadCellsCoeff   =1;
 710     //  Float_t etaBandBadCellsCoeff=1;
 711     //  Float_t phiBandBadCellsCoeff=1;
 712     //  GetIsolationCut()->GetCoeffNormBadCell(pCandidate,   GetReader(),coneBadCellsCoeff,etaBandBadCellsCoeff,phiBandBadCellsCoeff) ;
 713
 714     //coneptsumCluster=coneptsumCluster*coneBadCellsCoeff*correctConeSumClusterEta*correctConeSumClusterPhi;
 715
 716     coneptsumClusterSubPhi = coneptsumCluster - phiUEptsumClusterNorm;
 717     coneptsumClusterSubEta = coneptsumCluster - etaUEptsumClusterNorm;
 718
 719     etaBandptsumClusterNorm = etaUEptsumClusterNorm;
 720
 721     fhConeSumPtPhiUESubCluster           ->Fill(ptTrig ,          coneptsumClusterSubPhi);
 722     fhConeSumPtPhiUESubClusterTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumClusterSubPhi);
 723     fhConeSumPtEtaUESubCluster           ->Fill(ptTrig ,          coneptsumClusterSubEta);
 724     fhConeSumPtEtaUESubClusterTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumClusterSubEta);
 725
 726     fhFractionClusterOutConeEta          ->Fill(ptTrig ,          correctConeSumClusterEta-1);
 727     fhFractionClusterOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumClusterEta-1);
 728     fhFractionClusterOutConePhi          ->Fill(ptTrig ,          correctConeSumClusterPhi-1);
 729     fhFractionClusterOutConePhiTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumClusterPhi-1);
 730
 731     if(coneptsumCluster!=0)
 732     {
 733             coneptsumClusterSubPhiNorm = coneptsumClusterSubPhi/coneptsumCluster;
 734         coneptsumClusterSubEtaNorm = coneptsumClusterSubEta/coneptsumCluster;
 735     }
 736
 737     fhConeSumPtSubvsConeSumPtTotPhiCluster    ->Fill(coneptsumCluster,coneptsumClusterSubPhi);
 738     fhConeSumPtSubNormvsConeSumPtTotPhiCluster->Fill(coneptsumCluster,coneptsumClusterSubPhiNorm);
 739     fhConeSumPtSubvsConeSumPtTotEtaCluster    ->Fill(coneptsumCluster,coneptsumClusterSubEta);
 740     fhConeSumPtSubNormvsConeSumPtTotEtaCluster->Fill(coneptsumCluster,coneptsumClusterSubEtaNorm);
 741
 742     // ----------- //
 743     // EMCal Cells //
 744     // ----------- //
 745
 746     // Sum the pT in the phi or eta band for clusters or tracks
 747     CalculateCaloCellUEBand(pCandidate,etaUEptsumCell   ,phiUEptsumCell   );
 748
 749     // Move to AliIsolationCut the calculation not the histograms??
 750
 751     //Careful here if EMCal limits changed .. 2010 (4 SM) to 2011-12 (10 SM), for the moment consider 100 deg in phi
 752     Float_t emcEtaSize = 0.7*2; // TO FIX
 753     Float_t emcPhiSize = TMath::DegToRad()*100.; // TO FIX
 754
 755     if(((2*conesize*emcPhiSize)-coneA)!=0)phiUEptsumCellNorm = phiUEptsumCell*(coneA / ((2*conesize*emcPhiSize)-coneA));
 756     if(((2*conesize*emcEtaSize)-coneA)!=0)etaUEptsumCellNorm = etaUEptsumCell*(coneA / ((2*conesize*emcEtaSize)-coneA));
 757
 758     // Need to correct coneptsumCluster by the fraction of the cone out of the calorimeter cut acceptance!
 759
 760     Float_t correctConeSumCellEta = 1;
 761     if(TMath::Abs(etaTrig)+conesize > emcEtaSize/2.)
 762     {
 763       Float_t excess = TMath::Abs(etaTrig) + conesize - emcEtaSize/2.;
 764       correctConeSumCellEta = GetIsolationCut()->CalculateExcessAreaFraction(excess);
 765       //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);
 766       // Need to correct phi band surface if part of the cone falls out of track cut acceptance!
 767       if(((2*(conesize-excess)*emcPhiSize)-(coneA-correctConeSumCellEta))!=0)phiUEptsumCellNorm = phiUEptsumCell*(coneA / ((2*(conesize-excess)*emcPhiSize)-(coneA-correctConeSumCellEta)));
 768     }
 769
 770     Float_t correctConeSumCellPhi = 1;
 771     //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() );
 772     if((phiTrig+conesize > 180*TMath::DegToRad()) ||
 773        (phiTrig-conesize <  80*TMath::DegToRad()))
 774     {
 775       Float_t excess = 0;
 776       if( phiTrig+conesize > 180*TMath::DegToRad() ) excess = conesize + phiTrig - 180*TMath::DegToRad() ;
 777       else                                           excess = conesize - phiTrig +  80*TMath::DegToRad() ;
 778
 779       correctConeSumCellPhi = GetIsolationCut()->CalculateExcessAreaFraction(excess);
 780       //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);
 781
 782       // Need to correct eta band surface if part of the cone falls out of track cut acceptance!
 783       if(((2*(conesize-excess)*emcEtaSize)-(coneA-correctConeSumCellPhi))!=0)etaUEptsumCellNorm = etaUEptsumCell*(coneA / ((2*(conesize-excess)*emcEtaSize)-(coneA-correctConeSumCellPhi)));
 784
 785     }
 786
 787     // In case that cone is out of eta and phi side, we are over correcting, not too often with the current cuts ...
 788     coneptsumCellSubPhi = coneptsumCell*correctConeSumCellEta*correctConeSumCellPhi - phiUEptsumCellNorm;
 789     coneptsumCellSubEta = coneptsumCell*correctConeSumCellEta*correctConeSumCellPhi - etaUEptsumCellNorm;
 790
 791     fhConeSumPtPhiUESubCell           ->Fill(ptTrig ,          coneptsumCellSubPhi);
 792     fhConeSumPtPhiUESubCellTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumCellSubPhi);
 793     fhConeSumPtEtaUESubCell           ->Fill(ptTrig ,          coneptsumCellSubEta);
 794     fhConeSumPtEtaUESubCellTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumCellSubEta);
 795
 796     fhFractionCellOutConeEta          ->Fill(ptTrig ,          correctConeSumCellEta-1);
 797     fhFractionCellOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumCellEta-1);
 798     fhFractionCellOutConePhi          ->Fill(ptTrig ,          correctConeSumCellPhi-1);
 799     fhFractionCellOutConePhiTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumCellPhi-1);
 800     if(coneptsumCell!=0)
 801     {
 802       coneptsumCellSubPhiNorm = coneptsumCellSubPhi/coneptsumCell;
 803         coneptsumCellSubEtaNorm = coneptsumCellSubEta/coneptsumCell;
 804     }
 805
 806     fhConeSumPtSubvsConeSumPtTotPhiCell    ->Fill(coneptsumCell,coneptsumCellSubPhi);
 807     fhConeSumPtSubNormvsConeSumPtTotPhiCell->Fill(coneptsumCell,coneptsumCellSubPhiNorm);
 808     fhConeSumPtSubvsConeSumPtTotEtaCell    ->Fill(coneptsumCell,coneptsumCellSubEta);
 809     fhConeSumPtSubNormvsConeSumPtTotEtaCell->Fill(coneptsumCell,coneptsumCellSubEtaNorm);
 810   }
 811
 812   if( partTypeInCone==AliIsolationCut::kNeutralAndCharged )
 813   {
 814
 815     // --------------------------- //
 816     // Tracks and clusters in cone //
 817     // --------------------------- //
 818
 819     Double_t sumPhiUESub = coneptsumClusterSubPhi + coneptsumTrackSubPhi;
 820     Double_t sumEtaUESub = coneptsumClusterSubEta + coneptsumTrackSubEta;
 821
 822     fhConeSumPtPhiUESub          ->Fill(ptTrig ,          sumPhiUESub);
 823     fhConeSumPtPhiUESubTrigEtaPhi->Fill(etaTrig, phiTrig, sumPhiUESub);
 824     fhConeSumPtEtaUESub          ->Fill(ptTrig ,          sumEtaUESub);
 825     fhConeSumPtEtaUESubTrigEtaPhi->Fill(etaTrig, phiTrig, sumEtaUESub);
 826
 827     fhEtaBandClustervsTrack    ->Fill(etaUEptsumCluster    ,etaUEptsumTrack    );
 828     fhPhiBandClustervsTrack    ->Fill(phiUEptsumCluster    ,phiUEptsumTrack    );
 829     fhEtaBandNormClustervsTrack->Fill(etaUEptsumClusterNorm,etaUEptsumTrackNorm);
 830     fhPhiBandNormClustervsTrack->Fill(phiUEptsumClusterNorm,phiUEptsumTrackNorm);
 831
 832     fhConeSumPtEtaUESubClustervsTrack->Fill(coneptsumClusterSubEta,coneptsumTrackSubEta);
 833     fhConeSumPtPhiUESubClustervsTrack->Fill(coneptsumClusterSubPhi,coneptsumTrackSubPhi);
 834
 835     // ------------------------ //
 836     // Tracks and cells in cone //
 837     // ------------------------ //
 838
 839     Double_t sumPhiUESubTrackCell = coneptsumCellSubPhi + coneptsumTrackSubPhi;
 840     Double_t sumEtaUESubTrackCell = coneptsumCellSubEta + coneptsumTrackSubEta;
 841
 842     fhConeSumPtPhiUESubTrackCell          ->Fill(ptTrig ,          sumPhiUESubTrackCell);
 843     fhConeSumPtPhiUESubTrackCellTrigEtaPhi->Fill(etaTrig, phiTrig, sumPhiUESubTrackCell);
 844     fhConeSumPtEtaUESubTrackCell          ->Fill(ptTrig ,          sumEtaUESubTrackCell);
 845     fhConeSumPtEtaUESubTrackCellTrigEtaPhi->Fill(etaTrig, phiTrig, sumEtaUESubTrackCell);
 846
 847     fhEtaBandCellvsTrack    ->Fill(etaUEptsumCell    ,etaUEptsumTrack    );
 848     fhPhiBandCellvsTrack    ->Fill(phiUEptsumCell    ,phiUEptsumTrack    );
 849     fhEtaBandNormCellvsTrack->Fill(etaUEptsumCellNorm,etaUEptsumTrackNorm);
 850     fhPhiBandNormCellvsTrack->Fill(phiUEptsumCellNorm,phiUEptsumTrackNorm);
 851
 852     fhConeSumPtEtaUESubCellvsTrack->Fill(coneptsumCellSubEta,coneptsumTrackSubEta);
 853     fhConeSumPtPhiUESubCellvsTrack->Fill(coneptsumCellSubPhi,coneptsumTrackSubPhi);
 854
 855
 856   }
 857 }
 858
 859
 860 //__________________________________________________________________________________________________
 861 void AliAnaParticleIsolation::CalculateCaloSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
 862                                                         Float_t & coneptsumCluster)
 863 {
 864   // Get the cluster pT or sum of pT in isolation cone
 865
 866   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
 867
 868   //Recover reference arrays with clusters and tracks
 869   TObjArray * refclusters = aodParticle->GetObjArray(GetAODObjArrayName()+"Clusters");
 870   if(!refclusters) return ;
 871
 872   Float_t ptTrig = aodParticle->Pt();
 873
 874   //Get vertex for cluster momentum calculation
 875   Double_t vertex[] = {0,0,0} ; //vertex ;
 876   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 877     GetReader()->GetVertex(vertex);
 878
 879   TLorentzVector mom ;
 880   for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
 881   {
 882     AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
 883     calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
 884
 885     fhPtInCone       ->Fill(ptTrig, mom.Pt());
 886     fhPtClusterInCone->Fill(ptTrig, mom.Pt());
 887
 888     if(fFillPileUpHistograms)
 889     {
 890       if(GetReader()->IsPileUpFromSPD())               fhPtInConePileUp[0]->Fill(ptTrig,mom.Pt());
 891       if(GetReader()->IsPileUpFromEMCal())             fhPtInConePileUp[1]->Fill(ptTrig,mom.Pt());
 892       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtInConePileUp[2]->Fill(ptTrig,mom.Pt());
 893       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtInConePileUp[3]->Fill(ptTrig,mom.Pt());
 894       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtInConePileUp[4]->Fill(ptTrig,mom.Pt());
 895       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtInConePileUp[5]->Fill(ptTrig,mom.Pt());
 896       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(ptTrig,mom.Pt());
 897     }
 898
 899     fhPtInConeCent->Fill(GetEventCentrality(),mom.Pt());
 900     coneptsumCluster+=mom.Pt();
 901   }
 902
 903   fhConeSumPtCluster   ->Fill(ptTrig,     coneptsumCluster);
 904 }
 905
 906 //__________________________________________________________________________________________________
 907 void AliAnaParticleIsolation::CalculateCaloCellSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
 908                                                             Float_t & coneptsumCell)
 909 {
 910   // Get the cell amplityde or sum of amplitudes in isolation cone
 911   // Mising: Remove signal cells in cone in case the trigger is a cluster!
 912
 913   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
 914
 915   Float_t conesize = GetIsolationCut()->GetConeSize();
 916
 917   Float_t  ptTrig  = aodParticle->Pt();
 918   Float_t  phiTrig = aodParticle->Phi();
 919   if(phiTrig<0) phiTrig += TMath::TwoPi();
 920   Float_t  etaTrig = aodParticle->Eta();
 921
 922   if(aodParticle->GetDetector()=="EMCAL")
 923   {
 924     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
 925     Int_t absId = -999;
 926
 927     if (eGeom->GetAbsCellIdFromEtaPhi(etaTrig,phiTrig,absId))
 928     {
 929       if(!eGeom->CheckAbsCellId(absId)) return ;
 930
 931       // Get absolute (col,row) of trigger particle
 932       Int_t nSupMod = eGeom->GetSuperModuleNumber(absId);
 933       Int_t nModule = -1;
 934       Int_t imEta=-1, imPhi=-1;
 935       Int_t ieta =-1, iphi =-1;
 936
 937       if (eGeom->GetCellIndex(absId,nSupMod,nModule,imPhi,imEta))
 938       {
 939         Int_t iEta=-1, iPhi=-1;
 940         eGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,imPhi,imEta,iphi,ieta);
 941
 942         Int_t colTrig = iEta;
 943         if (nSupMod % 2) colTrig = AliEMCALGeoParams::fgkEMCALCols + iEta ;
 944         Int_t rowTrig = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
 945
 946         Int_t sqrSize = int(conesize/0.0143);
 947
 948         AliVCaloCells * cells = GetEMCALCells();
 949
 950         // Loop on cells in cone
 951         for(Int_t irow = rowTrig-sqrSize; irow < rowTrig+sqrSize; irow++)
 952         {
 953           for(Int_t icol = colTrig-sqrSize; icol < colTrig+sqrSize; icol++)
 954           {
 955             Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
 956      if(inSector==5) continue;
 957
 958                   Int_t inSupMod = -1;
 959             Int_t icolLoc  = -1;
 960             if(icol < AliEMCALGeoParams::fgkEMCALCols)
 961             {
 962               inSupMod = 2*inSector + 1;
 963               icolLoc  = icol;
 964             }
 965             else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
 966             {
 967               inSupMod = 2*inSector;
 968               icolLoc  = icol-AliEMCALGeoParams::fgkEMCALCols;
 969             }
 970
 971             Int_t irowLoc  = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;
 972
 973             Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
 974             if(!eGeom->CheckAbsCellId(iabsId)) continue;
 975
 976             fhPtCellInCone->Fill(ptTrig, cells->GetCellAmplitude(iabsId));
 977             coneptsumCell += cells->GetCellAmplitude(iabsId);
 978           }
 979         }
 980       }
 981     }
 982   }
 983
 984   fhConeSumPtCell->Fill(ptTrig,coneptsumCell);
 985
 986 }
 987
 988 //___________________________________________________________________________________________________
 989 void AliAnaParticleIsolation::CalculateTrackSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
 990                                                          Float_t & coneptsumTrack)
 991 {
 992   // Get the track pT or sum of pT in isolation cone
 993
 994   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyNeutral ) return ;
 995
 996   //Recover reference arrays with clusters and tracks
 997   TObjArray * reftracks   = aodParticle->GetObjArray(GetAODObjArrayName()+"Tracks");
 998   if(!reftracks) return ;
 999
1000   Float_t  ptTrig = aodParticle->Pt();
1001   Double_t bz     = GetReader()->GetInputEvent()->GetMagneticField();
1002
1003   for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
1004   {
1005     AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
1006     Float_t pTtrack = track->Pt();
1007
1008     fhPtInCone     ->Fill(ptTrig,pTtrack);
1009     fhPtTrackInCone->Fill(ptTrig,pTtrack);
1010
1011     if(fFillPileUpHistograms)
1012     {
1013       ULong_t status = track->GetStatus();
1014       Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1015       //Double32_t tof = track->GetTOFsignal()*1e-3;
1016       Int_t trackBC = track->GetTOFBunchCrossing(bz);
1017
1018       if     ( okTOF && trackBC!=0 ) fhPtTrackInConeOtherBC->Fill(ptTrig,pTtrack);
1019       else if( okTOF && trackBC==0 ) fhPtTrackInConeBC0    ->Fill(ptTrig,pTtrack);
1020
1021       Int_t vtxBC = GetReader()->GetVertexBC();
1022       if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA) fhPtTrackInConeVtxBC0->Fill(ptTrig,pTtrack);
1023
1024       if(GetReader()->IsPileUpFromSPD())             { fhPtInConePileUp[0]->Fill(ptTrig,pTtrack);
1025       if(okTOF && trackBC!=0 )                         fhPtTrackInConeOtherBCPileUpSPD->Fill(ptTrig,pTtrack);
1026       if(okTOF && trackBC==0 )                         fhPtTrackInConeBC0PileUpSPD    ->Fill(ptTrig,pTtrack); }
1027       if(GetReader()->IsPileUpFromEMCal())             fhPtInConePileUp[1]->Fill(ptTrig,pTtrack);
1028       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtInConePileUp[2]->Fill(ptTrig,pTtrack);
1029       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtInConePileUp[3]->Fill(ptTrig,pTtrack);
1030       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtInConePileUp[4]->Fill(ptTrig,pTtrack);
1031       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtInConePileUp[5]->Fill(ptTrig,pTtrack);
1032       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(ptTrig,pTtrack);
1033     }
1034
1035     fhPtInConeCent->Fill(GetEventCentrality(),pTtrack);
1036     coneptsumTrack+=pTtrack;
1037   }
1038
1039   fhConeSumPtTrack->Fill(ptTrig, coneptsumTrack);
1040
1041 }
1042
1043 //_________________________________________________________________
1044 void AliAnaParticleIsolation::FillPileUpHistograms(Int_t clusterID)
1045 {
1046   // Fill some histograms to understand pile-up
1047   if(!fFillPileUpHistograms) return;
1048
1049   if(clusterID < 0 )
1050   {
1051     printf("AliAnaParticleIsolation::FillPileUpHistograms(), ID of cluster = %d, not possible! ", clusterID);
1052     return;
1053   }
1054
1055   Int_t iclus = -1;
1056   TObjArray* clusters = 0x0;
1057   if     (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
1058   else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
1059
1060   Float_t energy = 0;
1061   Float_t time   = -1000;
1062
1063   if(clusters)
1064   {
1065     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
1066     energy = cluster->E();
1067     time   = cluster->GetTOF()*1e9;
1068   }
1069
1070   //printf("E %f, time %f\n",energy,time);
1071   AliVEvent * event = GetReader()->GetInputEvent();
1072
1073   fhTimeENoCut->Fill(energy,time);
1074   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
1075   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
1076
1077   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
1078
1079   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
1080   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
1081
1082   // N pile up vertices
1083   Int_t nVerticesSPD    = -1;
1084   Int_t nVerticesTracks = -1;
1085
1086   if      (esdEv)
1087   {
1088     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
1089     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
1090
1091   }//ESD
1092   else if (aodEv)
1093   {
1094     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
1095     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
1096   }//AOD
1097
1098   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
1099   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
1100
1101   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
1102   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
1103
1104   Int_t ncont = -1;
1105   Float_t z1 = -1, z2 = -1;
1106   Float_t diamZ = -1;
1107   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
1108   {
1109     if      (esdEv)
1110     {
1111       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
1112       ncont=pv->GetNContributors();
1113       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
1114       z2 = pv->GetZ();
1115       diamZ = esdEv->GetDiamondZ();
1116     }//ESD
1117     else if (aodEv)
1118     {
1119       AliAODVertex *pv=aodEv->GetVertex(iVert);
1120       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
1121       ncont=pv->GetNContributors();
1122       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
1123       z2=pv->GetZ();
1124       diamZ = aodEv->GetDiamondZ();
1125     }// AOD
1126
1127     Double_t distZ  = TMath::Abs(z2-z1);
1128     diamZ  = TMath::Abs(z2-diamZ);
1129
1130     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
1131     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
1132     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
1133
1134   }// loop
1135 }
1136
1137 //_____________________________________________________________________________________________________________________
1138 void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(AliAODPWG4ParticleCorrelation  *pCandidate,
1139                                                                             AliCaloPID * pid)
1140 {
1141   // Fill Track matching and Shower Shape control histograms
1142   if(!fFillTMHisto &&  !fFillSSHisto) return;
1143
1144   Int_t  clusterID = pCandidate->GetCaloLabel(0) ;
1145   Int_t  nMaxima   = pCandidate->GetFiducialArea(); // bad name, just place holder for the moment
1146   Int_t  mcTag     = pCandidate->GetTag() ;
1147   Bool_t isolated  = pCandidate->IsIsolated();
1148
1149   if(clusterID < 0 )
1150   {
1151     printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(), ID of cluster = %d, not possible! \n", clusterID);
1152     return;
1153   }
1154
1155   Int_t iclus = -1;
1156   TObjArray* clusters = 0x0;
1157   if     (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
1158   else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
1159
1160   if(clusters)
1161   {
1162     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
1163     Float_t energy = cluster->E();
1164
1165     if(fFillSSHisto)
1166     {
1167       fhELambda0[isolated]->Fill(energy, cluster->GetM02() );
1168       fhELambda1[isolated]->Fill(energy, cluster->GetM20() );
1169
1170       if(IsDataMC())
1171       {
1172         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||
1173                 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) fhELambda0MCPhoton    [isolated]->Fill(energy, cluster->GetM02());
1174         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))           fhELambda0MCPi0       [isolated]->Fill(energy, cluster->GetM02());
1175         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      fhELambda0MCPi0Decay  [isolated]->Fill(energy, cluster->GetM02());
1176         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      fhELambda0MCEtaDecay  [isolated]->Fill(energy, cluster->GetM02());
1177         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    fhELambda0MCOtherDecay[isolated]->Fill(energy, cluster->GetM02());
1178
1179         //        else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(energy, cluster->GetM02());
1180         else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))     fhELambda0MCHadron    [isolated]->Fill(energy, cluster->GetM02());
1181
1182       }
1183
1184       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5) // TO DO: CHANGE FOR 2012
1185       {
1186         fhELambda0TRD[isolated]->Fill(energy, cluster->GetM02() );
1187         fhELambda1TRD[isolated]->Fill(energy, cluster->GetM20() );
1188       }
1189
1190       fhNLocMax[isolated]->Fill(energy,nMaxima);
1191       if     (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax1[isolated]->Fill(energy,cluster->GetM20()); }
1192       else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax2[isolated]->Fill(energy,cluster->GetM20()); }
1193       else                { fhELambda0LocMaxN[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMaxN[isolated]->Fill(energy,cluster->GetM20()); }
1194
1195       if(isolated==0)
1196       {
1197         //Analyse non-isolated events
1198         Int_t   n         = 0;
1199         Int_t   nfrac     = 0;
1200         Bool_t  iso       = kFALSE ;
1201         Float_t coneptsum = 0 ;
1202
1203         TObjArray * plNe  = pCandidate->GetObjArray(GetAODObjArrayName()+"Clusters");
1204         TObjArray * plCTS = pCandidate->GetObjArray(GetAODObjArrayName()+"Tracks");
1205
1206         GetIsolationCut()->SetPtThresholdMax(1.);
1207         GetIsolationCut()->MakeIsolationCut(plCTS,   plNe,
1208                                             GetReader(), pid,
1209                                             kFALSE, pCandidate, "",
1210                                             n,nfrac,coneptsum, iso);
1211
1212         if (!iso) fhELambda0SSBkg->Fill(energy, cluster->GetM02());
1213
1214         if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
1215       }
1216
1217       GetIsolationCut()->SetPtThresholdMax(10000.);
1218
1219     } // SS histo fill
1220
1221
1222     if(fFillTMHisto)
1223     {
1224       Float_t dZ  = cluster->GetTrackDz();
1225       Float_t dR  = cluster->GetTrackDx();
1226
1227       if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1228       {
1229         dR = 2000., dZ = 2000.;
1230         GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1231       }
1232
1233       //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
1234       if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
1235       {
1236         fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
1237         fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
1238         if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
1239       }
1240
1241       // Check dEdx and E/p of matched clusters
1242
1243       if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1244       {
1245
1246         AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1247
1248         if(track)
1249         {
1250           Float_t dEdx = track->GetTPCsignal();
1251           fhdEdx[isolated]->Fill(cluster->E(), dEdx);
1252
1253           Float_t eOverp = cluster->E()/track->P();
1254           fhEOverP[isolated]->Fill(cluster->E(),  eOverp);
1255         }
1256         //else
1257         //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1258
1259
1260         if(IsDataMC())
1261         {
1262           if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
1263           {
1264             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
1265                        GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
1266             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
1267             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
1268             else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
1269
1270           }
1271           else
1272           {
1273             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
1274                        GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
1275             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
1276             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
1277             else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
1278           }
1279
1280         }  // MC
1281
1282       } // match window
1283
1284     }// TM histos fill
1285
1286   } // clusters array available
1287
1288 }
1289
1290 //______________________________________________________
1291 TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
1292 {
1293   //Save parameters used for analysis
1294   TString parList ; //this will be list of parameters used for this analysis.
1295   const Int_t buffersize = 255;
1296   char onePar[buffersize] ;
1297
1298   snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
1299   parList+=onePar ;
1300   snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1301   parList+=onePar ;
1302   snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
1303   parList+=onePar ;
1304   snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
1305   parList+=onePar ;
1306   snprintf(onePar, buffersize,"fFillTMHisto=%d (Flag for track matching histograms) \n",fFillTMHisto) ;
1307   parList+=onePar ;
1308   snprintf(onePar, buffersize,"fFillSSHisto=%d (Flag for shower shape histograms) \n",fFillSSHisto) ;
1309   parList+=onePar ;
1310
1311   if(fMakeSeveralIC)
1312   {
1313     snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
1314     parList+=onePar ;
1315     snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
1316     parList+=onePar ;
1317
1318     for(Int_t icone = 0; icone < fNCones ; icone++)
1319     {
1320       snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
1321       parList+=onePar ;
1322     }
1323     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1324     {
1325       snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
1326       parList+=onePar ;
1327     }
1328     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1329     {
1330       snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
1331       parList+=onePar ;
1332     }
1333     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1334     {
1335       snprintf(onePar, buffersize,"fSumPtThresholds[%d]=%1.2f (isolation sum pt threshold) \n",ipt, fSumPtThresholds[ipt]) ;
1336       parList+=onePar ;
1337     }
1338   }
1339
1340   //Get parameters set in base class.
1341   parList += GetBaseParametersList() ;
1342
1343   //Get parameters set in IC class.
1344   if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
1345
1346   return new TObjString(parList) ;
1347
1348 }
1349
1350 //________________________________________________________
1351 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
1352 {
1353   // Create histograms to be saved in output file and
1354   // store them in outputContainer
1355   TList * outputContainer = new TList() ;
1356   outputContainer->SetName("IsolatedParticleHistos") ;
1357
1358   Int_t   nptbins  = GetHistogramRanges()->GetHistoPtBins();
1359   Int_t   nphibins = GetHistogramRanges()->GetHistoPhiBins();
1360   Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
1361   Float_t ptmax    = GetHistogramRanges()->GetHistoPtMax();
1362   Float_t phimax   = GetHistogramRanges()->GetHistoPhiMax();
1363   Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
1364   Float_t ptmin    = GetHistogramRanges()->GetHistoPtMin();
1365   Float_t phimin   = GetHistogramRanges()->GetHistoPhiMin();
1366   Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
1367   Int_t   ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();
1368   Float_t ssmax    = GetHistogramRanges()->GetHistoShowerShapeMax();
1369   Float_t ssmin    = GetHistogramRanges()->GetHistoShowerShapeMin();
1370   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
1371   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
1372   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
1373
1374   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1375   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1376   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1377   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1378   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1379   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1380
1381   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
1382   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
1383   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1384   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1385   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
1386   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1387
1388   Int_t   nptsumbins    = fHistoNPtSumBins;
1389   Float_t ptsummax      = fHistoPtSumMax;
1390   Float_t ptsummin      = fHistoPtSumMin;
1391   Int_t   nptinconebins = fHistoNPtInConeBins;
1392   Float_t ptinconemax   = fHistoPtInConeMax;
1393   Float_t ptinconemin   = fHistoPtInConeMin;
1394
1395   Float_t ptthre = GetIsolationCut()->GetPtThreshold();
1396   Float_t ptfrac = GetIsolationCut()->GetPtFraction();
1397   Float_t r      = GetIsolationCut()->GetConeSize();
1398
1399   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1400
1401   if(!fMakeSeveralIC)
1402   {
1403     TString hName [] = {"NoIso",""};
1404     TString hTitle[] = {"Not isolated"  ,"isolated"};
1405
1406     fhEIso   = new TH1F("hE",
1407                         Form("Number of isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1408                         nptbins,ptmin,ptmax);
1409     fhEIso->SetYTitle("dN / dE");
1410     fhEIso->SetXTitle("E (GeV/c)");
1411     outputContainer->Add(fhEIso) ;
1412
1413     fhPtIso  = new TH1F("hPt",
1414                         Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1415                         nptbins,ptmin,ptmax);
1416     fhPtIso->SetYTitle("dN / p_{T}");
1417     fhPtIso->SetXTitle("p_{T} (GeV/c)");
1418     outputContainer->Add(fhPtIso) ;
1419
1420     fhPtCentralityIso  = new TH2F("hPtCentrality","centrality vs p_{T} for isolated particles",nptbins,ptmin,ptmax, 100,0,100);
1421     fhPtCentralityIso->SetYTitle("centrality");
1422     fhPtCentralityIso->SetXTitle("p_{T}(GeV/c)");
1423     outputContainer->Add(fhPtCentralityIso) ;
1424
1425     fhPtEventPlaneIso  = new TH2F("hPtEventPlane","event plane angle vs p_{T} for isolated particles",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1426     fhPtEventPlaneIso->SetYTitle("Event plane angle (rad)");
1427     fhPtEventPlaneIso->SetXTitle("p_{T} (GeV/c)");
1428     outputContainer->Add(fhPtEventPlaneIso) ;
1429
1430
1431     fhPtNLocMaxIso  = new TH2F("hPtNLocMax",
1432                                Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f vs NLM, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1433                                nptbins,ptmin,ptmax,10,0,10);
1434     fhPtNLocMaxIso->SetYTitle("NLM");
1435     fhPtNLocMaxIso->SetXTitle("p_{T} (GeV/c)");
1436     outputContainer->Add(fhPtNLocMaxIso) ;
1437
1438     fhPhiIso  = new TH2F("hPhi",
1439                          Form("Number of isolated particles vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1440                          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1441     fhPhiIso->SetYTitle("#phi");
1442     fhPhiIso->SetXTitle("p_{T} (GeV/c)");
1443     outputContainer->Add(fhPhiIso) ;
1444
1445     fhEtaIso  = new TH2F("hEta",
1446                          Form("Number of isolated particles vs #eta for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1447                          nptbins,ptmin,ptmax,netabins,etamin,etamax);
1448     fhEtaIso->SetYTitle("#eta");
1449     fhEtaIso->SetXTitle("p_{T} (GeV/c)");
1450     outputContainer->Add(fhEtaIso) ;
1451
1452     fhEtaPhiIso  = new TH2F("hEtaPhiIso",
1453                             Form("Number of isolated particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1454                             netabins,etamin,etamax,nphibins,phimin,phimax);
1455     fhEtaPhiIso->SetXTitle("#eta");
1456     fhEtaPhiIso->SetYTitle("#phi");
1457     outputContainer->Add(fhEtaPhiIso) ;
1458
1459     fhPtDecayIso  = new TH1F("hPtDecayIso",
1460                              Form("Number of isolated #pi^{0} decay particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1461                              nptbins,ptmin,ptmax);
1462     fhPtDecayIso->SetYTitle("N");
1463     fhPtDecayIso->SetXTitle("p_{T}(GeV/c)");
1464     outputContainer->Add(fhPtDecayIso) ;
1465
1466     fhEtaPhiDecayIso  = new TH2F("hEtaPhiDecayIso",
1467                                  Form("Number of isolated Pi0 decay particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1468                                  netabins,etamin,etamax,nphibins,phimin,phimax);
1469     fhEtaPhiDecayIso->SetXTitle("#eta");
1470     fhEtaPhiDecayIso->SetYTitle("#phi");
1471     outputContainer->Add(fhEtaPhiDecayIso) ;
1472
1473     fhConeSumPt  = new TH2F("hConePtSum",
1474                             Form("Track and Cluster #Sigma p_{T} in isolation cone for R = %2.2f",r),
1475                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1476     fhConeSumPt->SetYTitle("#Sigma p_{T}");
1477     fhConeSumPt->SetXTitle("p_{T, trigger} (GeV/c)");
1478     outputContainer->Add(fhConeSumPt) ;
1479
1480     fhConeSumPtTrigEtaPhi  = new TH2F("hConePtSumTrigEtaPhi",
1481                             Form("Trigger #eta vs #phi, #Sigma p_{T} in isolation cone for R = %2.2f",r),
1482                             netabins,etamin,etamax,nphibins,phimin,phimax);
1483     fhConeSumPtTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1484     fhConeSumPtTrigEtaPhi->SetXTitle("#eta_{trigger}");
1485     fhConeSumPtTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1486     outputContainer->Add(fhConeSumPtTrigEtaPhi) ;
1487
1488     fhPtInCone  = new TH2F("hPtInCone",
1489                            Form("p_{T} of clusters and tracks in isolation cone for R = %2.2f",r),
1490                            nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1491     fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
1492     fhPtInCone->SetXTitle("p_{T} (GeV/c)");
1493     outputContainer->Add(fhPtInCone) ;
1494
1495     fhPtInConeCent  = new TH2F("hPtInConeCent",
1496                                Form("p_{T} in isolation cone for R = %2.2f",r),
1497                                100,0,100,nptinconebins,ptinconemin,ptinconemax);
1498     fhPtInConeCent->SetYTitle("p_{T in cone} (GeV/c)");
1499     fhPtInConeCent->SetXTitle("centrality");
1500     outputContainer->Add(fhPtInConeCent) ;
1501
1502     // Cluster only histograms
1503     if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyCharged)
1504     {
1505       fhConeSumPtCluster  = new TH2F("hConePtSumCluster",
1506                                      Form("Cluster #Sigma p_{T} in isolation cone for R = %2.2f",r),
1507                                      nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1508       fhConeSumPtCluster->SetYTitle("#Sigma p_{T}");
1509       fhConeSumPtCluster->SetXTitle("p_{T, trigger} (GeV/c)");
1510       outputContainer->Add(fhConeSumPtCluster) ;
1511
1512       fhConeSumPtCell  = new TH2F("hConePtSumCell",
1513                                   Form("Cell #Sigma p_{T} in isolation cone for R = %2.2f",r),
1514                                   nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1515       fhConeSumPtCell->SetYTitle("#Sigma p_{T}");
1516       fhConeSumPtCell->SetXTitle("p_{T, trigger} (GeV/c)");
1517       outputContainer->Add(fhConeSumPtCell) ;
1518
1519       fhConeSumPtEtaBandUECluster  = new TH2F("hConePtSumEtaBandUECluster",
1520                                               "#Sigma cluster p_{T} in UE Eta Band",
1521                                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1522       fhConeSumPtEtaBandUECluster->SetYTitle("#Sigma p_{T}");
1523       fhConeSumPtEtaBandUECluster->SetXTitle("p_{T, trigger} (GeV/c)");
1524       outputContainer->Add(fhConeSumPtEtaBandUECluster) ;
1525
1526       fhConeSumPtPhiBandUECluster  = new TH2F("hConePtSumPhiBandUECluster",
1527                                               "#Sigma cluster p_{T} UE Phi Band",
1528                                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1529       fhConeSumPtPhiBandUECluster->SetYTitle("#Sigma p_{T}");
1530       fhConeSumPtPhiBandUECluster->SetXTitle("p_{T, trigger} (GeV/c)");
1531       outputContainer->Add(fhConeSumPtPhiBandUECluster) ;
1532
1533       fhConeSumPtEtaBandUEClusterTrigEtaPhi  = new TH2F("hConePtSumEtaBandUEClusterTrigEtaPhi",
1534                                                         "Trigger #eta vs #phi, #Sigma cluster p_{T} in UE Eta Band",
1535                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
1536       fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1537       fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1538       fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1539       outputContainer->Add(fhConeSumPtEtaBandUEClusterTrigEtaPhi) ;
1540
1541       fhConeSumPtPhiBandUEClusterTrigEtaPhi  = new TH2F("hConePtSumPhiBandUEClusterTrigEtaPhi",
1542                                                         "Trigger #eta vs #phi, #Sigma cluster p_{T} UE Phi Band",
1543                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
1544       fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1545       fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1546       fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1547       outputContainer->Add(fhConeSumPtPhiBandUEClusterTrigEtaPhi) ;
1548
1549
1550       fhConeSumPtEtaBandUECell  = new TH2F("hConePtSumEtaBandUECell",
1551                                            "#Sigma cell p_{T} in UE Eta Band",
1552                                            nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1553       fhConeSumPtEtaBandUECell->SetYTitle("#Sigma p_{T}");
1554       fhConeSumPtEtaBandUECell->SetXTitle("p_{T, trigger} (GeV/c)");
1555       outputContainer->Add(fhConeSumPtEtaBandUECell) ;
1556
1557       fhConeSumPtPhiBandUECell  = new TH2F("hConePtSumPhiBandUECell",
1558                                            "#Sigma cell p_{T} UE Phi Band",
1559                                            nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1560       fhConeSumPtPhiBandUECell->SetYTitle("#Sigma p_{T}");
1561       fhConeSumPtPhiBandUECell->SetXTitle("p_{T, trigger} (GeV/c)");
1562       outputContainer->Add(fhConeSumPtPhiBandUECell) ;
1563
1564       fhConeSumPtEtaBandUECellTrigEtaPhi  = new TH2F("hConePtSumEtaBandUECellTrigEtaPhi",
1565                                                      "Trigger #eta vs #phi, #Sigma cell p_{T} in UE Eta Band",
1566                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
1567       fhConeSumPtEtaBandUECellTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1568       fhConeSumPtEtaBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1569       fhConeSumPtEtaBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1570       outputContainer->Add(fhConeSumPtEtaBandUECellTrigEtaPhi) ;
1571
1572       fhConeSumPtPhiBandUECellTrigEtaPhi  = new TH2F("hConePtSumPhiBandUECellTrigEtaPhi",
1573                                                      "Trigger #eta vs #phi, #Sigma cell p_{T} UE Phi Band",
1574                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
1575       fhConeSumPtPhiBandUECellTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1576       fhConeSumPtPhiBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1577       fhConeSumPtPhiBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1578       outputContainer->Add(fhConeSumPtPhiBandUECellTrigEtaPhi) ;
1579
1580
1581       fhPtClusterInCone  = new TH2F("hPtClusterInCone",
1582                                     Form("p_{T} of clusters in isolation cone for R = %2.2f",r),
1583                                     nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1584       fhPtClusterInCone->SetYTitle("p_{T in cone} (GeV/c)");
1585       fhPtClusterInCone->SetXTitle("p_{T} (GeV/c)");
1586       outputContainer->Add(fhPtClusterInCone) ;
1587
1588       fhEtaBandCluster  = new TH2F("hEtaBandCluster",
1589                                    Form("#eta vs #phi of clusters in #eta band isolation cone for R = %2.2f",r),
1590                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
1591       fhEtaBandCluster->SetXTitle("#eta");
1592       fhEtaBandCluster->SetYTitle("#phi");
1593       outputContainer->Add(fhEtaBandCluster) ;
1594
1595       fhPhiBandCluster  = new TH2F("hPhiBandCluster",
1596                                    Form("#eta vs #phi of clusters in #phi band isolation cone for R = %2.2f",r),
1597                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
1598       fhPhiBandCluster->SetXTitle("#eta");
1599       fhPhiBandCluster->SetYTitle("#phi");
1600       outputContainer->Add(fhPhiBandCluster) ;
1601
1602         fhEtaPhiInConeCluster= new TH2F("hEtaPhiInConeCluster",
1603                                    Form("#eta vs #phi of clusters in cone for R = %2.2f",r),
1604                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
1605       fhEtaPhiInConeCluster->SetXTitle("#eta");
1606       fhEtaPhiInConeCluster->SetYTitle("#phi");
1607       outputContainer->Add(fhEtaPhiInConeCluster) ;
1608
1609         fhEtaPhiCluster= new TH2F("hEtaPhiCluster",
1610                                    Form("#eta vs #phi of all clusters"),
1611                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
1612       fhEtaPhiCluster->SetXTitle("#eta");
1613       fhEtaPhiCluster->SetYTitle("#phi");
1614       outputContainer->Add(fhEtaPhiCluster) ;
1615
1616         fhEtaPhiInConeTrack= new TH2F("hEtaPhiInConeTrack",
1617                                    Form("#eta vs #phi of Tracks in cone for R = %2.2f",r),
1618                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
1619       fhEtaPhiInConeTrack->SetXTitle("#eta");
1620       fhEtaPhiInConeTrack->SetYTitle("#phi");
1621       outputContainer->Add(fhEtaPhiInConeTrack) ;
1622
1623         fhEtaPhiTrack= new TH2F("hEtaPhiTrack",
1624                                    Form("#eta vs #phi of all Tracks"),
1625                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
1626       fhEtaPhiTrack->SetXTitle("#eta");
1627       fhEtaPhiTrack->SetYTitle("#phi");
1628       outputContainer->Add(fhEtaPhiTrack) ;
1629
1630     fhPtCellInCone  = new TH2F("hPtCellInCone",
1631                                  Form("p_{T} of cells in isolation cone for R = %2.2f",r),
1632                                  nptbins,ptmin,ptmax,1000,0,50);
1633       fhPtCellInCone->SetYTitle("p_{T in cone} (GeV/c)");
1634       fhPtCellInCone->SetXTitle("p_{T} (GeV/c)");
1635       outputContainer->Add(fhPtCellInCone) ;
1636
1637       fhEtaBandCell  = new TH2F("hEtaBandCell",
1638                                 Form("#col vs #row of cells in #eta band isolation cone for R = %2.2f",r),
1639                                 96,0,95,128,0,127);
1640       fhEtaBandCell->SetXTitle("#col");
1641       fhEtaBandCell->SetYTitle("#row");
1642       outputContainer->Add(fhEtaBandCell) ;
1643
1644       fhPhiBandCell  = new TH2F("hPhiBandCell",
1645                                 Form("#col vs #row of cells in #phi band isolation cone for R = %2.2f",r),
1646                                 96,0,95,128,0,127);
1647       fhPhiBandCell->SetXTitle("#col");
1648       fhPhiBandCell->SetYTitle("#row");
1649       outputContainer->Add(fhPhiBandCell) ;
1650
1651       fhConeSumPtEtaUESubCluster  = new TH2F("hConeSumPtEtaUESubCluster",
1652                                              Form("Clusters #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
1653                                              nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1654       fhConeSumPtEtaUESubCluster->SetYTitle("#Sigma p_{T}");
1655       fhConeSumPtEtaUESubCluster->SetXTitle("p_{T} (GeV/c)");
1656       outputContainer->Add(fhConeSumPtEtaUESubCluster) ;
1657
1658       fhConeSumPtPhiUESubCluster  = new TH2F("hConeSumPtPhiUESubCluster",
1659                                              Form("Clusters #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
1660                                              nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1661       fhConeSumPtPhiUESubCluster->SetYTitle("#Sigma p_{T}");
1662       fhConeSumPtPhiUESubCluster->SetXTitle("p_{T} (GeV/c)");
1663       outputContainer->Add(fhConeSumPtPhiUESubCluster) ;
1664
1665       fhConeSumPtEtaUESubClusterTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubClusterTrigEtaPhi",
1666                                                        Form("Trigger #eta vs #phi, Clusters #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
1667                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
1668       fhConeSumPtEtaUESubClusterTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1669       fhConeSumPtEtaUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1670       fhConeSumPtEtaUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1671       outputContainer->Add(fhConeSumPtEtaUESubClusterTrigEtaPhi) ;
1672
1673       fhConeSumPtPhiUESubClusterTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubClusterTrigEtaPhi",
1674                                                        Form("Trigger #eta vs #phi, Clusters #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
1675                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
1676       fhConeSumPtPhiUESubClusterTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1677       fhConeSumPtPhiUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1678       fhConeSumPtPhiUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1679       outputContainer->Add(fhConeSumPtPhiUESubClusterTrigEtaPhi) ;
1680
1681
1682       fhConeSumPtEtaUESubCell  = new TH2F("hConeSumPtEtaUESubCell",
1683                                           Form("Cells #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
1684                                           nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1685       fhConeSumPtEtaUESubCell->SetYTitle("#Sigma p_{T}");
1686       fhConeSumPtEtaUESubCell->SetXTitle("p_{T} (GeV/c)");
1687       outputContainer->Add(fhConeSumPtEtaUESubCell) ;
1688
1689       fhConeSumPtPhiUESubCell  = new TH2F("hConeSumPtPhiUESubCell",
1690                                           Form("Cells #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
1691                                           nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1692       fhConeSumPtPhiUESubCell->SetYTitle("#Sigma p_{T}");
1693       fhConeSumPtPhiUESubCell->SetXTitle("p_{T} (GeV/c)");
1694       outputContainer->Add(fhConeSumPtPhiUESubCell) ;
1695
1696       fhConeSumPtEtaUESubCellTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubCellTrigEtaPhi",
1697                                                     Form("Trigger #eta vs #phi, Cells #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
1698                                                     netabins,etamin,etamax,nphibins,phimin,phimax);
1699       fhConeSumPtEtaUESubCellTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1700       fhConeSumPtEtaUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1701       fhConeSumPtEtaUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1702       outputContainer->Add(fhConeSumPtEtaUESubCellTrigEtaPhi) ;
1703
1704       fhConeSumPtPhiUESubCellTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubCellTrigEtaPhi",
1705                                                     Form("Trigger #eta vs #phi, Cells #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
1706                                                     netabins,etamin,etamax,nphibins,phimin,phimax);
1707       fhConeSumPtPhiUESubCellTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1708       fhConeSumPtPhiUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1709       fhConeSumPtPhiUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1710       outputContainer->Add(fhConeSumPtPhiUESubCellTrigEtaPhi) ;
1711
1712
1713       fhFractionClusterOutConeEta  = new TH2F("hFractionClusterOutConeEta",
1714                                               Form("Fraction of the isolation cone R = %2.2f, out of clusters #eta acceptance",r),
1715                                               nptbins,ptmin,ptmax,100,0,1);
1716       fhFractionClusterOutConeEta->SetYTitle("fraction");
1717       fhFractionClusterOutConeEta->SetXTitle("p_{T,trigger} (GeV/c)");
1718       outputContainer->Add(fhFractionClusterOutConeEta) ;
1719
1720       fhFractionClusterOutConeEtaTrigEtaPhi  = new TH2F("hFractionClusterOutConeEtaTrigEtaPhi",
1721                                                         Form("Fraction of the isolation cone R = %2.2f, out of clusters #eta acceptance, in trigger #eta-#phi ",r),
1722                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
1723       fhFractionClusterOutConeEtaTrigEtaPhi->SetZTitle("fraction");
1724       fhFractionClusterOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
1725       fhFractionClusterOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1726       outputContainer->Add(fhFractionClusterOutConeEtaTrigEtaPhi) ;
1727
1728       fhFractionClusterOutConePhi  = new TH2F("hFractionClusterOutConePhi",
1729                                               Form("Fraction of the isolation cone R = %2.2f, out of clusters #phi acceptance",r),
1730                                               nptbins,ptmin,ptmax,100,0,1);
1731       fhFractionClusterOutConePhi->SetYTitle("fraction");
1732       fhFractionClusterOutConePhi->SetXTitle("p_{T,trigger} (GeV/c)");
1733       outputContainer->Add(fhFractionClusterOutConePhi) ;
1734
1735       fhFractionClusterOutConePhiTrigEtaPhi  = new TH2F("hFractionClusterOutConePhiTrigEtaPhi",
1736                                                         Form("Fraction of the isolation cone R = %2.2f, out of clusters #phi acceptance, in trigger #eta-#phi ",r),
1737                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
1738       fhFractionClusterOutConePhiTrigEtaPhi->SetZTitle("fraction");
1739       fhFractionClusterOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
1740       fhFractionClusterOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1741       outputContainer->Add(fhFractionClusterOutConePhiTrigEtaPhi) ;
1742
1743       fhConeSumPtSubvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCluster",
1744                                                         Form("#Sigma p_{T} in cone after bkg sub from #phi band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1745                                                         nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1746       fhConeSumPtSubvsConeSumPtTotPhiCluster->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1747       fhConeSumPtSubvsConeSumPtTotPhiCluster->SetYTitle("#Sigma p_{T, sub} (GeV/c)");
1748       outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCluster);
1749
1750       fhConeSumPtSubNormvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCluster",
1751                                                             Form("#Sigma p_{T, norm} in cone after bkg sub from #phi band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1752                                                             nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1753       fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1754       fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetYTitle("#Sigma p_{T, sub norm} (GeV/c)");
1755       outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCluster);
1756
1757       fhConeSumPtSubvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCluster",
1758                                                         Form("#Sigma p_{T} in cone after bkg sub from #eta band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1759                                                         nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1760       fhConeSumPtSubvsConeSumPtTotEtaCluster->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1761       fhConeSumPtSubvsConeSumPtTotEtaCluster->SetYTitle("#Sigma p_{T, sub} (GeV/c)");
1762       outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCluster);
1763
1764       fhConeSumPtSubNormvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCluster",
1765                                                             Form("#Sigma p_{T, norm} in cone after bkg sub from #eta band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1766                                                             nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1767       fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1768       fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetYTitle("#Sigma p_{T, sub norm} (GeV/c)");
1769       outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCluster);
1770
1771
1772       fhFractionCellOutConeEta  = new TH2F("hFractionCellOutConeEta",
1773                                            Form("Fraction of the isolation cone R = %2.2f, out of cells #eta acceptance",r),
1774                                            nptbins,ptmin,ptmax,100,0,1);
1775       fhFractionCellOutConeEta->SetYTitle("fraction");
1776       fhFractionCellOutConeEta->SetXTitle("p_{T,trigger} (GeV/c)");
1777       outputContainer->Add(fhFractionCellOutConeEta) ;
1778
1779       fhFractionCellOutConeEtaTrigEtaPhi  = new TH2F("hFractionCellOutConeEtaTrigEtaPhi",
1780                                                      Form("Fraction of the isolation cone R = %2.2f, out of cells #eta acceptance, in trigger #eta-#phi ",r),
1781                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
1782       fhFractionCellOutConeEtaTrigEtaPhi->SetZTitle("fraction");
1783       fhFractionCellOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
1784       fhFractionCellOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1785       outputContainer->Add(fhFractionCellOutConeEtaTrigEtaPhi) ;
1786
1787       fhFractionCellOutConePhi  = new TH2F("hFractionCellOutConePhi",
1788                                            Form("Fraction of the isolation cone R = %2.2f, out of cells #phi acceptance",r),
1789                                            nptbins,ptmin,ptmax,100,0,1);
1790       fhFractionCellOutConePhi->SetYTitle("fraction");
1791       fhFractionCellOutConePhi->SetXTitle("p_{T,trigger} (GeV/c)");
1792       outputContainer->Add(fhFractionCellOutConePhi) ;
1793
1794       fhFractionCellOutConePhiTrigEtaPhi  = new TH2F("hFractionCellOutConePhiTrigEtaPhi",
1795                                                      Form("Fraction of the isolation cone R = %2.2f, out of cells #phi acceptance, in trigger #eta-#phi ",r),
1796                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
1797       fhFractionCellOutConePhiTrigEtaPhi->SetZTitle("fraction");
1798       fhFractionCellOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
1799       fhFractionCellOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1800       outputContainer->Add(fhFractionCellOutConePhiTrigEtaPhi) ;
1801
1802
1803       fhConeSumPtSubvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCell",
1804                                                      Form("#Sigma p_{T} in cone after bkg sub from #phi band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1805                                                      nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1806       fhConeSumPtSubvsConeSumPtTotPhiCell->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1807       fhConeSumPtSubvsConeSumPtTotPhiCell->SetYTitle("#Sigma p_{T, sub} (GeV/c)");
1808       outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCell);
1809
1810       fhConeSumPtSubNormvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCell",
1811                                                          Form("#Sigma p_{T, norm} in cone after bkg sub from #phi band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1812                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1813       fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1814       fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetYTitle("#Sigma p_{T, sub norm} (GeV/c)");
1815       outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCell);
1816
1817       fhConeSumPtSubvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCell",
1818                                                      Form("#Sigma p_{T} in cone after bkg sub from #eta band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1819                                                      nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1820       fhConeSumPtSubvsConeSumPtTotEtaCell->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1821       fhConeSumPtSubvsConeSumPtTotEtaCell->SetYTitle("#Sigma p_{T, sub} (GeV/c)");
1822       outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCell);
1823
1824       fhConeSumPtSubNormvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCell",
1825                                                          Form("#Sigma p_{T, norm} in cone after bkg sub from #eta band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1826                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1827       fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1828       fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetYTitle("#Sigma p_{T, sub norm} (GeV/c)");
1829       outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCell);
1830
1831       fhConeSumPtVSUETracksEtaBand  = new TH2F("hConeSumPtVSUETracksEtaBand",
1832                                                          Form("#Sigma p_{T} in cone versus #Sigma p_{T} in eta band for tracks (before normalization), R=%2.2f",r),
1833                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
1834       fhConeSumPtVSUETracksEtaBand->SetXTitle("#Sigma p_{T} cone (GeV/c)");
1835       fhConeSumPtVSUETracksEtaBand->SetYTitle("#Sigma p_{T} UE (GeV/c)");
1836       outputContainer->Add(fhConeSumPtVSUETracksEtaBand);
1837
1838       fhConeSumPtVSUETracksPhiBand  = new TH2F("hConeSumPtVSUETracksPhiBand",
1839                                                          Form("#Sigma p_{T} in cone versus #Sigma p_{T} in phi band for tracks (before normalization), R=%2.2f",r),
1840                                                          nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
1841       fhConeSumPtVSUETracksPhiBand->SetXTitle("#Sigma p_{T} cone (GeV/c)");
1842       fhConeSumPtVSUETracksPhiBand->SetYTitle("#Sigma p_{T} UE (GeV/c)");
1843       outputContainer->Add(fhConeSumPtVSUETracksPhiBand);
1844
1845       fhConeSumPtVSUEClusterEtaBand  = new TH2F("hConeSumPtVSUEClusterEtaBand",
1846                                                          Form("#Sigma p_{T} in cone versus #Sigma p_{T} in eta band for cluster (before normalization), R=%2.2f",r),
1847                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
1848       fhConeSumPtVSUEClusterEtaBand->SetXTitle("#Sigma p_{T} cone (GeV/c)");
1849       fhConeSumPtVSUEClusterEtaBand->SetYTitle("#Sigma p_{T} UE (GeV/c)");
1850       outputContainer->Add(fhConeSumPtVSUEClusterEtaBand);
1851
1852       fhConeSumPtVSUEClusterPhiBand  = new TH2F("hConeSumPtVSUEClusterPhiBand",
1853                                                          Form("#Sigma p_{T} in cone versus #Sigma p_{T} in phi band for cluster (before normalization), R=%2.2f",r),
1854                                                          nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
1855       fhConeSumPtVSUEClusterPhiBand->SetXTitle("#Sigma p_{T} cone (GeV/c)");
1856       fhConeSumPtVSUEClusterPhiBand->SetYTitle("#Sigma p_{T} UE (GeV/c)");
1857       outputContainer->Add(fhConeSumPtVSUEClusterPhiBand);
1858
1859     }
1860
1861     // Track only histograms
1862     if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyNeutral)
1863     {
1864       fhConeSumPtTrack  = new TH2F("hConePtSumTrack",
1865                                    Form("Track #Sigma p_{T} in isolation cone for R = %2.2f",r),
1866                                    nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1867       fhConeSumPtTrack->SetYTitle("#Sigma p_{T}");
1868       fhConeSumPtTrack->SetXTitle("p_{T, trigger} (GeV/c)");
1869       outputContainer->Add(fhConeSumPtTrack) ;
1870
1871
1872       fhConeSumPtEtaBandUETrack  = new TH2F("hConePtSumEtaBandUETrack",
1873                                             "#Sigma track p_{T} in UE Eta Band",
1874                                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1875       fhConeSumPtEtaBandUETrack->SetYTitle("#Sigma p_{T}");
1876       fhConeSumPtEtaBandUETrack->SetXTitle("p_{T, trigger} (GeV/c)");
1877       outputContainer->Add(fhConeSumPtEtaBandUETrack) ;
1878
1879       fhConeSumPtPhiBandUETrack  = new TH2F("hConePtSumPhiBandUETrack",
1880                                             "#Sigma track p_{T} in UE Phi Band",
1881                                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax*8);
1882       fhConeSumPtPhiBandUETrack->SetYTitle("#Sigma p_{T}");
1883       fhConeSumPtPhiBandUETrack->SetXTitle("p_{T, trigger} (GeV/c)");
1884       outputContainer->Add(fhConeSumPtPhiBandUETrack) ;
1885
1886
1887       fhConeSumPtEtaBandUETrackTrigEtaPhi  = new TH2F("hConePtSumEtaBandUETrackTrigEtaPhi",
1888                                                       "Trigger #eta vs #phi, #Sigma track p_{T} in UE Eta Band",
1889                                                       netabins,etamin,etamax,nphibins,phimin,phimax);
1890       fhConeSumPtEtaBandUETrackTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1891       fhConeSumPtEtaBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
1892       fhConeSumPtEtaBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1893       outputContainer->Add(fhConeSumPtEtaBandUETrackTrigEtaPhi) ;
1894
1895       fhConeSumPtPhiBandUETrackTrigEtaPhi  = new TH2F("hConePtSumPhiBandUETrackTrigEtaPhi",
1896                                                       "Trigger #eta vs #phi, #Sigma track p_{T} in UE Phi Band",
1897                                                       netabins,etamin,etamax,nphibins,phimin,phimax);
1898       fhConeSumPtPhiBandUETrackTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1899       fhConeSumPtPhiBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
1900       fhConeSumPtPhiBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1901       outputContainer->Add(fhConeSumPtPhiBandUETrackTrigEtaPhi) ;
1902
1903
1904       fhPtTrackInCone  = new TH2F("hPtTrackInCone",
1905                                   Form("p_{T} of tracks in isolation cone for R = %2.2f",r),
1906                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1907       fhPtTrackInCone->SetYTitle("p_{T in cone} (GeV/c)");
1908       fhPtTrackInCone->SetXTitle("p_{T} (GeV/c)");
1909       outputContainer->Add(fhPtTrackInCone) ;
1910
1911
1912       fhEtaBandTrack  = new TH2F("hEtaBandTrack",
1913                                  Form("#eta vs #phi of tracks in #eta band isolation cone for R = %2.2f",r),
1914                                  netabins,-1,1,nphibins,0,TMath::TwoPi());
1915       fhEtaBandTrack->SetXTitle("#eta");
1916       fhEtaBandTrack->SetYTitle("#phi");
1917       outputContainer->Add(fhEtaBandTrack) ;
1918
1919       fhPhiBandTrack  = new TH2F("hPhiBandTrack",
1920                                  Form("#eta vs #phi of tracks in #phi band isolation cone for R = %2.2f",r),
1921                                  netabins,-1,1,nphibins,0,TMath::TwoPi());
1922       fhPhiBandTrack->SetXTitle("#eta");
1923       fhPhiBandTrack->SetYTitle("#phi");
1924       outputContainer->Add(fhPhiBandTrack) ;
1925
1926
1927       fhConeSumPtEtaUESubTrack  = new TH2F("hConeSumPtEtaUESubTrack",
1928                                            Form("Tracks #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
1929                                            nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1930       fhConeSumPtEtaUESubTrack->SetYTitle("#Sigma p_{T}");
1931       fhConeSumPtEtaUESubTrack->SetXTitle("p_{T} (GeV/c)");
1932       outputContainer->Add(fhConeSumPtEtaUESubTrack) ;
1933
1934       fhConeSumPtPhiUESubTrack  = new TH2F("hConeSumPtPhiUESubTrack",
1935                                            Form("Tracks #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
1936                                            nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1937       fhConeSumPtPhiUESubTrack->SetYTitle("#Sigma p_{T}");
1938       fhConeSumPtPhiUESubTrack->SetXTitle("p_{T} (GeV/c)");
1939       outputContainer->Add(fhConeSumPtPhiUESubTrack) ;
1940
1941       fhConeSumPtEtaUESubTrackTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrackTrigEtaPhi",
1942                                                      Form("Trigger #eta vs #phi, Tracks #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
1943                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
1944       fhConeSumPtEtaUESubTrackTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1945       fhConeSumPtEtaUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
1946       fhConeSumPtEtaUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1947       outputContainer->Add(fhConeSumPtEtaUESubTrackTrigEtaPhi) ;
1948
1949       fhConeSumPtPhiUESubTrackTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrackTrigEtaPhi",
1950                                                      Form("Trigger #eta vs #phi, Tracks #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
1951                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
1952       fhConeSumPtPhiUESubTrackTrigEtaPhi->SetZTitle("#Sigma p_{T}");
1953       fhConeSumPtPhiUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
1954       fhConeSumPtPhiUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1955       outputContainer->Add(fhConeSumPtPhiUESubTrackTrigEtaPhi) ;
1956
1957       fhFractionTrackOutConeEta  = new TH2F("hFractionTrackOutConeEta",
1958                                             Form("Fraction of the isolation cone R = %2.2f, out of tracks #eta acceptance",r),
1959                                             nptbins,ptmin,ptmax,100,0,1);
1960       fhFractionTrackOutConeEta->SetYTitle("fraction");
1961       fhFractionTrackOutConeEta->SetXTitle("p_{T,trigger} (GeV/c)");
1962       outputContainer->Add(fhFractionTrackOutConeEta) ;
1963
1964       fhFractionTrackOutConeEtaTrigEtaPhi  = new TH2F("hFractionTrackOutConeEtaTrigEtaPhi",
1965                                                       Form("Fraction of the isolation cone R = %2.2f, out of tracks #eta acceptance, in trigger #eta-#phi ",r),
1966                                                       netabins,etamin,etamax,nphibins,phimin,phimax);
1967       fhFractionTrackOutConeEtaTrigEtaPhi->SetZTitle("fraction");
1968       fhFractionTrackOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
1969       fhFractionTrackOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1970       outputContainer->Add(fhFractionTrackOutConeEtaTrigEtaPhi) ;
1971
1972       fhConeSumPtSubvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubvsConeSumPtTotPhiTrack",
1973                                                       Form("#Sigma p_{T} in cone after bkg sub from #phi band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1974                                                       nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1975       fhConeSumPtSubvsConeSumPtTotPhiTrack->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1976       fhConeSumPtSubvsConeSumPtTotPhiTrack->SetYTitle("#Sigma p_{T, sub} (GeV/c)");
1977       outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiTrack);
1978
1979       fhConeSumPtSubNormvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiTrack",
1980                                                           Form("#Sigma p_{T, norm} in cone after bkg sub from #phi band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1981                                                           nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1982       fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1983       fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetYTitle("#Sigma p_{T, sub norm} (GeV/c)");
1984       outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiTrack);
1985
1986       fhConeSumPtSubvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubvsConeSumPtTotEtaTrack",
1987                                                       Form("#Sigma p_{T} in cone after bkg sub from #eta band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1988                                                       nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1989       fhConeSumPtSubvsConeSumPtTotEtaTrack->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1990       fhConeSumPtSubvsConeSumPtTotEtaTrack->SetYTitle("#Sigma p_{T, sub} (GeV/c)");
1991       outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaTrack);
1992
1993       fhConeSumPtSubNormvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaTrack",
1994                                                           Form("#Sigma p_{T, norm} in cone after bkg sub from #eta band vs #Sigma p_{T} in cone before bkg sub, R=%2.2f",r),
1995                                                           nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
1996       fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetXTitle("#Sigma p_{T, tot} (GeV/c)");
1997       fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetYTitle("#Sigma p_{T, sub norm} (GeV/c)");
1998       outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaTrack);
1999
2000       // UE in perpendicular cone
2001       fhPerpConeSumPt  = new TH2F("hPerpConePtSum",
2002                                   Form("#Sigma p_{T} in isolation cone at #pm 45 degree phi from trigger particle, R = %2.2f",r),
2003                                   nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2004       fhPerpConeSumPt->SetYTitle("#Sigma p_{T}");
2005       fhPerpConeSumPt->SetXTitle("p_{T} (GeV/c)");
2006       outputContainer->Add(fhPerpConeSumPt) ;
2007
2008       fhPtInPerpCone  = new TH2F("hPtInPerpCone",
2009                                  Form("p_{T} in isolation cone at #pm 45 degree phi from trigger particle, R = %2.2f",r),
2010                                  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2011       fhPtInPerpCone->SetYTitle("p_{T in cone} (GeV/c)");
2012       fhPtInPerpCone->SetXTitle("p_{T} (GeV/c)");
2013       outputContainer->Add(fhPtInPerpCone) ;
2014     }
2015
2016     if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged)
2017     {
2018       fhConeSumPtEtaUESub  = new TH2F("hConeSumPtEtaUESub",
2019                                       Form("#Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
2020                                       nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2021       fhConeSumPtEtaUESub->SetYTitle("#Sigma p_{T}");
2022       fhConeSumPtEtaUESub->SetXTitle("p_{T} (GeV/c)");
2023       outputContainer->Add(fhConeSumPtEtaUESub) ;
2024
2025       fhConeSumPtPhiUESub  = new TH2F("hConeSumPtPhiUESub",
2026                                       Form("#Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
2027                                       nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2028       fhConeSumPtPhiUESub->SetYTitle("#Sigma p_{T}");
2029       fhConeSumPtPhiUESub->SetXTitle("p_{T} (GeV/c)");
2030       outputContainer->Add(fhConeSumPtPhiUESub) ;
2031
2032       fhConeSumPtEtaUESubTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrigEtaPhi",
2033                                                 Form("Trigger #eta vs #phi, #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
2034                                                 netabins,etamin,etamax,nphibins,phimin,phimax);
2035       fhConeSumPtEtaUESubTrigEtaPhi->SetZTitle("#Sigma p_{T}");
2036       fhConeSumPtEtaUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2037       fhConeSumPtEtaUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2038       outputContainer->Add(fhConeSumPtEtaUESubTrigEtaPhi) ;
2039
2040       fhConeSumPtPhiUESubTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrigEtaPhi",
2041                                                 Form("Trigger #eta vs #phi, #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
2042                                                 netabins,etamin,etamax,nphibins,phimin,phimax);
2043       fhConeSumPtPhiUESubTrigEtaPhi->SetZTitle("#Sigma p_{T}");
2044       fhConeSumPtPhiUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2045       fhConeSumPtPhiUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2046       outputContainer->Add(fhConeSumPtPhiUESubTrigEtaPhi) ;
2047
2048       fhConeSumPtClustervsTrack   = new TH2F("hConePtSumClustervsTrack",
2049                                              Form("Track vs Cluster #Sigma p_{T} in isolation cone for R = %2.2f",r),
2050                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2051       fhConeSumPtClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2052       fhConeSumPtClustervsTrack->SetYTitle("#Sigma p_{T} track");
2053       outputContainer->Add(fhConeSumPtClustervsTrack) ;
2054
2055       fhConeSumPtEtaUESubClustervsTrack   = new TH2F("hConePtSumEtaUESubClustervsTrack",
2056                                                      Form("Track vs Cluster #Sigma p_{T} UE sub eta band in isolation cone for R = %2.2f",r),
2057                                                      2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2058       fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2059       fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma p_{T} track");
2060       outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2061
2062       fhConeSumPtPhiUESubClustervsTrack   = new TH2F("hConePhiUESubPtSumClustervsTrack",
2063                                                      Form("Track vs Cluster #Sigma p_{T} UE sub phi band in isolation cone for R = %2.2f",r),
2064                                                      2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2065       fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2066       fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma p_{T} track");
2067       outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2068
2069       fhEtaBandClustervsTrack   = new TH2F("hEtaBandClustervsTrack",
2070                                            "Track vs Cluster #Sigma p_{T} in Eta band in isolation cone for R = %2.2f",
2071                                            nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2072       fhEtaBandClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2073       fhEtaBandClustervsTrack->SetYTitle("#Sigma p_{T} track");
2074       outputContainer->Add(fhEtaBandClustervsTrack) ;
2075
2076       fhPhiBandClustervsTrack   = new TH2F("hPhiBandClustervsTrack",
2077                                            "Track vs Cluster #Sigma p_{T} in Phi band in isolation cone for R = %2.2f",
2078                                            nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2079       fhPhiBandClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2080       fhPhiBandClustervsTrack->SetYTitle("#Sigma p_{T} track");
2081       outputContainer->Add(fhPhiBandClustervsTrack) ;
2082
2083       fhEtaBandNormClustervsTrack   = new TH2F("hEtaBandNormClustervsTrack",
2084                                                "Track vs Cluster #Sigma p_{T} in Eta band in isolation cone for R = %2.2f",
2085                                                nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2086       fhEtaBandNormClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2087       fhEtaBandNormClustervsTrack->SetYTitle("#Sigma p_{T} track");
2088       outputContainer->Add(fhEtaBandNormClustervsTrack) ;
2089
2090       fhPhiBandNormClustervsTrack   = new TH2F("hPhiBandNormClustervsTrack",
2091                                                "Track vs Cluster #Sigma p_{T} in Phi band in isolation cone for R = %2.2f",
2092                                                nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2093       fhPhiBandNormClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2094       fhPhiBandNormClustervsTrack->SetYTitle("#Sigma p_{T} track");
2095       outputContainer->Add(fhPhiBandNormClustervsTrack) ;
2096
2097
2098       fhConeSumPtCellTrack = new TH2F("hConePtSumCellTrack",
2099                                       Form("Track and Cell #Sigma p_{T} in isolation cone for R = %2.2f",r),
2100                                       nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2101       fhConeSumPtCellTrack->SetYTitle("#Sigma p_{T}");
2102       fhConeSumPtCellTrack->SetXTitle("p_{T, trigger} (GeV/c)");
2103       outputContainer->Add(fhConeSumPtCellTrack) ;
2104
2105       fhConeSumPtCellTrackTrigEtaPhi  = new TH2F("hConePtSumCellTrackTrigEtaPhi",
2106                                                  Form("Trigger #eta vs #phi, #Sigma p_{T} in isolation cone for R = %2.2f",r),
2107                                                  netabins,etamin,etamax,nphibins,phimin,phimax);
2108       fhConeSumPtCellTrackTrigEtaPhi->SetZTitle("#Sigma p_{T}");
2109       fhConeSumPtCellTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2110       fhConeSumPtCellTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2111       outputContainer->Add(fhConeSumPtCellTrackTrigEtaPhi) ;
2112
2113
2114       fhConeSumPtEtaUESubClustervsTrack   = new TH2F("hConePtSumEtaUESubClustervsTrack",
2115                                                      Form("Track vs Cluster #Sigma p_{T} UE sub eta band in isolation cone for R = %2.2f",r),
2116                                                      2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2117       fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2118       fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma p_{T} track");
2119       outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2120
2121       fhConeSumPtPhiUESubClustervsTrack   = new TH2F("hConePhiUESubPtSumClustervsTrack",
2122                                                      Form("Track vs Cluster #Sigma p_{T} UE sub phi band in isolation cone for R = %2.2f",r),
2123                                                      2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2124       fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma p_{T} cluster");
2125       fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma p_{T} track");
2126       outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2127
2128       fhConeSumPtCellvsTrack   = new TH2F("hConePtSumCellvsTrack",
2129                                              Form("Track vs cell #Sigma p_{T} in isolation cone for R = %2.2f",r),
2130                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2131       fhConeSumPtCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2132       fhConeSumPtCellvsTrack->SetYTitle("#Sigma p_{T} track");
2133       outputContainer->Add(fhConeSumPtCellvsTrack) ;
2134
2135       fhConeSumPtEtaUESubCellvsTrack   = new TH2F("hConePtSumEtaUESubCellvsTrack",
2136                                                   Form("Track vs Cell #Sigma p_{T} UE sub eta band in isolation cone for R = %2.2f",r),
2137                                                   2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2138       fhConeSumPtEtaUESubCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2139       fhConeSumPtEtaUESubCellvsTrack->SetYTitle("#Sigma p_{T} track");
2140       outputContainer->Add(fhConeSumPtEtaUESubCellvsTrack) ;
2141
2142       fhConeSumPtPhiUESubCellvsTrack   = new TH2F("hConePhiUESubPtSumCellvsTrack",
2143                                                   Form("Track vs Cell #Sigma p_{T} UE sub phi band in isolation cone for R = %2.2f",r),
2144                                                   2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2145       fhConeSumPtPhiUESubCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2146       fhConeSumPtPhiUESubCellvsTrack->SetYTitle("#Sigma p_{T} track");
2147       outputContainer->Add(fhConeSumPtPhiUESubCellvsTrack) ;
2148
2149       fhEtaBandCellvsTrack   = new TH2F("hEtaBandCellvsTrack",
2150                                         "Track vs Cell #Sigma p_{T} in Eta band in isolation cone for R = %2.2f",
2151                                         nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2152       fhEtaBandCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2153       fhEtaBandCellvsTrack->SetYTitle("#Sigma p_{T} track");
2154       outputContainer->Add(fhEtaBandCellvsTrack) ;
2155
2156       fhPhiBandCellvsTrack   = new TH2F("hPhiBandCellvsTrack",
2157                                         "Track vs Cell #Sigma p_{T} in Phi band in isolation cone for R = %2.2f",
2158                                         nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2159       fhPhiBandCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2160       fhPhiBandCellvsTrack->SetYTitle("#Sigma p_{T} track");
2161       outputContainer->Add(fhPhiBandCellvsTrack) ;
2162
2163       fhEtaBandNormCellvsTrack   = new TH2F("hEtaBandNormCellvsTrack",
2164                                             "Track vs Cell #Sigma p_{T} in Eta band in isolation cone for R = %2.2f",
2165                                             nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2166       fhEtaBandNormCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2167       fhEtaBandNormCellvsTrack->SetYTitle("#Sigma p_{T} track");
2168       outputContainer->Add(fhEtaBandNormCellvsTrack) ;
2169
2170       fhPhiBandNormCellvsTrack   = new TH2F("hPhiBandNormCellvsTrack",
2171                                             "Track vs Cell #Sigma p_{T} in Phi band in isolation cone for R = %2.2f",
2172                                             nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2173       fhPhiBandNormCellvsTrack->SetXTitle("#Sigma p_{T} cell");
2174       fhPhiBandNormCellvsTrack->SetYTitle("#Sigma p_{T} track");
2175       outputContainer->Add(fhPhiBandNormCellvsTrack) ;
2176
2177       fhConeSumPtEtaUESubTrackCell  = new TH2F("hConeSumPtEtaUESubTrackCell",
2178                                            Form("Tracks #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
2179                                            nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2180       fhConeSumPtEtaUESubTrackCell->SetYTitle("#Sigma p_{T}");
2181       fhConeSumPtEtaUESubTrackCell->SetXTitle("p_{T} (GeV/c)");
2182       outputContainer->Add(fhConeSumPtEtaUESubTrackCell) ;
2183
2184       fhConeSumPtPhiUESubTrackCell  = new TH2F("hConeSumPtPhiUESubTrackCell",
2185                                            Form("Tracks #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
2186                                            nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2187       fhConeSumPtPhiUESubTrackCell->SetYTitle("#Sigma p_{T}");
2188       fhConeSumPtPhiUESubTrackCell->SetXTitle("p_{T} (GeV/c)");
2189       outputContainer->Add(fhConeSumPtPhiUESubTrackCell) ;
2190
2191       fhConeSumPtEtaUESubTrackCellTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrackCellTrigEtaPhi",
2192                                                      Form("Trigger #eta vs #phi, Tracks #Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
2193                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
2194       fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma p_{T}");
2195       fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2196       fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2197       outputContainer->Add(fhConeSumPtEtaUESubTrackCellTrigEtaPhi) ;
2198
2199       fhConeSumPtPhiUESubTrackCellTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrackCellTrigEtaPhi",
2200                                                      Form("Trigger #eta vs #phi, Tracks #Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
2201                                                      netabins,etamin,etamax,nphibins,phimin,phimax);
2202       fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma p_{T}");
2203       fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2204       fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2205       outputContainer->Add(fhConeSumPtPhiUESubTrackCellTrigEtaPhi) ;
2206
2207     }
2208
2209     if(fFillSSHisto)
2210     {
2211         fhELambda0SSBkg  = new TH2F
2212       ("hELambda0SSBkg","Non isolated clusters : E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2213       fhELambda0SSBkg->SetYTitle("#lambda_{0}^{2}");
2214       fhELambda0SSBkg->SetXTitle("E (GeV)");
2215       outputContainer->Add(fhELambda0SSBkg) ;
2216     }
2217
2218     for(Int_t iso = 0; iso < 2; iso++)
2219     {
2220       if(fFillTMHisto)
2221       {
2222         fhTrackMatchedDEta[iso]  = new TH2F
2223         (Form("hTrackMatchedDEta%s",hName[iso].Data()),
2224          Form("%s - d#eta of cluster-track vs cluster energy for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
2225          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2226         fhTrackMatchedDEta[iso]->SetYTitle("d#eta");
2227         fhTrackMatchedDEta[iso]->SetXTitle("E_{cluster} (GeV)");
2228
2229         fhTrackMatchedDPhi[iso]  = new TH2F
2230         (Form("hTrackMatchedDPhi%s",hName[iso].Data()),
2231          Form("%s - d#phi of cluster-track vs cluster energy for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
2232          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2233         fhTrackMatchedDPhi[iso]->SetYTitle("d#phi (rad)");
2234         fhTrackMatchedDPhi[iso]->SetXTitle("E_{cluster} (GeV)");
2235
2236         fhTrackMatchedDEtaDPhi[iso]  = new TH2F
2237         (Form("hTrackMatchedDEtaDPhi%s",hName[iso].Data()),
2238          Form("%s - d#eta vs d#phi of cluster-track for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
2239          nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2240         fhTrackMatchedDEtaDPhi[iso]->SetYTitle("d#phi (rad)");
2241         fhTrackMatchedDEtaDPhi[iso]->SetXTitle("d#eta");
2242
2243         outputContainer->Add(fhTrackMatchedDEta[iso]) ;
2244         outputContainer->Add(fhTrackMatchedDPhi[iso]) ;
2245         outputContainer->Add(fhTrackMatchedDEtaDPhi[iso]) ;
2246
2247         fhdEdx[iso]  = new TH2F
2248         (Form("hdEdx%s",hName[iso].Data()),
2249          Form("%s - Matched track <dE/dx> vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
2250          nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2251         fhdEdx[iso]->SetXTitle("E (GeV)");
2252         fhdEdx[iso]->SetYTitle("<dE/dx>");
2253         outputContainer->Add(fhdEdx[iso]);
2254
2255         fhEOverP[iso]  = new TH2F
2256         (Form("hEOverP%s",hName[iso].Data()),
2257          Form("%s - Matched track E/p vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
2258          nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2259         fhEOverP[iso]->SetXTitle("E (GeV)");
2260         fhEOverP[iso]->SetYTitle("E/p");
2261         outputContainer->Add(fhEOverP[iso]);
2262
2263         if(IsDataMC())
2264         {
2265           fhTrackMatchedMCParticle[iso]  = new TH2F
2266           (Form("hTrackMatchedMCParticle%s",hName[iso].Data()),
2267            Form("%s - Origin of particle vs energy vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
2268            nptbins,ptmin,ptmax,8,0,8);
2269           fhTrackMatchedMCParticle[iso]->SetXTitle("E (GeV)");
2270           //fhTrackMatchedMCParticle[iso]->SetYTitle("Particle type");
2271
2272           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(1 ,"Photon");
2273           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(2 ,"Electron");
2274           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
2275           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(4 ,"Rest");
2276           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
2277           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
2278           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
2279           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
2280
2281           outputContainer->Add(fhTrackMatchedMCParticle[iso]);
2282         }
2283       }
2284
2285       if(fFillSSHisto)
2286       {
2287         fhELambda0[iso]  = new TH2F
2288         (Form("hELambda0%s",hName[iso].Data()),
2289          Form("%s cluster : E vs #lambda_{0}",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2290         fhELambda0[iso]->SetYTitle("#lambda_{0}^{2}");
2291         fhELambda0[iso]->SetXTitle("E (GeV)");
2292         outputContainer->Add(fhELambda0[iso]) ;
2293
2294         if(IsDataMC())
2295         {
2296           fhELambda0MCPhoton[iso]  = new TH2F
2297           (Form("hELambda0%s_MCPhoton",hName[iso].Data()),
2298            Form("%s cluster : E vs #lambda_{0}: Origin is final state photon",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2299           fhELambda0MCPhoton[iso]->SetYTitle("#lambda_{0}^{2}");
2300           fhELambda0MCPhoton[iso]->SetXTitle("E (GeV)");
2301           outputContainer->Add(fhELambda0MCPhoton[iso]) ;
2302
2303           fhELambda0MCPi0[iso]  = new TH2F
2304           (Form("hELambda0%s_MCPi0",hName[iso].Data()),
2305            Form("%s cluster : E vs #lambda_{0}: Origin is pi0 (2 #gamma)",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2306           fhELambda0MCPi0[iso]->SetYTitle("#lambda_{0}^{2}");
2307           fhELambda0MCPi0[iso]->SetXTitle("E (GeV)");
2308           outputContainer->Add(fhELambda0MCPi0[iso]) ;
2309
2310           fhELambda0MCPi0Decay[iso]  = new TH2F
2311           (Form("hELambda0%s_MCPi0Decay",hName[iso].Data()),
2312            Form("%s cluster : E vs #lambda_{0}: Origin is pi0 decay",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2313           fhELambda0MCPi0Decay[iso]->SetYTitle("#lambda_{0}^{2}");
2314           fhELambda0MCPi0Decay[iso]->SetXTitle("E (GeV)");
2315           outputContainer->Add(fhELambda0MCPi0Decay[iso]) ;
2316
2317           fhELambda0MCEtaDecay[iso]  = new TH2F
2318           (Form("hELambda0%s_MCEtaDecay",hName[iso].Data()),
2319            Form("%s cluster : E vs #lambda_{0}: Origin is eta decay",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2320           fhELambda0MCEtaDecay[iso]->SetYTitle("#lambda_{0}^{2}");
2321           fhELambda0MCEtaDecay[iso]->SetXTitle("E (GeV)");
2322           outputContainer->Add(fhELambda0MCEtaDecay[iso]) ;
2323
2324           fhELambda0MCOtherDecay[iso]  = new TH2F
2325           (Form("hELambda0%s_MCOtherDecay",hName[iso].Data()),
2326            Form("%s cluster : E vs #lambda_{0}: Origin is other decay",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2327           fhELambda0MCOtherDecay[iso]->SetYTitle("#lambda_{0}^{2}");
2328           fhELambda0MCOtherDecay[iso]->SetXTitle("E (GeV)");
2329           outputContainer->Add(fhELambda0MCOtherDecay[iso]) ;
2330
2331           fhELambda0MCHadron[iso]  = new TH2F
2332           (Form("hELambda0%s_MCHadron",hName[iso].Data()),
2333            Form("%s cluster : E vs #lambda_{0}: Origin is hadron",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2334           fhELambda0MCHadron[iso]->SetYTitle("#lambda_{0}^{2}");
2335           fhELambda0MCHadron[iso]->SetXTitle("E (GeV)");
2336           outputContainer->Add(fhELambda0MCHadron[iso]) ;
2337         }
2338
2339         fhELambda1[iso]  = new TH2F
2340         (Form("hELambda1%s",hName[iso].Data()),
2341          Form("%s cluster: E vs #lambda_{1}",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2342         fhELambda1[iso]->SetYTitle("#lambda_{1}^{2}");
2343         fhELambda1[iso]->SetXTitle("E (GeV)");
2344         outputContainer->Add(fhELambda1[iso]) ;
2345
2346         if(fCalorimeter=="EMCAL")
2347         {
2348           fhELambda0TRD[iso]  = new TH2F
2349           (Form("hELambda0TRD%s",hName[iso].Data()),
2350            Form("%s cluster: E vs #lambda_{0}, SM behind TRD",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2351           fhELambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
2352           fhELambda0TRD[iso]->SetXTitle("E (GeV)");
2353           outputContainer->Add(fhELambda0TRD[iso]) ;
2354
2355           fhELambda1TRD[iso]  = new TH2F
2356           (Form("hELambda1TRD%s",hName[iso].Data()),
2357            Form("%s cluster: E vs #lambda_{1}, SM behind TRD",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2358           fhELambda1TRD[iso]->SetYTitle("#lambda_{1}^{2}");
2359           fhELambda1TRD[iso]->SetXTitle("E (GeV)");
2360           outputContainer->Add(fhELambda1TRD[iso]) ;
2361         }
2362
2363         fhNLocMax[iso] = new TH2F
2364         (Form("hNLocMax%s",hName[iso].Data()),
2365          Form("%s - Number of local maxima in cluster",hTitle[iso].Data()),
2366          nptbins,ptmin,ptmax,10,0,10);
2367         fhNLocMax[iso]->SetYTitle("N maxima");
2368         fhNLocMax[iso]->SetXTitle("E (GeV)");
2369         outputContainer->Add(fhNLocMax[iso]) ;
2370
2371         fhELambda0LocMax1[iso]  = new TH2F
2372         (Form("hELambda0LocMax1%s",hName[iso].Data()),
2373          Form("%s cluster (#eta) pairs: E vs #lambda_{0}, 1 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2374         fhELambda0LocMax1[iso]->SetYTitle("#lambda_{0}^{2}");
2375         fhELambda0LocMax1[iso]->SetXTitle("E (GeV)");
2376         outputContainer->Add(fhELambda0LocMax1[iso]) ;
2377
2378         fhELambda1LocMax1[iso]  = new TH2F
2379         (Form("hELambda1LocMax1%s",hName[iso].Data()),
2380          Form("%s cluster (#eta) pairs: E vs #lambda_{1}, 1 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2381         fhELambda1LocMax1[iso]->SetYTitle("#lambda_{1}^{2}");
2382         fhELambda1LocMax1[iso]->SetXTitle("E (GeV)");
2383         outputContainer->Add(fhELambda1LocMax1[iso]) ;
2384
2385         fhELambda0LocMax2[iso]  = new TH2F
2386         (Form("hELambda0LocMax2%s",hName[iso].Data()),
2387          Form("%s cluster (#eta) pairs: E vs #lambda_{0}, 2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2388         fhELambda0LocMax2[iso]->SetYTitle("#lambda_{0}^{2}");
2389         fhELambda0LocMax2[iso]->SetXTitle("E (GeV)");
2390         outputContainer->Add(fhELambda0LocMax2[iso]) ;
2391
2392         fhELambda1LocMax2[iso]  = new TH2F
2393         (Form("hELambda1LocMax2%s",hName[iso].Data()),
2394          Form("%s cluster (#eta) pairs: E vs #lambda_{1}, 2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2395         fhELambda1LocMax2[iso]->SetYTitle("#lambda_{1}^{2}");
2396         fhELambda1LocMax2[iso]->SetXTitle("E (GeV)");
2397         outputContainer->Add(fhELambda1LocMax2[iso]) ;
2398
2399         fhELambda0LocMaxN[iso]  = new TH2F
2400         ( Form("hELambda0LocMaxN%s",hName[iso].Data()),
2401          Form("%s cluster (#eta) pairs: E vs #lambda_{0}, N>2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2402         fhELambda0LocMaxN[iso]->SetYTitle("#lambda_{0}^{2}");
2403         fhELambda0LocMaxN[iso]->SetXTitle("E (GeV)");
2404         outputContainer->Add(fhELambda0LocMaxN[iso]) ;
2405
2406         fhELambda1LocMaxN[iso]  = new TH2F
2407         (Form("hELambda1LocMaxN%s",hName[iso].Data()),
2408          Form("%s cluster (#eta) pairs: E vs #lambda_{1}, N>2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2409         fhELambda1LocMaxN[iso]->SetYTitle("#lambda_{1}^{2}");
2410         fhELambda1LocMaxN[iso]->SetXTitle("E (GeV)");
2411         outputContainer->Add(fhELambda1LocMaxN[iso]) ;
2412
2413       }
2414     } // control histograms for isolated and non isolated objects
2415
2416
2417     if(fFillPileUpHistograms)
2418     {
2419       fhPtTrackInConeOtherBC  = new TH2F("hPtTrackInConeOtherBC",
2420                                          Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC!=0",r),
2421                                          nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2422       fhPtTrackInConeOtherBC->SetYTitle("p_{T in cone} (GeV/c)");
2423       fhPtTrackInConeOtherBC->SetXTitle("p_{T} (GeV/c)");
2424       outputContainer->Add(fhPtTrackInConeOtherBC) ;
2425
2426       fhPtTrackInConeOtherBCPileUpSPD  = new TH2F("hPtTrackInConeOtherBCPileUpSPD",
2427                                                   Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC!=0, pile-up from SPD",r),
2428                                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2429       fhPtTrackInConeOtherBCPileUpSPD->SetYTitle("p_{T in cone} (GeV/c)");
2430       fhPtTrackInConeOtherBCPileUpSPD->SetXTitle("p_{T} (GeV/c)");
2431       outputContainer->Add(fhPtTrackInConeOtherBCPileUpSPD) ;
2432
2433       fhPtTrackInConeBC0  = new TH2F("hPtTrackInConeBC0",
2434                                      Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC==0",r),
2435                                      nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2436       fhPtTrackInConeBC0->SetYTitle("p_{T in cone} (GeV/c)");
2437       fhPtTrackInConeBC0->SetXTitle("p_{T} (GeV/c)");
2438       outputContainer->Add(fhPtTrackInConeBC0) ;
2439
2440       fhPtTrackInConeVtxBC0  = new TH2F("hPtTrackInConeVtxBC0",
2441                                         Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC==0",r),
2442                                         nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2443       fhPtTrackInConeVtxBC0->SetYTitle("p_{T in cone} (GeV/c)");
2444       fhPtTrackInConeVtxBC0->SetXTitle("p_{T} (GeV/c)");
2445       outputContainer->Add(fhPtTrackInConeVtxBC0) ;
2446
2447
2448       fhPtTrackInConeBC0PileUpSPD  = new TH2F("hPtTrackInConeBC0PileUpSPD",
2449                                               Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC==0, pile-up from SPD",r),
2450                                               nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2451       fhPtTrackInConeBC0PileUpSPD->SetYTitle("p_{T in cone} (GeV/c)");
2452       fhPtTrackInConeBC0PileUpSPD->SetXTitle("p_{T} (GeV/c)");
2453       outputContainer->Add(fhPtTrackInConeBC0PileUpSPD) ;
2454
2455
2456       for (Int_t i = 0; i < 7 ; i++)
2457       {
2458         fhPtInConePileUp[i]  = new TH2F(Form("hPtInConePileUp%s",pileUpName[i].Data()),
2459                                         Form("p_{T} in isolation cone for R = %2.2f, from pile-up (%s)",r,pileUpName[i].Data()),
2460                                         nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2461         fhPtInConePileUp[i]->SetYTitle("p_{T in cone} (GeV/c)");
2462         fhPtInConePileUp[i]->SetXTitle("p_{T} (GeV/c)");
2463         outputContainer->Add(fhPtInConePileUp[i]) ;
2464       }
2465     }
2466
2467     if(IsDataMC())
2468     {
2469       fhPtIsoPrompt  = new TH1F("hPtMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax);
2470       fhPtIsoPrompt->SetYTitle("N");
2471       fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
2472       outputContainer->Add(fhPtIsoPrompt) ;
2473
2474       fhPhiIsoPrompt  = new TH2F
2475       ("hPhiMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2476       fhPhiIsoPrompt->SetYTitle("#phi");
2477       fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
2478       outputContainer->Add(fhPhiIsoPrompt) ;
2479
2480       fhEtaIsoPrompt  = new TH2F
2481       ("hEtaMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2482       fhEtaIsoPrompt->SetYTitle("#eta");
2483       fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
2484       outputContainer->Add(fhEtaIsoPrompt) ;
2485
2486       fhPtIsoFragmentation  = new TH1F("hPtMCFragmentation","Number of isolated #gamma",nptbins,ptmin,ptmax);
2487       fhPtIsoFragmentation->SetYTitle("N");
2488       fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
2489       outputContainer->Add(fhPtIsoFragmentation) ;
2490
2491       fhPhiIsoFragmentation  = new TH2F
2492       ("hPhiMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2493       fhPhiIsoFragmentation->SetYTitle("#phi");
2494       fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
2495       outputContainer->Add(fhPhiIsoFragmentation) ;
2496
2497       fhEtaIsoFragmentation  = new TH2F
2498       ("hEtaMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2499       fhEtaIsoFragmentation->SetYTitle("#eta");
2500       fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
2501       outputContainer->Add(fhEtaIsoFragmentation) ;
2502
2503       fhPtIsoPi0  = new TH1F("hPtMCPi0","Number of isolated #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax);
2504       fhPtIsoPi0->SetYTitle("N");
2505       fhPtIsoPi0->SetXTitle("p_{T #gamma}(GeV/c)");
2506       outputContainer->Add(fhPtIsoPi0) ;
2507
2508       fhPhiIsoPi0  = new TH2F
2509       ("hPhiMCPi0","Number of isolated #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2510       fhPhiIsoPi0->SetYTitle("#phi");
2511       fhPhiIsoPi0->SetXTitle("p_{T #gamma} (GeV/c)");
2512       outputContainer->Add(fhPhiIsoPi0) ;
2513
2514       fhEtaIsoPi0  = new TH2F
2515       ("hEtaMCPi0","Number of isolated #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2516       fhEtaIsoPi0->SetYTitle("#eta");
2517       fhEtaIsoPi0->SetXTitle("p_{T #gamma} (GeV/c)");
2518       outputContainer->Add(fhEtaIsoPi0) ;
2519
2520       fhPtIsoPi0Decay  = new TH1F("hPtMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
2521       fhPtIsoPi0Decay->SetYTitle("N");
2522       fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
2523       outputContainer->Add(fhPtIsoPi0Decay) ;
2524
2525       fhPhiIsoPi0Decay  = new TH2F
2526       ("hPhiMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2527       fhPhiIsoPi0Decay->SetYTitle("#phi");
2528       fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
2529       outputContainer->Add(fhPhiIsoPi0Decay) ;
2530
2531       fhEtaIsoPi0Decay  = new TH2F
2532       ("hEtaMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2533       fhEtaIsoPi0Decay->SetYTitle("#eta");
2534       fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
2535       outputContainer->Add(fhEtaIsoPi0Decay) ;
2536
2537       fhPtIsoEtaDecay  = new TH1F("hPtMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax);
2538       fhPtIsoEtaDecay->SetYTitle("N");
2539       fhPtIsoEtaDecay->SetXTitle("p_{T #gamma}(GeV/c)");
2540       outputContainer->Add(fhPtIsoEtaDecay) ;
2541
2542       fhPhiIsoEtaDecay  = new TH2F
2543       ("hPhiMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2544       fhPhiIsoEtaDecay->SetYTitle("#phi");
2545       fhPhiIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
2546       outputContainer->Add(fhPhiIsoEtaDecay) ;
2547
2548       fhEtaIsoEtaDecay  = new TH2F
2549       ("hEtaMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2550       fhEtaIsoEtaDecay->SetYTitle("#eta");
2551       fhEtaIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
2552       outputContainer->Add(fhEtaIsoEtaDecay) ;
2553
2554       fhPtIsoOtherDecay  = new TH1F("hPtMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax);
2555       fhPtIsoOtherDecay->SetYTitle("N");
2556       fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
2557       outputContainer->Add(fhPtIsoOtherDecay) ;
2558
2559       fhPhiIsoOtherDecay  = new TH2F
2560       ("hPhiMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2561       fhPhiIsoOtherDecay->SetYTitle("#phi");
2562       fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
2563       outputContainer->Add(fhPhiIsoOtherDecay) ;
2564
2565       fhEtaIsoOtherDecay  = new TH2F
2566       ("hEtaMCOtherDecay","Number of isolated #gamma non-#pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2567       fhEtaIsoOtherDecay->SetYTitle("#eta");
2568       fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
2569       outputContainer->Add(fhEtaIsoOtherDecay) ;
2570
2571       //      fhPtIsoConversion  = new TH1F("hPtMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax);
2572       //      fhPtIsoConversion->SetYTitle("N");
2573       //      fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
2574       //      outputContainer->Add(fhPtIsoConversion) ;
2575       //
2576       //      fhPhiIsoConversion  = new TH2F
2577       //      ("hPhiMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2578       //      fhPhiIsoConversion->SetYTitle("#phi");
2579       //      fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
2580       //      outputContainer->Add(fhPhiIsoConversion) ;
2581       //
2582       //      fhEtaIsoConversion  = new TH2F
2583       //      ("hEtaMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2584       //      fhEtaIsoConversion->SetYTitle("#eta");
2585       //      fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
2586       //      outputContainer->Add(fhEtaIsoConversion) ;
2587
2588       fhPtIsoHadron  = new TH1F("hPtMCHadron","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax);
2589       fhPtIsoHadron->SetYTitle("N");
2590       fhPtIsoHadron->SetXTitle("p_{T}(GeV/c)");
2591       outputContainer->Add(fhPtIsoHadron) ;
2592
2593
2594       fhPhiIsoHadron  = new TH2F
2595       ("hPhiMCHadron","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2596       fhPhiIsoHadron->SetYTitle("#phi");
2597       fhPhiIsoHadron->SetXTitle("p_{T} (GeV/c)");
2598       outputContainer->Add(fhPhiIsoHadron) ;
2599
2600       fhEtaIsoHadron  = new TH2F
2601       ("hEtaMCHadron","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2602       fhEtaIsoHadron->SetYTitle("#eta");
2603       fhEtaIsoHadron->SetXTitle("p_{T} (GeV/c)");
2604       outputContainer->Add(fhEtaIsoHadron) ;
2605
2606     }//Histos with MC
2607
2608   }
2609
2610   // Not Isolated histograms, reference histograms
2611
2612   fhENoIso  = new TH1F("hENoIso",
2613                         Form("Number of not isolated leading particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
2614                         nptbins,ptmin,ptmax);
2615   fhENoIso->SetYTitle("N");
2616   fhENoIso->SetXTitle("E (GeV/c)");
2617   outputContainer->Add(fhENoIso) ;
2618
2619   fhPtNoIso  = new TH1F("hPtNoIso",
2620                         Form("Number of not isolated leading particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
2621                         nptbins,ptmin,ptmax);
2622   fhPtNoIso->SetYTitle("N");
2623   fhPtNoIso->SetXTitle("p_{T} (GeV/c)");
2624   outputContainer->Add(fhPtNoIso) ;
2625
2626   fhPtNLocMaxNoIso  = new TH2F("hPtNLocMaxNoIso",
2627                                Form("Number of not isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f vs NLM, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
2628                                nptbins,ptmin,ptmax,10,0,10);
2629   fhPtNLocMaxNoIso->SetYTitle("NLM");
2630   fhPtNLocMaxNoIso->SetXTitle("p_{T} (GeV/c)");
2631   outputContainer->Add(fhPtNLocMaxNoIso) ;
2632
2633   fhEtaPhiNoIso  = new TH2F("hEtaPhiNoIso",
2634                             Form("Number of not isolated leading particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
2635                             netabins,etamin,etamax,nphibins,phimin,phimax);
2636   fhEtaPhiNoIso->SetXTitle("#eta");
2637   fhEtaPhiNoIso->SetYTitle("#phi");
2638   outputContainer->Add(fhEtaPhiNoIso) ;
2639
2640   fhPtDecayNoIso  = new TH1F("hPtDecayNoIso",
2641                              Form("Number of not isolated leading pi0 decay particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
2642                              nptbins,ptmin,ptmax);
2643   fhPtDecayNoIso->SetYTitle("N");
2644   fhPtDecayNoIso->SetXTitle("p_{T} (GeV/c)");
2645   outputContainer->Add(fhPtDecayNoIso) ;
2646
2647   fhEtaPhiDecayNoIso  = new TH2F("hEtaPhiDecayNoIso",
2648                                  Form("Number of not isolated leading Pi0 decay particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
2649                                  netabins,etamin,etamax,nphibins,phimin,phimax);
2650   fhEtaPhiDecayNoIso->SetXTitle("#eta");
2651   fhEtaPhiDecayNoIso->SetYTitle("#phi");
2652   outputContainer->Add(fhEtaPhiDecayNoIso) ;
2653
2654
2655   if(IsDataMC())
2656   {
2657     fhPtNoIsoPi0  = new TH1F
2658     ("hPtNoIsoPi0","Number of not isolated leading #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax);
2659     fhPtNoIsoPi0->SetYTitle("N");
2660     fhPtNoIsoPi0->SetXTitle("p_{T} (GeV/c)");
2661     outputContainer->Add(fhPtNoIsoPi0) ;
2662
2663     fhPtNoIsoPi0Decay  = new TH1F
2664     ("hPtNoIsoPi0Decay","Number of not isolated leading #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
2665     fhPtNoIsoPi0Decay->SetYTitle("N");
2666     fhPtNoIsoPi0Decay->SetXTitle("p_{T} (GeV/c)");
2667     outputContainer->Add(fhPtNoIsoPi0Decay) ;
2668
2669     fhPtNoIsoEtaDecay  = new TH1F
2670     ("hPtNoIsoEtaDecay","Number of not isolated leading #gamma from eta decay",nptbins,ptmin,ptmax);
2671     fhPtNoIsoEtaDecay->SetYTitle("N");
2672     fhPtNoIsoEtaDecay->SetXTitle("p_{T} (GeV/c)");
2673     outputContainer->Add(fhPtNoIsoEtaDecay) ;
2674
2675     fhPtNoIsoOtherDecay  = new TH1F
2676     ("hPtNoIsoOtherDecay","Number of not isolated leading #gamma from other decay",nptbins,ptmin,ptmax);
2677     fhPtNoIsoOtherDecay->SetYTitle("N");
2678     fhPtNoIsoOtherDecay->SetXTitle("p_{T} (GeV/c)");
2679     outputContainer->Add(fhPtNoIsoOtherDecay) ;
2680
2681     fhPtNoIsoPrompt  = new TH1F
2682     ("hPtNoIsoPrompt","Number of not isolated leading prompt #gamma",nptbins,ptmin,ptmax);
2683     fhPtNoIsoPrompt->SetYTitle("N");
2684     fhPtNoIsoPrompt->SetXTitle("p_{T} (GeV/c)");
2685     outputContainer->Add(fhPtNoIsoPrompt) ;
2686
2687     fhPtIsoMCPhoton  = new TH1F
2688     ("hPtIsoMCPhoton","Number of isolated leading  #gamma",nptbins,ptmin,ptmax);
2689     fhPtIsoMCPhoton->SetYTitle("N");
2690     fhPtIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
2691     outputContainer->Add(fhPtIsoMCPhoton) ;
2692
2693     fhPtNoIsoMCPhoton  = new TH1F
2694     ("hPtNoIsoMCPhoton","Number of not isolated leading #gamma",nptbins,ptmin,ptmax);
2695     fhPtNoIsoMCPhoton->SetYTitle("N");
2696     fhPtNoIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
2697     outputContainer->Add(fhPtNoIsoMCPhoton) ;
2698
2699 //    fhPtNoIsoConversion  = new TH1F
2700 //    ("hPtNoIsoConversion","Number of not isolated leading conversion #gamma",nptbins,ptmin,ptmax);
2701 //    fhPtNoIsoConversion->SetYTitle("N");
2702 //    fhPtNoIsoConversion->SetXTitle("p_{T} (GeV/c)");
2703 //    outputContainer->Add(fhPtNoIsoConversion) ;
2704
2705     fhPtNoIsoFragmentation  = new TH1F
2706     ("hPtNoIsoFragmentation","Number of not isolated leading fragmentation #gamma",nptbins,ptmin,ptmax);
2707     fhPtNoIsoFragmentation->SetYTitle("N");
2708     fhPtNoIsoFragmentation->SetXTitle("p_{T} (GeV/c)");
2709     outputContainer->Add(fhPtNoIsoFragmentation) ;
2710
2711     fhPtNoIsoHadron  = new TH1F
2712     ("hPtNoIsoHadron","Number of not isolated leading hadrons",nptbins,ptmin,ptmax);
2713     fhPtNoIsoHadron->SetYTitle("N");
2714     fhPtNoIsoHadron->SetXTitle("p_{T} (GeV/c)");
2715     outputContainer->Add(fhPtNoIsoHadron) ;
2716
2717   }//Histos with MC
2718
2719
2720   if(fMakeSeveralIC)
2721   {
2722     const Int_t buffersize = 255;
2723     char name[buffersize];
2724     char title[buffersize];
2725     for(Int_t icone = 0; icone<fNCones; icone++)
2726     {
2727       // sum pt in cone vs. pt leading
2728       snprintf(name, buffersize,"hSumPtLeadingPt_Cone_%d",icone);
2729       snprintf(title, buffersize,"#Sigma p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
2730       fhSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2731       fhSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}p_{T} (GeV/c)");//#Sigma p_{T}
2732       fhSumPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
2733       outputContainer->Add(fhSumPtLeadingPt[icone]) ;
2734
2735       // pt in cone vs. pt leading
2736       snprintf(name, buffersize,"hPtLeadingPt_Cone_%d",icone);
2737       snprintf(title, buffersize,"p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
2738       fhPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2739       fhPtLeadingPt[icone] ->SetYTitle("p_{T}^{cone} (GeV/c)");
2740       fhPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
2741       outputContainer->Add(fhPtLeadingPt[icone]) ;
2742
2743        // sum pt in cone vs. pt leading in the forward region (for background subtraction studies)
2744         snprintf(name, buffersize,"hPerpSumPtLeadingPt_Cone_%d",icone);
2745       snprintf(title, buffersize,"#Sigma p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
2746       fhPerpSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2747       fhPerpSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}p_{T} (GeV/c)");//#Sigma p_{T}
2748       fhPerpSumPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
2749       outputContainer->Add(fhPerpSumPtLeadingPt[icone]) ;
2750
2751       // pt in cone vs. pt leading in the forward region (for background subtraction studies)
2752         snprintf(name, buffersize,"hPerpPtLeadingPt_Cone_%d",icone);
2753       snprintf(title, buffersize,"p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
2754       fhPerpPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2755       fhPerpPtLeadingPt[icone] ->SetYTitle("p_{T}^{cone} (GeV/c)");
2756       fhPerpPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
2757       outputContainer->Add(fhPerpPtLeadingPt[icone]) ;
2758
2759
2760       if(IsDataMC())
2761       {
2762         snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);
2763         snprintf(title, buffersize,"Candidate Prompt cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2764         fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2765         fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2766         fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
2767         outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
2768
2769         snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);
2770         snprintf(title, buffersize,"Candidate Fragmentation cone sum p_{T} for R = %2.2fvs candidate p_{T}",fConeSizes[icone]);
2771         fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2772         fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2773         fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
2774         outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
2775
2776         snprintf(name, buffersize,"hPtSumPi0_Cone_%d",icone);
2777         snprintf(title, buffersize,"Candidate Pi0 cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2778         fhPtSumIsolatedPi0[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2779         fhPtSumIsolatedPi0[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2780         fhPtSumIsolatedPi0[icone]->SetXTitle("p_{T} (GeV/c)");
2781         outputContainer->Add(fhPtSumIsolatedPi0[icone]) ;
2782
2783         snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);
2784         snprintf(title, buffersize,"Candidate Pi0Decay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2785         fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2786         fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2787         fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
2788         outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
2789
2790         snprintf(name, buffersize,"hPtSumEtaDecay_Cone_%d",icone);
2791         snprintf(title, buffersize,"Candidate EtaDecay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2792         fhPtSumIsolatedEtaDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2793         fhPtSumIsolatedEtaDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2794         fhPtSumIsolatedEtaDecay[icone]->SetXTitle("p_{T} (GeV/c)");
2795         outputContainer->Add(fhPtSumIsolatedEtaDecay[icone]) ;
2796
2797         snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);
2798         snprintf(title, buffersize,"Candidate OtherDecay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2799         fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2800         fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2801         fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
2802         outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
2803
2804 //        snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);
2805 //        snprintf(title, buffersize,"Candidate Conversion cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2806 //        fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2807 //        fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2808 //        fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
2809 //        outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
2810
2811         snprintf(name, buffersize,"hPtSumHadron_Cone_%d",icone);
2812         snprintf(title, buffersize,"Candidate Hadron cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
2813         fhPtSumIsolatedHadron[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2814         fhPtSumIsolatedHadron[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2815         fhPtSumIsolatedHadron[icone]->SetXTitle("p_{T} (GeV/c)");
2816         outputContainer->Add(fhPtSumIsolatedHadron[icone]) ;
2817
2818       }//Histos with MC
2819
2820       for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++)
2821       {
2822
2823         snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
2824         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
2825         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2826         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2827         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
2828
2829         snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
2830         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2831         fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2832         fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2833         outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
2834
2835
2836         snprintf(name, buffersize,"hPtSum_Cone_%d_Pt%d",icone,ipt);
2837         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2838         fhPtSumIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2839         // fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
2840         fhPtSumIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2841         outputContainer->Add(fhPtSumIsolated[icone][ipt]) ;
2842
2843         snprintf(name, buffersize,"hPtSumDensity_Cone_%d_Pt%d",icone,ipt);
2844         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for density in R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2845         fhPtSumDensityIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2846         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
2847         fhPtSumDensityIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2848         outputContainer->Add(fhPtSumDensityIso[icone][ipt]) ;
2849
2850         snprintf(name, buffersize,"hPtFracPtSum_Cone_%d_Pt%d",icone,ipt);
2851         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for PtFracPtSum in R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2852         fhPtFracPtSumIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2853         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
2854         fhPtFracPtSumIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2855         outputContainer->Add(fhPtFracPtSumIso[icone][ipt]) ;
2856
2857         // pt decays isolated
2858         snprintf(name, buffersize,"hPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
2859         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
2860         fhPtPtThresDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2861         fhPtPtThresDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2862         outputContainer->Add(fhPtPtThresDecayIso[icone][ipt]) ;
2863
2864         snprintf(name, buffersize,"hPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
2865         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2866         fhPtPtFracDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2867         fhPtPtFracDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2868         outputContainer->Add(fhPtPtFracDecayIso[icone][ipt]) ;
2869
2870         snprintf(name, buffersize,"hPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2871         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2872         fhPtPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2873         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2874         fhPtPtSumDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2875         outputContainer->Add(fhPtPtSumDecayIso[icone][ipt]) ;
2876
2877         snprintf(name, buffersize,"hPtSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
2878         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for density in R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2879         fhPtSumDensityDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2880         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2881         fhPtSumDensityDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2882         outputContainer->Add(fhPtSumDensityDecayIso[icone][ipt]) ;
2883
2884         snprintf(name, buffersize,"hPtFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2885         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for PtFracPtSum in R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2886         fhPtFracPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2887         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
2888         fhPtFracPtSumDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2889         outputContainer->Add(fhPtFracPtSumDecayIso[icone][ipt]) ;
2890
2891
2892         // eta:phi
2893         snprintf(name, buffersize,"hEtaPhiPtThres_Cone_%d_Pt%d",icone,ipt);
2894         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
2895         fhEtaPhiPtThresIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2896         fhEtaPhiPtThresIso[icone][ipt]->SetXTitle("#eta");
2897         fhEtaPhiPtThresIso[icone][ipt]->SetYTitle("#phi");
2898         outputContainer->Add(fhEtaPhiPtThresIso[icone][ipt]) ;
2899
2900         snprintf(name, buffersize,"hEtaPhiPtFrac_Cone_%d_Pt%d",icone,ipt);
2901         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2902         fhEtaPhiPtFracIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2903         fhEtaPhiPtFracIso[icone][ipt]->SetXTitle("#eta");
2904         fhEtaPhiPtFracIso[icone][ipt]->SetYTitle("#phi");
2905         outputContainer->Add(fhEtaPhiPtFracIso[icone][ipt]) ;
2906
2907         snprintf(name, buffersize,"hEtaPhiPtSum_Cone_%d_Pt%d",icone,ipt);
2908         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2909         fhEtaPhiPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2910         fhEtaPhiPtSumIso[icone][ipt]->SetXTitle("#eta");
2911         fhEtaPhiPtSumIso[icone][ipt]->SetYTitle("#phi");
2912         outputContainer->Add(fhEtaPhiPtSumIso[icone][ipt]) ;
2913
2914         snprintf(name, buffersize,"hEtaPhiSumDensity_Cone_%d_Pt%d",icone,ipt);
2915         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for density R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2916         fhEtaPhiSumDensityIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2917         fhEtaPhiSumDensityIso[icone][ipt]->SetXTitle("#eta");
2918         fhEtaPhiSumDensityIso[icone][ipt]->SetYTitle("#phi");
2919         outputContainer->Add(fhEtaPhiSumDensityIso[icone][ipt]) ;
2920
2921         snprintf(name, buffersize,"hEtaPhiFracPtSum_Cone_%d_Pt%d",icone,ipt);
2922         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for FracPtSum R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2923         fhEtaPhiFracPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2924         fhEtaPhiFracPtSumIso[icone][ipt]->SetXTitle("#eta");
2925         fhEtaPhiFracPtSumIso[icone][ipt]->SetYTitle("#phi");
2926         outputContainer->Add(fhEtaPhiFracPtSumIso[icone][ipt]) ;
2927
2928         // eta:phi decays
2929         snprintf(name, buffersize,"hEtaPhiPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
2930         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
2931         fhEtaPhiPtThresDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2932         fhEtaPhiPtThresDecayIso[icone][ipt]->SetXTitle("#eta");
2933         fhEtaPhiPtThresDecayIso[icone][ipt]->SetYTitle("#phi");
2934         outputContainer->Add(fhEtaPhiPtThresDecayIso[icone][ipt]) ;
2935
2936         snprintf(name, buffersize,"hEtaPhiPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
2937         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2938         fhEtaPhiPtFracDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2939         fhEtaPhiPtFracDecayIso[icone][ipt]->SetXTitle("#eta");
2940         fhEtaPhiPtFracDecayIso[icone][ipt]->SetYTitle("#phi");
2941         outputContainer->Add(fhEtaPhiPtFracDecayIso[icone][ipt]) ;
2942
2943
2944         snprintf(name, buffersize,"hEtaPhiPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2945         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2946         fhEtaPhiPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2947         fhEtaPhiPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
2948         fhEtaPhiPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
2949         outputContainer->Add(fhEtaPhiPtSumDecayIso[icone][ipt]) ;
2950
2951         snprintf(name, buffersize,"hEtaPhiSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
2952         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for density R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
2953         fhEtaPhiSumDensityDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2954         fhEtaPhiSumDensityDecayIso[icone][ipt]->SetXTitle("#eta");
2955         fhEtaPhiSumDensityDecayIso[icone][ipt]->SetYTitle("#phi");
2956         outputContainer->Add(fhEtaPhiSumDensityDecayIso[icone][ipt]) ;
2957
2958         snprintf(name, buffersize,"hEtaPhiFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2959         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for FracPtSum R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
2960         fhEtaPhiFracPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2961         fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
2962         fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
2963         outputContainer->Add(fhEtaPhiFracPtSumDecayIso[icone][ipt]) ;
2964
2965
2966         if(IsDataMC())
2967         {
2968           snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);
2969           snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
2970           fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2971           fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2972           outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
2973
2974           snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);
2975           snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
2976           fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2977           fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2978           outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
2979
2980           snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);
2981           snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
2982           fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2983           fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2984           outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
2985
2986           snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);
2987           snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
2988           fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2989           fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2990           outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
2991
2992           snprintf(name, buffersize,"hPtThresMCPi0_Cone_%d_Pt%d",icone,ipt);
2993           snprintf(title, buffersize,"Isolated candidate Pi0 p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
2994           fhPtThresIsolatedPi0[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2995           fhPtThresIsolatedPi0[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
2996           outputContainer->Add(fhPtThresIsolatedPi0[icone][ipt]) ;
2997
2998           snprintf(name, buffersize,"hPtFracMCPi0_Cone_%d_Pt%d",icone,ipt);
2999           snprintf(title, buffersize,"Isolated candidate Pi0 p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3000           fhPtFracIsolatedPi0[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3001           fhPtFracIsolatedPi0[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3002           outputContainer->Add(fhPtFracIsolatedPi0[icone][ipt]) ;
3003
3004           snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
3005           snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3006           fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3007           fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3008           outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
3009
3010           snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
3011           snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3012           fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3013           fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3014           outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
3015
3016           snprintf(name, buffersize,"hPtThresMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
3017           snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3018           fhPtThresIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3019           fhPtThresIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3020           outputContainer->Add(fhPtThresIsolatedEtaDecay[icone][ipt]) ;
3021
3022           snprintf(name, buffersize,"hPtFracMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
3023           snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3024           fhPtFracIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3025           fhPtFracIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3026           outputContainer->Add(fhPtFracIsolatedEtaDecay[icone][ipt]) ;
3027
3028
3029           snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
3030           snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3031           fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3032           fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3033           outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
3034
3035           snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
3036           snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3037           fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3038           fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3039           outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
3040
3041 //          snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);
3042 //          snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3043 //          fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3044 //          fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3045 //          outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
3046
3047 //          snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);
3048 //          snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3049 //          fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3050 //          fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3051 //          outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
3052
3053           snprintf(name, buffersize,"hPtThresMCHadron_Cone_%d_Pt%d",icone,ipt);
3054           snprintf(title, buffersize,"Isolated candidate Hadron p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3055           fhPtThresIsolatedHadron[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3056           fhPtThresIsolatedHadron[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3057           outputContainer->Add(fhPtThresIsolatedHadron[icone][ipt]) ;
3058
3059           snprintf(name, buffersize,"hPtFracMCHadron_Cone_%d_Pt%d",icone,ipt);
3060           snprintf(title, buffersize,"Isolated candidate Hadron p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
3061           fhPtFracIsolatedHadron[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
3062           fhPtFracIsolatedHadron[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
3063           outputContainer->Add(fhPtFracIsolatedHadron[icone][ipt]) ;
3064
3065         }//Histos with MC
3066       }//icone loop
3067     }//ipt loop
3068   }
3069
3070   if(fFillPileUpHistograms)
3071   {
3072     for (Int_t i = 0; i < 7 ; i++)
3073     {
3074       fhEIsoPileUp[i]   = new TH1F(Form("hEPileUp%s",pileUpName[i].Data()),
3075                                 Form("Number of isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f, pile-up event by %s",r,ptthre,ptfrac,pileUpName[i].Data()),
3076                                 nptbins,ptmin,ptmax);
3077       fhEIsoPileUp[i]->SetYTitle("dN / dE");
3078       fhEIsoPileUp[i]->SetXTitle("E (GeV)");
3079       outputContainer->Add(fhEIsoPileUp[i]) ;
3080
3081       fhPtIsoPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
3082                                 Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr}= %2.2f, pile-up event by %s ",r,ptthre,ptfrac,pileUpName[i].Data()),
3083                                 nptbins,ptmin,ptmax);
3084       fhPtIsoPileUp[i]->SetYTitle("dN / p_{T}");
3085       fhPtIsoPileUp[i]->SetXTitle("p_{T} (GeV/c)");
3086       outputContainer->Add(fhPtIsoPileUp[i]) ;
3087
3088       fhENoIsoPileUp[i]   = new TH1F(Form("hENoIsoPileUp%s",pileUpName[i].Data()),
3089                                   Form("Number of not isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f, pile-up event by %s",r,ptthre,ptfrac,pileUpName[i].Data()),
3090                                   nptbins,ptmin,ptmax);
3091       fhENoIsoPileUp[i]->SetYTitle("dN / dE");
3092       fhENoIsoPileUp[i]->SetXTitle("E (GeV)");
3093       outputContainer->Add(fhENoIsoPileUp[i]) ;
3094
3095       fhPtNoIsoPileUp[i]  = new TH1F(Form("hPtNoIsoPileUp%s",pileUpName[i].Data()),
3096                                   Form("Number of not isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr}= %2.2f, pile-up event by %s ",r,ptthre,ptfrac,pileUpName[i].Data()),
3097                                   nptbins,ptmin,ptmax);
3098       fhPtNoIsoPileUp[i]->SetYTitle("dN / p_{T}");
3099       fhPtNoIsoPileUp[i]->SetXTitle("p_{T} (GeV/c)");
3100       outputContainer->Add(fhPtNoIsoPileUp[i]) ;
3101     }
3102
3103     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3104     fhTimeENoCut->SetXTitle("E (GeV)");
3105     fhTimeENoCut->SetYTitle("time (ns)");
3106     outputContainer->Add(fhTimeENoCut);
3107
3108     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3109     fhTimeESPD->SetXTitle("E (GeV)");
3110     fhTimeESPD->SetYTitle("time (ns)");
3111     outputContainer->Add(fhTimeESPD);
3112
3113     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3114     fhTimeESPDMulti->SetXTitle("E (GeV)");
3115     fhTimeESPDMulti->SetYTitle("time (ns)");
3116     outputContainer->Add(fhTimeESPDMulti);
3117
3118     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3119     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
3120     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
3121     outputContainer->Add(fhTimeNPileUpVertSPD);
3122
3123     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
3124     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
3125     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
3126     outputContainer->Add(fhTimeNPileUpVertTrack);
3127
3128     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3129     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
3130     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
3131     outputContainer->Add(fhTimeNPileUpVertContributors);
3132
3133     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);
3134     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
3135     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
3136     outputContainer->Add(fhTimePileUpMainVertexZDistance);
3137
3138     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
3139     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
3140     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
3141     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
3142   }
3143
3144   return outputContainer ;
3145
3146 }
3147
3148 //__________________________________
3149 void AliAnaParticleIsolation::Init()
3150 {
3151   // Do some checks and init stuff
3152
3153   // In case of several cone and thresholds analysis, open the cuts for the filling of the
3154   // track and cluster reference arrays in cone when done in the MakeAnalysisFillAOD().
3155   // The different cones, thresholds are tested for this list of tracks, clusters.
3156   if(fMakeSeveralIC)
3157   {
3158     printf("AliAnaParticleIsolation::Init() - Open default isolation cuts for multiple Isolation analysis\n");
3159     GetIsolationCut()->SetPtThreshold(100);
3160     GetIsolationCut()->SetPtFraction(100);
3161     GetIsolationCut()->SetConeSize(1);
3162   }
3163 }
3164
3165 //____________________________________________
3166 void AliAnaParticleIsolation::InitParameters()
3167 {
3168
3169   //Initialize the parameters of the analysis.
3170   SetInputAODName("PWG4Particle");
3171   SetAODObjArrayName("IsolationCone");
3172   AddToHistogramsName("AnaIsolation_");
3173
3174   fCalorimeter = "PHOS" ;
3175   fReMakeIC = kFALSE ;
3176   fMakeSeveralIC = kFALSE ;
3177
3178   //----------- Several IC-----------------
3179   fNCones             = 5 ;
3180   fNPtThresFrac       = 5 ;
3181   fConeSizes      [0] = 0.1;     fConeSizes      [1] = 0.2;   fConeSizes      [2] = 0.3; fConeSizes      [3] = 0.4;  fConeSizes      [4] = 0.5;
3182   fPtThresholds   [0] = 1.;      fPtThresholds   [1] = 2.;    fPtThresholds   [2] = 3.;  fPtThresholds   [3] = 4.;   fPtThresholds   [4] = 5.;
3183   fPtFractions    [0] = 0.05;    fPtFractions    [1] = 0.075; fPtFractions    [2] = 0.1; fPtFractions    [3] = 1.25; fPtFractions    [4] = 1.5;
3184   fSumPtThresholds[0] = 1.;      fSumPtThresholds[1] = 2.;    fSumPtThresholds[2] = 3.;  fSumPtThresholds[3] = 4.;   fSumPtThresholds[4] = 5.;
3185
3186   //------------- Histograms settings -------
3187   fHistoNPtSumBins = 100 ;
3188   fHistoPtSumMax   = 50 ;
3189   fHistoPtSumMin   = 0.  ;
3190
3191   fHistoNPtInConeBins = 100 ;
3192   fHistoPtInConeMax   = 50 ;
3193   fHistoPtInConeMin   = 0.  ;
3194
3195 }
3196
3197 //__________________________________________________
3198 void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
3199 {
3200   //Do analysis and fill aods
3201   //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
3202
3203   if(!GetInputAODBranch())
3204   {
3205     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
3206     abort();
3207   }
3208
3209   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
3210   {
3211     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
3212     abort();
3213   }
3214
3215   Int_t n = 0, nfrac = 0;
3216   Bool_t  isolated  = kFALSE ;
3217   Float_t coneptsum = 0 ;
3218   TObjArray * pl    = 0x0; ;
3219
3220   //Select the calorimeter for candidate isolation with neutral particles
3221   if      (fCalorimeter == "PHOS" )
3222     pl = GetPHOSClusters();
3223   else if (fCalorimeter == "EMCAL")
3224     pl = GetEMCALClusters();
3225
3226   //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
3227   Double_t ptLeading = 0. ;
3228   Int_t    idLeading = -1 ;
3229   TLorentzVector mom ;
3230   Int_t naod = GetInputAODBranch()->GetEntriesFast();
3231   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
3232
3233   for(Int_t iaod = 0; iaod < naod; iaod++)
3234   {
3235     AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3236
3237     //If too small or too large pt, skip
3238     if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ;
3239
3240     //check if it is low pt trigger particle
3241     if((aodinput->Pt() < GetIsolationCut()->GetPtThreshold() ||
3242         aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold()) &&
3243        !fMakeSeveralIC)
3244     {
3245       continue ; //trigger should not come from underlying event
3246     }
3247
3248     //vertex cut in case of mixing
3249     Int_t check = CheckMixedEventVertex(aodinput->GetCaloLabel(0), aodinput->GetTrackLabel(0));
3250     if(check ==  0) continue;
3251     if(check == -1) return;
3252
3253     //find the leading particles with highest momentum
3254     if ( aodinput->Pt() > ptLeading )
3255     {
3256       ptLeading = aodinput->Pt() ;
3257       idLeading = iaod ;
3258     }
3259
3260     aodinput->SetLeadingParticle(kFALSE);
3261
3262   }//finish searching for leading trigger particle
3263
3264   // Check isolation of leading particle
3265   if(idLeading < 0) return;
3266
3267   AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(idLeading));
3268   aodinput->SetLeadingParticle(kTRUE);
3269
3270   // Check isolation only of clusters in fiducial region
3271   if(IsFiducialCutOn())
3272   {
3273     Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodinput->Momentum(),aodinput->GetDetector()) ;
3274     if(! in ) return ;
3275   }
3276
3277   //After cuts, study isolation
3278   n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
3279   GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,
3280                                       GetReader(), GetCaloPID(),
3281                                       kTRUE, aodinput, GetAODObjArrayName(),
3282                                       n,nfrac,coneptsum, isolated);
3283
3284   if(!fMakeSeveralIC) aodinput->SetIsolated(isolated);
3285
3286   if(GetDebug() > 1)
3287   {
3288     if(isolated)printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
3289     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
3290   }
3291
3292 }
3293
3294 //_________________________________________________________
3295 void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
3296 {
3297   //Do analysis and fill histograms
3298
3299   //Loop on stored AOD
3300   Int_t naod = GetInputAODBranch()->GetEntriesFast();
3301   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
3302
3303   for(Int_t iaod = 0; iaod < naod ; iaod++)
3304   {
3305     AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3306
3307     if(!aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
3308
3309     // Check isolation only of clusters in fiducial region
3310     if(IsFiducialCutOn())
3311     {
3312       Bool_t in = GetFiducialCut()->IsInFiducialCut(*aod->Momentum(),aod->GetDetector()) ;
3313       if(! in ) continue ;
3314     }
3315
3316     Float_t pt         = aod->Pt();
3317     //If too small or too large pt, skip
3318     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3319
3320     Bool_t  isolated   = aod->IsIsolated();
3321     Bool_t  decay      = aod->IsTagged();
3322     Float_t energy     = aod->E();
3323     Float_t phi        = aod->Phi();
3324     Float_t eta        = aod->Eta();
3325
3326     // --- In case of redoing isolation from delta AOD ----
3327
3328     if     (fMakeSeveralIC)
3329     {
3330       //Analysis of multiple IC at same time
3331       MakeSeveralICAnalysis(aod);
3332       continue;
3333     }
3334     else if(fReMakeIC)
3335     {
3336       //In case a more strict IC is needed in the produced AOD
3337       isolated = kFALSE;
3338       Int_t   n = 0, nfrac = 0;
3339       Float_t coneptsum = 0 ;
3340
3341       //Recover reference arrays with clusters and tracks
3342       TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
3343       TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
3344
3345       GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
3346                                           GetReader(), GetCaloPID(),
3347                                           kFALSE, aod, "",
3348                                           n,nfrac,coneptsum, isolated);
3349
3350       fhConeSumPt          ->Fill(pt,     coneptsum);
3351       fhConeSumPtTrigEtaPhi->Fill(eta,phi,coneptsum);
3352
3353       if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
3354     }
3355     else
3356     {
3357       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f\n",pt, eta, phi);
3358
3359       FillTrackMatchingShowerShapeControlHistograms(aod,GetCaloPID());
3360
3361       //Fill pt/sum pT distribution of particles in cone or in UE band
3362       Float_t coneptsumCluster = 0;
3363       Float_t coneptsumTrack   = 0;
3364       Float_t coneptsumCell    = 0;
3365       Float_t etaBandptsumClusterNorm = 0;
3366       Float_t etaBandptsumTrackNorm   = 0;
3367
3368       CalculateTrackSignalInCone   (aod,coneptsumTrack  );
3369       CalculateCaloSignalInCone    (aod,coneptsumCluster);
3370       CalculateCaloCellSignalInCone(aod,coneptsumCell   );
3371
3372       if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged)
3373       {
3374         fhConeSumPtClustervsTrack     ->Fill(coneptsumCluster,coneptsumTrack);
3375         fhConeSumPtCellvsTrack        ->Fill(coneptsumCell,   coneptsumTrack);
3376         fhConeSumPtCellTrack          ->Fill(pt,     coneptsumTrack+coneptsumCell);
3377         fhConeSumPtCellTrackTrigEtaPhi->Fill(eta,phi,coneptsumTrack+coneptsumCell);
3378       }
3379
3380       fhConeSumPt              ->Fill(pt,     coneptsumTrack+coneptsumCluster);
3381       fhConeSumPtTrigEtaPhi    ->Fill(eta,phi,coneptsumTrack+coneptsumCluster);
3382
3383       if(GetDebug() > 1)
3384         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsumTrack+coneptsumCluster);
3385
3386       //normalize phi/eta band per area unit
3387       CalculateNormalizeUEBandPerUnitArea(aod, coneptsumCluster, coneptsumCell, coneptsumTrack, etaBandptsumTrackNorm, etaBandptsumClusterNorm) ;
3388
3389     //  printf("Histograms analysis : cluster pt = %f, etaBandTrack = %f, etaBandCluster = %f, isolation = %d\n",aod->Pt(),etaBandptsumTrackNorm,etaBandptsumClusterNorm,aod->IsIsolated());
3390
3391     }
3392
3393     Int_t mcTag = aod->GetTag() ;
3394
3395     if(isolated)
3396     {
3397       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED: fill histograms\n", iaod);
3398
3399       // Fill histograms to undertand pile-up before other cuts applied
3400       // Remember to relax time cuts in the reader
3401       FillPileUpHistograms(aod->GetCaloLabel(0));
3402
3403       fhEIso      ->Fill(energy);
3404       fhPtIso     ->Fill(pt);
3405       fhPhiIso    ->Fill(pt,phi);
3406       fhEtaIso    ->Fill(pt,eta);
3407       fhEtaPhiIso ->Fill(eta,phi);
3408
3409       fhPtNLocMaxIso    ->Fill(pt,aod->GetFiducialArea()) ;
3410       fhPtCentralityIso ->Fill(pt,GetEventCentrality()) ;
3411       fhPtEventPlaneIso ->Fill(pt,GetEventPlaneAngle() ) ;
3412
3413       if(decay)
3414       {
3415         fhPtDecayIso    ->Fill(pt);
3416         fhEtaPhiDecayIso->Fill(eta,phi);
3417       }
3418
3419       if(fFillPileUpHistograms)
3420       {
3421         if(GetReader()->IsPileUpFromSPD())               { fhEIsoPileUp[0] ->Fill(energy) ; fhPtIsoPileUp[0]->Fill(pt) ; }
3422         if(GetReader()->IsPileUpFromEMCal())             { fhEIsoPileUp[1] ->Fill(energy) ; fhPtIsoPileUp[1]->Fill(pt) ; }
3423         if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhEIsoPileUp[2] ->Fill(energy) ; fhPtIsoPileUp[2]->Fill(pt) ; }
3424         if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhEIsoPileUp[3] ->Fill(energy) ; fhPtIsoPileUp[3]->Fill(pt) ; }
3425         if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhEIsoPileUp[4] ->Fill(energy) ; fhPtIsoPileUp[4]->Fill(pt) ; }
3426         if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhEIsoPileUp[5] ->Fill(energy) ; fhPtIsoPileUp[5]->Fill(pt) ; }
3427         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhEIsoPileUp[6] ->Fill(energy) ; fhPtIsoPileUp[6]->Fill(pt) ; }
3428       }
3429
3430       if(IsDataMC())
3431       {
3432
3433         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3434         {
3435           fhPtIsoMCPhoton  ->Fill(pt);
3436         }
3437
3438         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))
3439         {
3440           fhPtIsoPrompt  ->Fill(pt);
3441           fhPhiIsoPrompt ->Fill(pt,phi);
3442           fhEtaIsoPrompt ->Fill(pt,eta);
3443         }
3444         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation))
3445         {
3446           fhPtIsoFragmentation  ->Fill(pt);
3447           fhPhiIsoFragmentation ->Fill(pt,phi);
3448           fhEtaIsoFragmentation ->Fill(pt,eta);
3449         }
3450         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))
3451         {
3452           fhPtIsoPi0  ->Fill(pt);
3453           fhPhiIsoPi0 ->Fill(pt,phi);
3454           fhEtaIsoPi0 ->Fill(pt,eta);
3455         }
3456         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))
3457         {
3458           fhPtIsoPi0Decay  ->Fill(pt);
3459           fhPhiIsoPi0Decay ->Fill(pt,phi);
3460           fhEtaIsoPi0Decay ->Fill(pt,eta);
3461         }
3462         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))
3463         {
3464           fhPtIsoEtaDecay  ->Fill(pt);
3465           fhPhiIsoEtaDecay ->Fill(pt,phi);
3466           fhEtaIsoEtaDecay ->Fill(pt,eta);
3467         }
3468         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))
3469         {
3470           fhPtIsoOtherDecay  ->Fill(pt);
3471           fhPhiIsoOtherDecay ->Fill(pt,phi);
3472           fhEtaIsoOtherDecay ->Fill(pt,eta);
3473         }
3474         //        else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))
3475         //        {
3476         //          fhPtIsoConversion  ->Fill(pt);
3477         //          fhPhiIsoConversion ->Fill(pt,phi);
3478         //          fhEtaIsoConversion ->Fill(pt,eta);
3479         //        }
3480         else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))// anything else but electrons
3481         {
3482           fhPtIsoHadron  ->Fill(pt);
3483           fhPhiIsoHadron ->Fill(pt,phi);
3484           fhEtaIsoHadron ->Fill(pt,eta);
3485         }
3486       }//Histograms with MC
3487
3488     }//Isolated histograms
3489     else // NON isolated
3490     {
3491       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d NOT ISOLATED, fill histograms\n", iaod);
3492
3493       fhENoIso        ->Fill(energy);
3494       fhPtNoIso       ->Fill(pt);
3495       fhEtaPhiNoIso   ->Fill(eta,phi);
3496       fhPtNLocMaxNoIso->Fill(pt,aod->GetFiducialArea());
3497
3498       if(fFillPileUpHistograms)
3499       {
3500         if(GetReader()->IsPileUpFromSPD())                { fhENoIsoPileUp[0] ->Fill(energy) ; fhPtNoIsoPileUp[0]->Fill(pt) ; }
3501         if(GetReader()->IsPileUpFromEMCal())              { fhENoIsoPileUp[1] ->Fill(energy) ; fhPtNoIsoPileUp[1]->Fill(pt) ; }
3502         if(GetReader()->IsPileUpFromSPDOrEMCal())         { fhENoIsoPileUp[2] ->Fill(energy) ; fhPtNoIsoPileUp[2]->Fill(pt) ; }
3503         if(GetReader()->IsPileUpFromSPDAndEMCal())        { fhENoIsoPileUp[3] ->Fill(energy) ; fhPtNoIsoPileUp[3]->Fill(pt) ; }
3504         if(GetReader()->IsPileUpFromSPDAndNotEMCal())     { fhENoIsoPileUp[4] ->Fill(energy) ; fhPtNoIsoPileUp[4]->Fill(pt) ; }
3505         if(GetReader()->IsPileUpFromEMCalAndNotSPD())     { fhENoIsoPileUp[5] ->Fill(energy) ; fhPtNoIsoPileUp[5]->Fill(pt) ; }
3506         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhENoIsoPileUp[6] ->Fill(energy) ; fhPtNoIsoPileUp[6]->Fill(pt) ; }
3507       }
3508
3509       if(decay)
3510       {
3511         fhPtDecayNoIso    ->Fill(pt);
3512         fhEtaPhiDecayNoIso->Fill(eta,phi);
3513       }
3514
3515       if(IsDataMC())
3516       {
3517         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))        fhPtNoIsoMCPhoton     ->Fill(pt);
3518         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))           fhPtNoIsoPi0          ->Fill(pt);
3519         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtNoIsoPi0Decay     ->Fill(pt);
3520         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtNoIsoEtaDecay     ->Fill(pt);
3521         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtNoIsoOtherDecay   ->Fill(pt);
3522         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))        fhPtNoIsoPrompt       ->Fill(pt);
3523         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) fhPtNoIsoFragmentation->Fill(pt);
3524         //        else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(pt);
3525         else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))     fhPtNoIsoHadron      ->Fill(pt);
3526       }
3527     }
3528   }// aod loop
3529
3530 }
3531
3532 //_____________________________________________________________________________________
3533 void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph)
3534 {
3535
3536   //Isolation Cut Analysis for both methods and different pt cuts and cones
3537   Float_t ptC   = ph->Pt();
3538   Float_t etaC  = ph->Eta();
3539   Float_t phiC  = ph->Phi();
3540   Int_t   tag   = ph->GetTag();
3541   Bool_t  decay = ph->IsTagged();
3542
3543   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f\n",ptC);
3544
3545   //Keep original setting used when filling AODs, reset at end of analysis
3546   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
3547   Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();
3548   Float_t rorg       = GetIsolationCut()->GetConeSize();
3549
3550   Float_t coneptsum = 0 ;
3551   Int_t   n    [10][10];//[fNCones][fNPtThresFrac];
3552   Int_t   nfrac[10][10];//[fNCones][fNPtThresFrac];
3553   Bool_t  isolated  = kFALSE;
3554   Int_t   nCone     = 0;
3555   Int_t   nFracCone = 0;
3556
3557   // Fill hist with all particles before isolation criteria
3558   fhPtNoIso    ->Fill(ptC);
3559   fhEtaPhiNoIso->Fill(etaC,phiC);
3560
3561   if(IsDataMC())
3562   {
3563     if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))        fhPtNoIsoMCPhoton     ->Fill(ptC);
3564     if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtNoIsoPi0          ->Fill(ptC);
3565     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtNoIsoPi0Decay     ->Fill(ptC);
3566     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtNoIsoEtaDecay     ->Fill(ptC);
3567     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtNoIsoOtherDecay   ->Fill(ptC);
3568     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtNoIsoPrompt       ->Fill(ptC);
3569     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtNoIsoFragmentation->Fill(ptC);
3570 //    else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(ptC);
3571     else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))    fhPtNoIsoHadron      ->Fill(ptC);
3572   }
3573
3574   if(decay)
3575   {
3576     fhPtDecayNoIso    ->Fill(ptC);
3577     fhEtaPhiDecayNoIso->Fill(etaC,phiC);
3578   }
3579   //Get vertex for photon momentum calculation
3580   Double_t vertex[] = {0,0,0} ; //vertex ;
3581   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3582       GetReader()->GetVertex(vertex);
3583
3584   //Loop on cone sizes
3585   for(Int_t icone = 0; icone<fNCones; icone++)
3586   {
3587     //Recover reference arrays with clusters and tracks
3588     TObjArray * refclusters = ph->GetObjArray(GetAODObjArrayName()+"Clusters");
3589     TObjArray * reftracks   = ph->GetObjArray(GetAODObjArrayName()+"Tracks");
3590
3591     //If too small or too large pt, skip
3592     if(ptC < GetMinPt() || ptC > GetMaxPt() ) continue ;
3593
3594    //In case a more strict IC is needed in the produced AOD
3595
3596     nCone=0; nFracCone = 0; isolated = kFALSE; coneptsum = 0;
3597
3598   GetIsolationCut()->SetSumPtThreshold(100);
3599   GetIsolationCut()->SetPtThreshold(100);
3600   GetIsolationCut()->SetPtFraction(100);
3601   GetIsolationCut()->SetConeSize(fConeSizes[icone]);
3602   GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
3603                                           GetReader(), GetCaloPID(),
3604                                           kFALSE, ph, "",
3605                                           nCone,nFracCone,coneptsum, isolated);
3606
3607     fhSumPtLeadingPt[icone]->Fill(ptC,coneptsum);
3608
3609     // retreive pt tracks to fill histo vs. pt leading
3610     //Fill pt distribution of particles in cone
3611     //fhPtLeadingPt(),fhPerpSumPtLeadingPt(),fhPerpPtLeadingPt(),
3612
3613     //Tracks
3614     coneptsum = 0;
3615     Double_t sumptPerp = 0. ;
3616     TObjArray * trackList   = GetCTSTracks() ;
3617     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
3618     {
3619       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
3620       //fill the histograms at forward range
3621       if(!track)
3622       {
3623         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
3624         continue;
3625       }
3626
3627       Double_t dPhi = phiC - track->Phi() + TMath::PiOver2();
3628       Double_t dEta = etaC - track->Eta();
3629       Double_t arg  = dPhi*dPhi + dEta*dEta;
3630       if(TMath::Sqrt(arg) < fConeSizes[icone])
3631       {
3632         fhPerpPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
3633         sumptPerp+=track->Pt();
3634       }
3635
3636       dPhi = phiC - track->Phi() - TMath::PiOver2();
3637       arg  = dPhi*dPhi + dEta*dEta;
3638       if(TMath::Sqrt(arg) < fConeSizes[icone])
3639       {
3640         fhPerpPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
3641         sumptPerp+=track->Pt();
3642       }
3643     }
3644
3645     fhPerpSumPtLeadingPt[icone]->Fill(ptC,sumptPerp);
3646
3647     if(reftracks)
3648     {
3649       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
3650       {
3651         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
3652         fhPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
3653         coneptsum+=track->Pt();
3654       }
3655     }
3656
3657     //CaloClusters
3658     if(refclusters)
3659     {
3660       TLorentzVector mom ;
3661       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
3662       {
3663         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
3664         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
3665
3666         fhPtLeadingPt[icone]->Fill(ptC, mom.Pt());
3667         coneptsum+=mom.Pt();
3668       }
3669     }
3670     ///////////////////
3671
3672
3673     //Loop on ptthresholds
3674     for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++)
3675     {
3676       n    [icone][ipt]=0;
3677       nfrac[icone][ipt]=0;
3678       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
3679       GetIsolationCut()->SetPtFraction(fPtFractions[ipt]) ;
3680       GetIsolationCut()->SetSumPtThreshold(fSumPtThresholds[ipt]);
3681
3682       GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
3683                                           GetReader(), GetCaloPID(),
3684                                           kFALSE, ph, "",
3685                                           n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
3686
3687       if(!isolated) continue;
3688       //Normal ptThreshold cut
3689
3690       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - cone size %1.1f, ptThres  %1.1f, sumptThresh  %1.1f, n %d, nfrac %d, coneptsum %2.2f, isolated %d\n",
3691                                 fConeSizes[icone],fPtThresholds[ipt],fSumPtThresholds[ipt],n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
3692       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - pt %1.1f, eta %1.1f, phi %1.1f\n",ptC, etaC, phiC);
3693
3694       if(n[icone][ipt] == 0)
3695       {
3696         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling pt threshold loop\n");
3697         fhPtThresIsolated[icone][ipt]->Fill(ptC);
3698         fhEtaPhiPtThresIso[icone][ipt]->Fill(etaC,phiC);
3699
3700         if(decay)
3701         {
3702           fhPtPtThresDecayIso[icone][ipt]->Fill(ptC);
3703           //      fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
3704         }
3705
3706         if(IsDataMC())
3707         {
3708           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtThresIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
3709 //          else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtThresIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
3710           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
3711           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtThresIsolatedPi0[icone][ipt]           ->Fill(ptC) ;
3712           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtThresIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
3713           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtThresIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
3714           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtThresIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
3715           else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))      fhPtThresIsolatedHadron[icone][ipt]      ->Fill(ptC) ;
3716         }
3717       }
3718
3719       // pt in cone fraction
3720       if(nfrac[icone][ipt] == 0)
3721       {
3722         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling frac loop\n");
3723         fhPtFracIsolated[icone][ipt]->Fill(ptC);
3724         fhEtaPhiPtFracIso[icone][ipt]->Fill(etaC,phiC);
3725
3726         if(decay)
3727         {
3728           fhPtPtFracDecayIso[icone][ipt]->Fill(ptC);
3729           fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
3730         }
3731
3732         if(IsDataMC())
3733         {
3734           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtFracIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
3735 //          else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtFracIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
3736           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
3737           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtFracIsolatedPi0[icone][ipt]          ->Fill(ptC) ;
3738           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtFracIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
3739           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtFracIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
3740           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtFracIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
3741           else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))      fhPtFracIsolatedHadron[icone][ipt]->Fill(ptC) ;
3742         }
3743       }
3744
3745       if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - checking IC method : %i\n",GetIsolationCut()->GetICMethod());
3746
3747       //Pt threshold on pt cand/ sum in cone histograms
3748       if(coneptsum<fSumPtThresholds[ipt])
3749       {//      if((GetIsolationCut()->GetICMethod())==1){//kSumPtIC){
3750         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling sum loop\n");
3751         fhPtSumIsolated[icone][ipt]->Fill(ptC) ;
3752         fhEtaPhiPtSumIso[icone][ipt]->Fill(etaC, phiC) ;
3753         if(decay)
3754         {
3755           fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
3756           fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
3757         }
3758       }
3759
3760     // pt sum pt frac method
3761 //    if( ((fPtFractions[ipt]*ptC < fSumPtThresholds[ipt]) && (coneptsum < fSumPtThresholds[ipt])) || ((fPtFractions[ipt]*ptC > fSumPtThresholds[ipt]) && (coneptsum < fPtFractions[ipt]*ptC)) )
3762
3763       if(coneptsum < fPtFractions[ipt]*ptC)
3764        {
3765         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling PtFrac PtSum loop\n");
3766         fhPtFracPtSumIso[icone][ipt]->Fill(ptC) ;
3767         fhEtaPhiFracPtSumIso[icone][ipt]->Fill(etaC,phiC) ;
3768
3769         if(decay)
3770         {
3771           fhPtFracPtSumDecayIso[icone][ipt]->Fill(ptC);
3772           fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
3773         }
3774       }
3775
3776       // density method
3777       Float_t cellDensity = GetIsolationCut()->GetCellDensity( ph, GetReader());
3778       if(coneptsum<fSumPtThresholds[ipt]*cellDensity)
3779       {//(GetIsolationCut()->GetICMethod())==4){//kSumDensityIC) {
3780         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling density loop\n");
3781         fhPtSumDensityIso[icone][ipt]->Fill(ptC) ;
3782         fhEtaPhiSumDensityIso[icone][ipt]->Fill(etaC,phiC) ;
3783
3784         if(decay)
3785         {
3786           fhPtSumDensityDecayIso[icone][ipt]->Fill(ptC);
3787           fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
3788         }
3789
3790       }
3791     }//pt thresh loop
3792
3793     if(IsDataMC())
3794     {
3795       if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtSumIsolatedPrompt[icone]       ->Fill(ptC,coneptsum) ;
3796 //      else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtSumIsolatedConversion[icone]   ->Fill(ptC,coneptsum) ;
3797       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
3798       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtSumIsolatedPi0[icone]          ->Fill(ptC,coneptsum) ;
3799       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtSumIsolatedPi0Decay[icone]     ->Fill(ptC,coneptsum) ;
3800       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtSumIsolatedEtaDecay[icone]     ->Fill(ptC,coneptsum) ;
3801       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtSumIsolatedOtherDecay[icone]   ->Fill(ptC,coneptsum) ;
3802       else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))      fhPtSumIsolatedHadron[icone]->Fill(ptC,coneptsum) ;
3803     }
3804
3805   }//cone size loop
3806
3807   //Reset original parameters for AOD analysis
3808   GetIsolationCut()->SetPtThreshold(ptthresorg);
3809   GetIsolationCut()->SetPtFraction(ptfracorg);
3810   GetIsolationCut()->SetConeSize(rorg);
3811
3812 }
3813
3814 //_____________________________________________________________
3815 void AliAnaParticleIsolation::Print(const Option_t * opt) const
3816 {
3817
3818   //Print some relevant parameters set for the analysis
3819   if(! opt)
3820     return;
3821
3822   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3823   AliAnaCaloTrackCorrBaseClass::Print(" ");
3824
3825   printf("ReMake Isolation          = %d \n",  fReMakeIC) ;
3826   printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;
3827   printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;
3828
3829   if(fMakeSeveralIC)
3830   {
3831     printf("N Cone Sizes       =     %d\n", fNCones) ;
3832     printf("Cone Sizes          =    \n") ;
3833     for(Int_t i = 0; i < fNCones; i++)
3834       printf("  %1.2f;",  fConeSizes[i]) ;
3835     printf("    \n") ;
3836
3837     printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
3838     printf(" pT thresholds         =    \n") ;
3839     for(Int_t i = 0; i < fNPtThresFrac; i++)
3840       printf("   %2.2f;",  fPtThresholds[i]) ;
3841
3842     printf("    \n") ;
3843
3844     printf(" pT fractions         =    \n") ;
3845     for(Int_t i = 0; i < fNPtThresFrac; i++)
3846       printf("   %2.2f;",  fPtFractions[i]) ;
3847
3848     printf("    \n") ;
3849
3850     printf("sum pT thresholds         =    \n") ;
3851     for(Int_t i = 0; i < fNPtThresFrac; i++)
3852       printf("   %2.2f;",  fSumPtThresholds[i]) ;
3853
3854
3855   }
3856
3857   printf("Histograms: %3.1f < pT sum < %3.1f,  Nbin = %d\n",    fHistoPtSumMin,    fHistoPtSumMax,    fHistoNPtSumBins   );
3858   printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);
3859
3860   printf("    \n") ;
3861
3862 }
3863