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