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