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