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