PWGGA/CaloTrackCorrelations/AliAnaParticleIsolation.cxx

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