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