PWGGA/CaloTrackCorrelations/AliAnaParticleIsolation.cxx

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