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