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( GetDebug() > 1 && leadptBin == 0 )
1200         printf("No track/clusters in isolation cone: cand pt %2.2f GeV/c, track multiplicity %d, N clusters %d\n",
1201                pt, GetTrackMultiplicity(),GetEMCALClusters()->GetEntriesFast());
1202
1203       if(IsDataMC())
1204       {
1205         Int_t leadptBinMC = leadptBin+mcIndex*fNBkgBin;
1206         Int_t  ptsumBinMC =  ptsumBin+mcIndex*fNBkgBin;
1207         if( leadptBin >=0 ) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,cluster->GetM02());
1208         if( ptsumBin  >=0 ) fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,cluster->GetM02());
1209         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1210         {
1211           leadptBinMC = leadptBin+kmcPhoton*fNBkgBin;
1212           ptsumBinMC  =  ptsumBin+kmcPhoton*fNBkgBin;
1213           if( leadptBin >=0 ) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,cluster->GetM02());
1214           if( ptsumBin  >=0 ) fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,cluster->GetM02());
1215         }
1216       }
1217     }
1218
1219
1220     if(fFillSSHisto)
1221     {
1222       fhELambda0 [isolated]->Fill(energy, cluster->GetM02() );
1223       fhPtLambda0[isolated]->Fill(pt,     cluster->GetM02() );
1224       fhELambda1 [isolated]->Fill(energy, cluster->GetM20() );
1225
1226       if(IsDataMC())
1227       {
1228         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1229           fhPtLambda0MC[kmcPhoton][isolated]->Fill(pt, cluster->GetM02());
1230
1231         fhPtLambda0MC[mcIndex][isolated]->Fill(pt, cluster->GetM02());
1232       }
1233
1234       if(fCalorimeter == "EMCAL" &&  GetFirstSMCoveredByTRD() >= 0 &&
1235          GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD()  )
1236       {
1237         fhELambda0TRD [isolated]->Fill(energy, cluster->GetM02() );
1238         fhPtLambda0TRD[isolated]->Fill(pt    , cluster->GetM02() );
1239         fhELambda1TRD [isolated]->Fill(energy, cluster->GetM20() );
1240       }
1241
1242       if(fFillNLMHistograms)
1243       {
1244         fhNLocMax[isolated]->Fill(energy,nMaxima);
1245         if     (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax1[isolated]->Fill(energy,cluster->GetM20()); }
1246         else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax2[isolated]->Fill(energy,cluster->GetM20()); }
1247         else                { fhELambda0LocMaxN[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMaxN[isolated]->Fill(energy,cluster->GetM20()); }
1248       }
1249     } // SS histo fill
1250
1251     if(fFillTMHisto)
1252     {
1253       Float_t dZ  = cluster->GetTrackDz();
1254       Float_t dR  = cluster->GetTrackDx();
1255
1256       if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1257       {
1258         dR = 2000., dZ = 2000.;
1259         GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1260       }
1261
1262       //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
1263       if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
1264       {
1265         fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
1266         fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
1267         if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
1268       }
1269
1270       // Check dEdx and E/p of matched clusters
1271
1272       if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1273       {
1274
1275         AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1276
1277         if(track)
1278         {
1279           Float_t dEdx = track->GetTPCsignal();
1280           fhdEdx[isolated]->Fill(cluster->E(), dEdx);
1281
1282           Float_t eOverp = cluster->E()/track->P();
1283           fhEOverP[isolated]->Fill(cluster->E(),  eOverp);
1284         }
1285         //else
1286         //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1287
1288
1289         if(IsDataMC())
1290         {
1291           if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
1292           {
1293             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
1294                        GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
1295             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
1296             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
1297             else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
1298
1299           }
1300           else
1301           {
1302             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
1303                        GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
1304             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
1305             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
1306             else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
1307           }
1308
1309         }  // MC
1310
1311       } // match window
1312
1313     }// TM histos fill
1314
1315   } // clusters array available
1316
1317 }
1318
1319 //______________________________________________________
1320 TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
1321 {
1322   //Save parameters used for analysis
1323   TString parList ; //this will be list of parameters used for this analysis.
1324   const Int_t buffersize = 255;
1325   char onePar[buffersize] ;
1326
1327   snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
1328   parList+=onePar ;
1329   snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1330   parList+=onePar ;
1331   snprintf(onePar, buffersize,"Isolation Cand Detector: %s\n",fIsoDetector.Data()) ;
1332   parList+=onePar ;
1333   snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
1334   parList+=onePar ;
1335   snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
1336   parList+=onePar ;
1337   snprintf(onePar, buffersize,"fFillTMHisto=%d (Flag for track matching histograms) \n",fFillTMHisto) ;
1338   parList+=onePar ;
1339   snprintf(onePar, buffersize,"fFillSSHisto=%d (Flag for shower shape histograms) \n",fFillSSHisto) ;
1340   parList+=onePar ;
1341
1342   if(fMakeSeveralIC)
1343   {
1344     snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
1345     parList+=onePar ;
1346     snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
1347     parList+=onePar ;
1348
1349     for(Int_t icone = 0; icone < fNCones ; icone++)
1350     {
1351       snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
1352       parList+=onePar ;
1353     }
1354     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1355     {
1356       snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
1357       parList+=onePar ;
1358     }
1359     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1360     {
1361       snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
1362       parList+=onePar ;
1363     }
1364     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1365     {
1366       snprintf(onePar, buffersize,"fSumPtThresholds[%d]=%1.2f (isolation sum pt threshold) \n",ipt, fSumPtThresholds[ipt]) ;
1367       parList+=onePar ;
1368     }
1369   }
1370
1371   //Get parameters set in base class.
1372   parList += GetBaseParametersList() ;
1373
1374   //Get parameters set in IC class.
1375   if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
1376
1377   return new TObjString(parList) ;
1378
1379 }
1380
1381 //________________________________________________________
1382 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
1383 {
1384   // Create histograms to be saved in output file and
1385   // store them in outputContainer
1386   TList * outputContainer = new TList() ;
1387   outputContainer->SetName("IsolatedParticleHistos") ;
1388
1389   Int_t   nptbins  = GetHistogramRanges()->GetHistoPtBins();
1390   Int_t   nphibins = GetHistogramRanges()->GetHistoPhiBins();
1391   Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
1392   Float_t ptmax    = GetHistogramRanges()->GetHistoPtMax();
1393   Float_t phimax   = GetHistogramRanges()->GetHistoPhiMax();
1394   Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
1395   Float_t ptmin    = GetHistogramRanges()->GetHistoPtMin();
1396   Float_t phimin   = GetHistogramRanges()->GetHistoPhiMin();
1397   Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
1398   Int_t   ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();
1399   Float_t ssmax    = GetHistogramRanges()->GetHistoShowerShapeMax();
1400   Float_t ssmin    = GetHistogramRanges()->GetHistoShowerShapeMin();
1401   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
1402   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
1403   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
1404
1405   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1406   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1407   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1408   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1409   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1410   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1411
1412   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
1413   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
1414   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1415   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1416   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
1417   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1418
1419   Int_t   nptsumbins    = GetHistogramRanges()->GetHistoNPtSumBins();
1420   Float_t ptsummax      = GetHistogramRanges()->GetHistoPtSumMax();
1421   Float_t ptsummin      = GetHistogramRanges()->GetHistoPtSumMin();
1422   Int_t   nptinconebins = GetHistogramRanges()->GetHistoNPtInConeBins();
1423   Float_t ptinconemax   = GetHistogramRanges()->GetHistoPtInConeMax();
1424   Float_t ptinconemin   = GetHistogramRanges()->GetHistoPtInConeMin();
1425
1426   //Float_t ptthre    = GetIsolationCut()->GetPtThreshold();
1427   //Float_t ptsumthre = GetIsolationCut()->GetSumPtThreshold();
1428   //Float_t ptfrac    = GetIsolationCut()->GetPtFraction();
1429   Float_t r         = GetIsolationCut()->GetConeSize();
1430   Int_t   method    = GetIsolationCut()->GetICMethod() ;
1431   Int_t   particle  = GetIsolationCut()->GetParticleTypeInCone() ;
1432
1433   TString sThreshold = "";
1434   if      ( method == AliIsolationCut::kSumPtIC )
1435   {
1436     sThreshold = Form(", %2.2f < #Sigma #it{p}_{T}^{in cone} < %2.2f GeV/#it{c}",
1437                       GetIsolationCut()->GetSumPtThreshold(), GetIsolationCut()->GetSumPtThresholdMax());
1438     if(GetIsolationCut()->GetSumPtThresholdMax() > 200)
1439       sThreshold = Form(", #Sigma #it{p}_{T}^{in cone} = %2.2f GeV/#it{c}",
1440                         GetIsolationCut()->GetSumPtThreshold());
1441   }
1442   else if ( method == AliIsolationCut::kPtThresIC)
1443   {
1444     sThreshold = Form(", %2.2f < #it{p}_{T}^{th} < %2.2f GeV/#it{c}",
1445                       GetIsolationCut()->GetPtThreshold(),GetIsolationCut()->GetPtThresholdMax());
1446     if(GetIsolationCut()->GetSumPtThreshold() > 200)
1447       sThreshold = Form(", #it{p}_{T}^{th} = %2.2f GeV/#it{c}",
1448                         GetIsolationCut()->GetPtThreshold());
1449   }
1450   else if ( method == AliIsolationCut::kPtFracIC)
1451     sThreshold = Form(", #Sigma #it{p}_{T}^{in cone}/#it{p}_{T}^{trig} = %2.2f" ,
1452                       GetIsolationCut()->GetPtFraction());
1453
1454   TString sParticle = ", x^{0,#pm}";
1455   if      ( particle == AliIsolationCut::kOnlyNeutral )  sParticle = ", x^{0}";
1456   else if ( particle == AliIsolationCut::kOnlyCharged )  sParticle = ", x^{#pm}";
1457
1458   TString parTitle = Form("#it{R} = %2.2f%s%s",GetIsolationCut()->GetConeSize(), sThreshold.Data(),sParticle.Data());
1459
1460   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1461
1462   // MC histograms title and name
1463   TString mcPartType[] = { "#gamma", "#gamma_{prompt}", "#gamma_{fragmentation}",
1464     "#pi^{0} (merged #gamma)","#gamma_{#pi decay}",
1465     "#gamma_{#eta decay}","#gamma_{other decay}",
1466     "e^{#pm}","hadrons?"} ;
1467
1468   TString mcPartName[] = { "Photon","PhotonPrompt","PhotonFrag",
1469     "Pi0","Pi0Decay","EtaDecay","OtherDecay",
1470     "Electron","Hadron"} ;
1471
1472   // Primary MC histograms title and name
1473   TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}",
1474     "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
1475
1476   TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay",
1477     "PhotonPrompt","PhotonFrag","PhotonISR"} ;
1478
1479   // Not Isolated histograms, reference histograms
1480
1481   fhENoIso  = new TH1F("hENoIso",
1482                        Form("Number of not isolated leading particles vs #it{p}_{T}, %s",parTitle.Data()),
1483                        nptbins,ptmin,ptmax);
1484   fhENoIso->SetYTitle("#it{counts}");
1485   fhENoIso->SetXTitle("E (GeV/#it{c})");
1486   outputContainer->Add(fhENoIso) ;
1487
1488   fhPtNoIso  = new TH1F("hPtNoIso",
1489                         Form("Number of not isolated leading particles vs #it{p}_{T}, %s",parTitle.Data()),
1490                         nptbins,ptmin,ptmax);
1491   fhPtNoIso->SetYTitle("#it{counts}");
1492   fhPtNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1493   outputContainer->Add(fhPtNoIso) ;
1494
1495   fhEtaPhiNoIso  = new TH2F("hEtaPhiNoIso",
1496                             Form("Number of not isolated leading particles #eta vs #phi, %s",parTitle.Data()),
1497                             netabins,etamin,etamax,nphibins,phimin,phimax);
1498   fhEtaPhiNoIso->SetXTitle("#eta");
1499   fhEtaPhiNoIso->SetYTitle("#phi");
1500   outputContainer->Add(fhEtaPhiNoIso) ;
1501
1502   if(IsDataMC())
1503   {
1504     // For histograms in arrays, index in the array, corresponding to any particle origin
1505
1506     for(Int_t imc = 0; imc < 9; imc++)
1507     {
1508
1509       fhPtNoIsoMC[imc]  = new TH1F(Form("hPtNoIsoMC%s",mcPartName[imc].Data()),
1510                                    Form("#it{p}_{T} of NOT isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1511                                    nptbins,ptmin,ptmax);
1512       fhPtNoIsoMC[imc]->SetYTitle("#it{counts}");
1513       fhPtNoIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1514       outputContainer->Add(fhPtNoIsoMC[imc]) ;
1515
1516       fhPtIsoMC[imc]  = new TH1F(Form("hPtMC%s",mcPartName[imc].Data()),
1517                                  Form("#it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1518                                  nptbins,ptmin,ptmax);
1519       fhPtIsoMC[imc]->SetYTitle("#it{counts}");
1520       fhPtIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1521       outputContainer->Add(fhPtIsoMC[imc]) ;
1522
1523       fhPhiIsoMC[imc]  = new TH2F(Form("hPhiMC%s",mcPartName[imc].Data()),
1524                                   Form("#phi vs #it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1525                                   nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1526       fhPhiIsoMC[imc]->SetYTitle("#phi");
1527       fhPhiIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1528       outputContainer->Add(fhPhiIsoMC[imc]) ;
1529
1530       fhEtaIsoMC[imc]  = new TH2F(Form("hEtaMC%s",mcPartName[imc].Data()),
1531                                   Form("#phi vs #it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1532                                   nptbins,ptmin,ptmax,netabins,etamin,etamax);
1533       fhEtaIsoMC[imc]->SetYTitle("#eta");
1534       fhEtaIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1535       outputContainer->Add(fhEtaIsoMC[imc]) ;
1536     }
1537   }
1538
1539   // Histograms for tagged candidates as decay
1540   if(fFillTaggedDecayHistograms)
1541   {
1542     for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1543     {
1544       fhPtDecayNoIso[ibit]  =
1545       new TH1F(Form("hPtDecayNoIso_bit%d",fDecayBits[ibit]),
1546                Form("Number of not isolated leading pi0 decay particles vs #it{p}_{T}, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1547                nptbins,ptmin,ptmax);
1548       fhPtDecayNoIso[ibit]->SetYTitle("#it{counts}");
1549       fhPtDecayNoIso[ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1550       outputContainer->Add(fhPtDecayNoIso[ibit]) ;
1551
1552       fhEtaPhiDecayNoIso[ibit]  =
1553       new TH2F(Form("hEtaPhiDecayNoIso_bit%d",fDecayBits[ibit]),
1554                Form("Number of not isolated leading Pi0 decay particles #eta vs #phi, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1555                netabins,etamin,etamax,nphibins,phimin,phimax);
1556       fhEtaPhiDecayNoIso[ibit]->SetXTitle("#eta");
1557       fhEtaPhiDecayNoIso[ibit]->SetYTitle("#phi");
1558       outputContainer->Add(fhEtaPhiDecayNoIso[ibit]) ;
1559
1560       if(!fMakeSeveralIC)
1561       {
1562         fhPtDecayIso[ibit]  =
1563         new TH1F(Form("hPtDecayIso_bit%d",fDecayBits[ibit]),
1564                  Form("Number of isolated #pi^{0} decay particles vs #it{p}_{T}, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1565                  nptbins,ptmin,ptmax);
1566         fhPtDecayIso[ibit]->SetYTitle("#it{counts}");
1567         fhPtDecayIso[ibit]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1568         outputContainer->Add(fhPtDecayIso[ibit]) ;
1569
1570         fhEtaPhiDecayIso[ibit]  =
1571         new TH2F(Form("hEtaPhiDecayIso_bit%d",fDecayBits[ibit]),
1572                  Form("Number of isolated Pi0 decay particles #eta vs #phi, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1573                  netabins,etamin,etamax,nphibins,phimin,phimax);
1574         fhEtaPhiDecayIso[ibit]->SetXTitle("#eta");
1575         fhEtaPhiDecayIso[ibit]->SetYTitle("#phi");
1576         outputContainer->Add(fhEtaPhiDecayIso[ibit]) ;
1577       }
1578
1579       if(IsDataMC())
1580       {
1581         for(Int_t imc = 0; imc < 9; imc++)
1582         {
1583
1584           fhPtDecayNoIsoMC[ibit][imc]  =
1585           new TH1F(Form("hPtDecayNoIso_bit%d_MC%s",fDecayBits[ibit],mcPartName[imc].Data()),
1586                    Form("#it{p}_{T} of NOT isolated, decay bit %d,  %s, %s",fDecayBits[ibit],mcPartType[imc].Data(),parTitle.Data()),
1587                    nptbins,ptmin,ptmax);
1588           fhPtDecayNoIsoMC[ibit][imc]->SetYTitle("#it{counts}");
1589           fhPtDecayNoIsoMC[ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1590           outputContainer->Add(fhPtDecayNoIsoMC[ibit][imc]) ;
1591
1592           if(!fMakeSeveralIC)
1593           {
1594             fhPtDecayIsoMC[ibit][imc]  =
1595             new TH1F(Form("hPtDecay_bit%d_MC%s",fDecayBits[ibit],mcPartName[imc].Data()),
1596                      Form("#it{p}_{T} of isolated %s, decay bit %d, %s",mcPartType[imc].Data(),fDecayBits[ibit],parTitle.Data()),
1597                      nptbins,ptmin,ptmax);
1598             fhPtDecayIsoMC[ibit][imc]->SetYTitle("#it{counts}");
1599             fhPtDecayIsoMC[ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1600             outputContainer->Add(fhPtDecayIsoMC[ibit][imc]) ;
1601           }
1602         }// MC particle loop
1603       }// MC
1604     } // bit loop
1605   }// decay
1606
1607   if(!fMakeSeveralIC)
1608   {
1609     TString isoName [] = {"NoIso",""};
1610     TString isoTitle[] = {"Not isolated"  ,"isolated"};
1611
1612     fhEIso   = new TH1F("hE",
1613                         Form("Number of isolated particles vs E, %s",parTitle.Data()),
1614                         nptbins,ptmin,ptmax);
1615     fhEIso->SetYTitle("d#it{N} / d#it{E}");
1616     fhEIso->SetXTitle("#it{E} (GeV/#it{c})");
1617     outputContainer->Add(fhEIso) ;
1618
1619     fhPtIso  = new TH1F("hPt",
1620                         Form("Number of isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1621                         nptbins,ptmin,ptmax);
1622     fhPtIso->SetYTitle("d#it{N} / #it{p}_{T}");
1623     fhPtIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1624     outputContainer->Add(fhPtIso) ;
1625
1626     fhPhiIso  = new TH2F("hPhi",
1627                          Form("Number of isolated particles vs #phi, %s",parTitle.Data()),
1628                          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1629     fhPhiIso->SetYTitle("#phi");
1630     fhPhiIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1631     outputContainer->Add(fhPhiIso) ;
1632
1633     fhEtaIso  = new TH2F("hEta",
1634                          Form("Number of isolated particles vs #eta, %s",parTitle.Data()),
1635                          nptbins,ptmin,ptmax,netabins,etamin,etamax);
1636     fhEtaIso->SetYTitle("#eta");
1637     fhEtaIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1638     outputContainer->Add(fhEtaIso) ;
1639
1640     fhEtaPhiIso  = new TH2F("hEtaPhiIso",
1641                             Form("Number of isolated particles #eta vs #phi, %s",parTitle.Data()),
1642                             netabins,etamin,etamax,nphibins,phimin,phimax);
1643     fhEtaPhiIso->SetXTitle("#eta");
1644     fhEtaPhiIso->SetYTitle("#phi");
1645     outputContainer->Add(fhEtaPhiIso) ;
1646
1647     if(fFillHighMultHistograms)
1648     {
1649       fhPtCentralityIso  = new TH2F("hPtCentrality",
1650                                     Form("centrality vs #it{p}_{T} for isolated particles, %s",parTitle.Data()),
1651                                     nptbins,ptmin,ptmax, 100,0,100);
1652       fhPtCentralityIso->SetYTitle("centrality");
1653       fhPtCentralityIso->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1654       outputContainer->Add(fhPtCentralityIso) ;
1655
1656       fhPtEventPlaneIso  = new TH2F("hPtEventPlane",
1657                                     Form("event plane angle vs #it{p}_{T} for isolated particles, %s",parTitle.Data()),
1658                                     nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1659       fhPtEventPlaneIso->SetYTitle("Event plane angle (rad)");
1660       fhPtEventPlaneIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1661       outputContainer->Add(fhPtEventPlaneIso) ;
1662     }
1663
1664     if(fFillNLMHistograms)
1665     {
1666       fhPtNLocMaxIso  = new TH2F("hPtNLocMax",
1667                                  Form("Number of isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1668                                  nptbins,ptmin,ptmax,10,0,10);
1669       fhPtNLocMaxIso->SetYTitle("#it{NLM}");
1670       fhPtNLocMaxIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1671
1672       fhPtNLocMaxNoIso  = new TH2F("hPtNLocMaxNoIso",
1673                                    Form("Number of not isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1674                                    nptbins,ptmin,ptmax,10,0,10);
1675       fhPtNLocMaxNoIso->SetYTitle("#it{NLM}");
1676       fhPtNLocMaxNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1677       outputContainer->Add(fhPtNLocMaxNoIso) ;
1678     }
1679
1680     fhConePtLead  = new TH2F("hConePtLead",
1681                             Form("Track or Cluster  leading #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1682                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1683     fhConePtLead->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
1684     fhConePtLead->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1685     outputContainer->Add(fhConePtLead) ;
1686
1687
1688     fhConeSumPt  = new TH2F("hConePtSum",
1689                             Form("Track and Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1690                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1691     fhConeSumPt->SetYTitle("#Sigma #it{p}_{T}");
1692     fhConeSumPt->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1693     outputContainer->Add(fhConeSumPt) ;
1694
1695     fhConeSumPtTrigEtaPhi  = new TH2F("hConePtSumTrigEtaPhi",
1696                                       Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1697                                       netabins,etamin,etamax,nphibins,phimin,phimax);
1698     fhConeSumPtTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1699     fhConeSumPtTrigEtaPhi->SetXTitle("#eta_{trigger}");
1700     fhConeSumPtTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1701     outputContainer->Add(fhConeSumPtTrigEtaPhi) ;
1702
1703     fhPtInCone  = new TH2F("hPtInCone",
1704                            Form("#it{p}_{T} of clusters and tracks in isolation cone for #it{R} =  %2.2f",r),
1705                            nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1706     fhPtInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1707     fhPtInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1708     outputContainer->Add(fhPtInCone) ;
1709
1710     if(fFillBackgroundBinHistograms)
1711     {
1712       fhPtLeadConeBinLambda0 = new TH2F*[fNBkgBin];
1713       fhSumPtConeBinLambda0  = new TH2F*[fNBkgBin];
1714
1715       if(IsDataMC())
1716       {
1717         fhPtLeadConeBinLambda0MC = new TH2F*[fNBkgBin*9];
1718         fhSumPtConeBinLambda0MC  = new TH2F*[fNBkgBin*9];
1719       }
1720
1721       for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1722       {
1723         fhPtLeadConeBinLambda0[ibin]  = new TH2F
1724         (Form("hPtLeadConeLambda0_Bin%d",ibin),
1725          Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f (GeV/#it{c}), %s",
1726               fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1727         fhPtLeadConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
1728         fhPtLeadConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1729         outputContainer->Add(fhPtLeadConeBinLambda0[ibin]) ;
1730
1731         fhSumPtConeBinLambda0[ibin]  = new TH2F
1732         (Form("hSumPtConeLambda0_Bin%d",ibin),
1733          Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f (GeV/#it{c}), %s",
1734               fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1735         fhSumPtConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
1736         fhSumPtConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1737         outputContainer->Add(fhSumPtConeBinLambda0[ibin]) ;
1738
1739         if(IsDataMC())
1740         {
1741           for(Int_t imc = 0; imc < 9; imc++)
1742           {
1743             Int_t binmc = ibin+imc*fNBkgBin;
1744             fhPtLeadConeBinLambda0MC[binmc]  = new TH2F
1745             (Form("hPtLeadConeLambda0_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
1746              Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f (GeV/#it{c}), MC %s, %s",
1747                   fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1748             fhPtLeadConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
1749             fhPtLeadConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1750             outputContainer->Add(fhPtLeadConeBinLambda0MC[binmc]) ;
1751
1752             fhSumPtConeBinLambda0MC[binmc]  = new TH2F
1753             (Form("hSumPtConeLambda0_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
1754              Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f (GeV/#it{c}), MC %s, %s",
1755                   fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1756             fhSumPtConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
1757             fhSumPtConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1758             outputContainer->Add(fhSumPtConeBinLambda0MC[binmc]) ;
1759           } // MC particle loop
1760         }
1761
1762       }// pt bin loop
1763     } // bkg cone pt bin histograms
1764
1765     if(fFillHighMultHistograms)
1766     {
1767       fhPtInConeCent  = new TH2F("hPtInConeCent",
1768                                  Form("#it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1769                                  100,0,100,nptinconebins,ptinconemin,ptinconemax);
1770       fhPtInConeCent->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1771       fhPtInConeCent->SetXTitle("centrality");
1772       outputContainer->Add(fhPtInConeCent) ;
1773     }
1774
1775     // Cluster only histograms
1776     if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyCharged)
1777     {
1778       fhConeSumPtCluster  = new TH2F("hConePtSumCluster",
1779                                      Form("Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1780                                      nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1781       fhConeSumPtCluster->SetYTitle("#Sigma #it{p}_{T}");
1782       fhConeSumPtCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1783       outputContainer->Add(fhConeSumPtCluster) ;
1784
1785       fhConePtLeadCluster  = new TH2F("hConeLeadPtCluster",
1786                                     Form("Cluster leading in isolation cone for #it{R} =  %2.2f",r),
1787                                     nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1788       fhConePtLeadCluster->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
1789       fhConePtLeadCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1790       outputContainer->Add(fhConePtLeadCluster) ;
1791
1792
1793       if(fFillCellHistograms)
1794       {
1795         fhConeSumPtCell  = new TH2F("hConePtSumCell",
1796                                     Form("Cell #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
1797                                     nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1798         fhConeSumPtCell->SetYTitle("#Sigma #it{p}_{T}");
1799         fhConeSumPtCell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1800         outputContainer->Add(fhConeSumPtCell) ;
1801       }
1802
1803       if(fFillUEBandSubtractHistograms)
1804       {
1805         fhConeSumPtEtaBandUECluster  = new TH2F("hConePtSumEtaBandUECluster",
1806                                                 "#Sigma cluster #it{p}_{T} in UE Eta Band",
1807                                                 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1808         fhConeSumPtEtaBandUECluster->SetYTitle("#Sigma #it{p}_{T}");
1809         fhConeSumPtEtaBandUECluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1810         outputContainer->Add(fhConeSumPtEtaBandUECluster) ;
1811
1812         fhConeSumPtPhiBandUECluster  = new TH2F("hConePtSumPhiBandUECluster",
1813                                                 "#Sigma cluster #it{p}_{T} UE Phi Band",
1814                                                 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1815         fhConeSumPtPhiBandUECluster->SetYTitle("#Sigma #it{p}_{T}");
1816         fhConeSumPtPhiBandUECluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1817         outputContainer->Add(fhConeSumPtPhiBandUECluster) ;
1818
1819         fhConeSumPtEtaBandUEClusterTrigEtaPhi  = new TH2F("hConePtSumEtaBandUEClusterTrigEtaPhi",
1820                                                           "Trigger #eta vs #phi, #Sigma cluster #it{p}_{T} in UE Eta Band",
1821                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
1822         fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1823         fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1824         fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1825         outputContainer->Add(fhConeSumPtEtaBandUEClusterTrigEtaPhi) ;
1826
1827         fhConeSumPtPhiBandUEClusterTrigEtaPhi  = new TH2F("hConePtSumPhiBandUEClusterTrigEtaPhi",
1828                                                           "Trigger #eta vs #phi, #Sigma cluster #it{p}_{T} UE Phi Band",
1829                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
1830         fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1831         fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1832         fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1833         outputContainer->Add(fhConeSumPtPhiBandUEClusterTrigEtaPhi) ;
1834         if(fFillCellHistograms)
1835         {
1836
1837           fhConeSumPtEtaBandUECell  = new TH2F("hConePtSumEtaBandUECell",
1838                                                "#Sigma cell #it{p}_{T} in UE Eta Band",
1839                                                nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1840           fhConeSumPtEtaBandUECell->SetYTitle("#Sigma #it{p}_{T}");
1841           fhConeSumPtEtaBandUECell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1842           outputContainer->Add(fhConeSumPtEtaBandUECell) ;
1843
1844           fhConeSumPtPhiBandUECell  = new TH2F("hConePtSumPhiBandUECell",
1845                                                "#Sigma cell #it{p}_{T} UE Phi Band",
1846                                                nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1847           fhConeSumPtPhiBandUECell->SetYTitle("#Sigma #it{p}_{T}");
1848           fhConeSumPtPhiBandUECell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1849           outputContainer->Add(fhConeSumPtPhiBandUECell) ;
1850
1851           fhConeSumPtEtaBandUECellTrigEtaPhi  = new TH2F("hConePtSumEtaBandUECellTrigEtaPhi",
1852                                                          "Trigger #eta vs #phi, #Sigma cell #it{p}_{T} in UE Eta Band",
1853                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
1854           fhConeSumPtEtaBandUECellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1855           fhConeSumPtEtaBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1856           fhConeSumPtEtaBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1857           outputContainer->Add(fhConeSumPtEtaBandUECellTrigEtaPhi) ;
1858
1859           fhConeSumPtPhiBandUECellTrigEtaPhi  = new TH2F("hConePtSumPhiBandUECellTrigEtaPhi",
1860                                                          "Trigger #eta vs #phi, #Sigma cell #it{p}_{T} UE Phi Band",
1861                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
1862           fhConeSumPtPhiBandUECellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1863           fhConeSumPtPhiBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1864           fhConeSumPtPhiBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1865           outputContainer->Add(fhConeSumPtPhiBandUECellTrigEtaPhi) ;
1866         }
1867
1868         fhEtaBandCluster  = new TH2F("hEtaBandCluster",
1869                                      Form("#eta vs #phi of clusters in #eta band isolation cone for #it{R} =  %2.2f",r),
1870                                      netabins,-1,1,nphibins,0,TMath::TwoPi());
1871         fhEtaBandCluster->SetXTitle("#eta");
1872         fhEtaBandCluster->SetYTitle("#phi");
1873         outputContainer->Add(fhEtaBandCluster) ;
1874
1875         fhPhiBandCluster  = new TH2F("hPhiBandCluster",
1876                                      Form("#eta vs #phi of clusters in #phi band isolation cone for #it{R} =  %2.2f",r),
1877                                      netabins,-1,1,nphibins,0,TMath::TwoPi());
1878         fhPhiBandCluster->SetXTitle("#eta");
1879         fhPhiBandCluster->SetYTitle("#phi");
1880         outputContainer->Add(fhPhiBandCluster) ;
1881
1882         fhEtaPhiInConeCluster= new TH2F("hEtaPhiInConeCluster",
1883                                         Form("#eta vs #phi of clusters in cone for #it{R} =  %2.2f",r),
1884                                         netabins,-1,1,nphibins,0,TMath::TwoPi());
1885         fhEtaPhiInConeCluster->SetXTitle("#eta");
1886         fhEtaPhiInConeCluster->SetYTitle("#phi");
1887         outputContainer->Add(fhEtaPhiInConeCluster) ;
1888
1889         fhEtaPhiCluster= new TH2F("hEtaPhiCluster",
1890                                   Form("#eta vs #phi of all clusters"),
1891                                   netabins,-1,1,nphibins,0,TMath::TwoPi());
1892         fhEtaPhiCluster->SetXTitle("#eta");
1893         fhEtaPhiCluster->SetYTitle("#phi");
1894         outputContainer->Add(fhEtaPhiCluster) ;
1895
1896       }
1897
1898       fhPtClusterInCone  = new TH2F("hPtClusterInCone",
1899                                     Form("#it{p}_{T} of clusters in isolation cone for #it{R} =  %2.2f",r),
1900                                     nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1901       fhPtClusterInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1902       fhPtClusterInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1903       outputContainer->Add(fhPtClusterInCone) ;
1904
1905       if(fFillCellHistograms)
1906       {
1907         fhPtCellInCone  = new TH2F("hPtCellInCone",
1908                                    Form("#it{p}_{T} of cells in isolation cone for #it{R} =  %2.2f",r),
1909                                    nptbins,ptmin,ptmax,1000,0,50);
1910         fhPtCellInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1911         fhPtCellInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1912         outputContainer->Add(fhPtCellInCone) ;
1913
1914         fhEtaBandCell  = new TH2F("hEtaBandCell",
1915                                   Form("#col vs #row of cells in #eta band isolation cone for #it{R} =  %2.2f",r),
1916                                   96,0,95,128,0,127);
1917         fhEtaBandCell->SetXTitle("#col");
1918         fhEtaBandCell->SetYTitle("#row");
1919         outputContainer->Add(fhEtaBandCell) ;
1920
1921         fhPhiBandCell  = new TH2F("hPhiBandCell",
1922                                   Form("#col vs #row of cells in #phi band isolation cone for #it{R} =  %2.2f",r),
1923                                   96,0,95,128,0,127);
1924         fhPhiBandCell->SetXTitle("#col");
1925         fhPhiBandCell->SetYTitle("#row");
1926         outputContainer->Add(fhPhiBandCell) ;
1927       }
1928
1929       if(fFillUEBandSubtractHistograms)
1930       {
1931         fhConeSumPtEtaUESubCluster  = new TH2F("hConeSumPtEtaUESubCluster",
1932                                                Form("Clusters #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
1933                                                nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1934         fhConeSumPtEtaUESubCluster->SetYTitle("#Sigma #it{p}_{T}");
1935         fhConeSumPtEtaUESubCluster->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1936         outputContainer->Add(fhConeSumPtEtaUESubCluster) ;
1937
1938         fhConeSumPtPhiUESubCluster  = new TH2F("hConeSumPtPhiUESubCluster",
1939                                                Form("Clusters #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
1940                                                nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1941         fhConeSumPtPhiUESubCluster->SetYTitle("#Sigma #it{p}_{T}");
1942         fhConeSumPtPhiUESubCluster->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1943         outputContainer->Add(fhConeSumPtPhiUESubCluster) ;
1944
1945         fhConeSumPtEtaUESubClusterTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubClusterTrigEtaPhi",
1946                                                          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),
1947                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
1948         fhConeSumPtEtaUESubClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1949         fhConeSumPtEtaUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1950         fhConeSumPtEtaUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1951         outputContainer->Add(fhConeSumPtEtaUESubClusterTrigEtaPhi) ;
1952
1953         fhConeSumPtPhiUESubClusterTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubClusterTrigEtaPhi",
1954                                                          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),
1955                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
1956         fhConeSumPtPhiUESubClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1957         fhConeSumPtPhiUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1958         fhConeSumPtPhiUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1959         outputContainer->Add(fhConeSumPtPhiUESubClusterTrigEtaPhi) ;
1960
1961         if(fFillCellHistograms)
1962         {
1963           fhConeSumPtEtaUESubCell  = new TH2F("hConeSumPtEtaUESubCell",
1964                                               Form("Cells #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
1965                                               nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1966           fhConeSumPtEtaUESubCell->SetYTitle("#Sigma #it{p}_{T}");
1967           fhConeSumPtEtaUESubCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1968           outputContainer->Add(fhConeSumPtEtaUESubCell) ;
1969
1970           fhConeSumPtPhiUESubCell  = new TH2F("hConeSumPtPhiUESubCell",
1971                                               Form("Cells #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
1972                                               nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1973           fhConeSumPtPhiUESubCell->SetYTitle("#Sigma #it{p}_{T}");
1974           fhConeSumPtPhiUESubCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1975           outputContainer->Add(fhConeSumPtPhiUESubCell) ;
1976
1977           fhConeSumPtEtaUESubCellTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubCellTrigEtaPhi",
1978                                                         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),
1979                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
1980           fhConeSumPtEtaUESubCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1981           fhConeSumPtEtaUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1982           fhConeSumPtEtaUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1983           outputContainer->Add(fhConeSumPtEtaUESubCellTrigEtaPhi) ;
1984
1985           fhConeSumPtPhiUESubCellTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubCellTrigEtaPhi",
1986                                                         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),
1987                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
1988           fhConeSumPtPhiUESubCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1989           fhConeSumPtPhiUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1990           fhConeSumPtPhiUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1991           outputContainer->Add(fhConeSumPtPhiUESubCellTrigEtaPhi) ;
1992         }
1993
1994         fhFractionClusterOutConeEta  = new TH2F("hFractionClusterOutConeEta",
1995                                                 Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #eta acceptance",r),
1996                                                 nptbins,ptmin,ptmax,100,0,1);
1997         fhFractionClusterOutConeEta->SetYTitle("#it{fraction}");
1998         fhFractionClusterOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
1999         outputContainer->Add(fhFractionClusterOutConeEta) ;
2000
2001         fhFractionClusterOutConeEtaTrigEtaPhi  = new TH2F("hFractionClusterOutConeEtaTrigEtaPhi",
2002                                                           Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #eta acceptance, in trigger #eta-#phi ",r),
2003                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
2004         fhFractionClusterOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2005         fhFractionClusterOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2006         fhFractionClusterOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2007         outputContainer->Add(fhFractionClusterOutConeEtaTrigEtaPhi) ;
2008
2009         fhFractionClusterOutConePhi  = new TH2F("hFractionClusterOutConePhi",
2010                                                 Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #phi acceptance",r),
2011                                                 nptbins,ptmin,ptmax,100,0,1);
2012         fhFractionClusterOutConePhi->SetYTitle("#it{fraction}");
2013         fhFractionClusterOutConePhi->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2014         outputContainer->Add(fhFractionClusterOutConePhi) ;
2015
2016         fhFractionClusterOutConePhiTrigEtaPhi  = new TH2F("hFractionClusterOutConePhiTrigEtaPhi",
2017                                                           Form("Fraction of the isolation cone #it{R} =  %2.2f, out of clusters #phi acceptance, in trigger #eta-#phi ",r),
2018                                                           netabins,etamin,etamax,nphibins,phimin,phimax);
2019         fhFractionClusterOutConePhiTrigEtaPhi->SetZTitle("#it{fraction}");
2020         fhFractionClusterOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
2021         fhFractionClusterOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2022         outputContainer->Add(fhFractionClusterOutConePhiTrigEtaPhi) ;
2023
2024         fhConeSumPtSubvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCluster",
2025                                                           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),
2026                                                           nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2027         fhConeSumPtSubvsConeSumPtTotPhiCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2028         fhConeSumPtSubvsConeSumPtTotPhiCluster->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2029         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCluster);
2030
2031         fhConeSumPtSubNormvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCluster",
2032                                                               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),
2033                                                               nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2034         fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2035         fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2036         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCluster);
2037
2038         fhConeSumPtSubvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCluster",
2039                                                           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),
2040                                                           nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2041         fhConeSumPtSubvsConeSumPtTotEtaCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2042         fhConeSumPtSubvsConeSumPtTotEtaCluster->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2043         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCluster);
2044
2045         fhConeSumPtSubNormvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCluster",
2046                                                               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),
2047                                                               nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2048         fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2049         fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2050         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCluster);
2051
2052         fhConeSumPtVSUEClusterEtaBand  = new TH2F("hConeSumPtVSUEClusterEtaBand",
2053                                                   Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in eta band for cluster (before normalization), R=%2.2f",r),
2054                                                   nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
2055         fhConeSumPtVSUEClusterEtaBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2056         fhConeSumPtVSUEClusterEtaBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2057         outputContainer->Add(fhConeSumPtVSUEClusterEtaBand);
2058
2059         fhConeSumPtVSUEClusterPhiBand  = new TH2F("hConeSumPtVSUEClusterPhiBand",
2060                                                   Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in phi band for cluster (before normalization), R=%2.2f",r),
2061                                                   nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
2062         fhConeSumPtVSUEClusterPhiBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2063         fhConeSumPtVSUEClusterPhiBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2064         outputContainer->Add(fhConeSumPtVSUEClusterPhiBand);
2065
2066         if(fFillCellHistograms)
2067         {
2068           fhFractionCellOutConeEta  = new TH2F("hFractionCellOutConeEta",
2069                                                Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #eta acceptance",r),
2070                                                nptbins,ptmin,ptmax,100,0,1);
2071           fhFractionCellOutConeEta->SetYTitle("#it{fraction}");
2072           fhFractionCellOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2073           outputContainer->Add(fhFractionCellOutConeEta) ;
2074
2075           fhFractionCellOutConeEtaTrigEtaPhi  = new TH2F("hFractionCellOutConeEtaTrigEtaPhi",
2076                                                          Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #eta acceptance, in trigger #eta-#phi ",r),
2077                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2078           fhFractionCellOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2079           fhFractionCellOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2080           fhFractionCellOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2081           outputContainer->Add(fhFractionCellOutConeEtaTrigEtaPhi) ;
2082
2083           fhFractionCellOutConePhi  = new TH2F("hFractionCellOutConePhi",
2084                                                Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #phi acceptance",r),
2085                                                nptbins,ptmin,ptmax,100,0,1);
2086           fhFractionCellOutConePhi->SetYTitle("#it{fraction}");
2087           fhFractionCellOutConePhi->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2088           outputContainer->Add(fhFractionCellOutConePhi) ;
2089
2090           fhFractionCellOutConePhiTrigEtaPhi  = new TH2F("hFractionCellOutConePhiTrigEtaPhi",
2091                                                          Form("Fraction of the isolation cone #it{R} =  %2.2f, out of cells #phi acceptance, in trigger #eta-#phi ",r),
2092                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2093           fhFractionCellOutConePhiTrigEtaPhi->SetZTitle("#it{fraction}");
2094           fhFractionCellOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
2095           fhFractionCellOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2096           outputContainer->Add(fhFractionCellOutConePhiTrigEtaPhi) ;
2097
2098
2099           fhConeSumPtSubvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCell",
2100                                                          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),
2101                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2102           fhConeSumPtSubvsConeSumPtTotPhiCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2103           fhConeSumPtSubvsConeSumPtTotPhiCell->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2104           outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCell);
2105
2106           fhConeSumPtSubNormvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCell",
2107                                                              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),
2108                                                              nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2109           fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2110           fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2111           outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCell);
2112
2113           fhConeSumPtSubvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCell",
2114                                                          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),
2115                                                          nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2116           fhConeSumPtSubvsConeSumPtTotEtaCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2117           fhConeSumPtSubvsConeSumPtTotEtaCell->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2118           outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCell);
2119
2120           fhConeSumPtSubNormvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCell",
2121                                                              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),
2122                                                              nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2123           fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2124           fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2125           outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCell);
2126         }
2127       }
2128     }
2129
2130     // Track only histograms
2131     if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyNeutral)
2132     {
2133       fhConeSumPtTrack  = new TH2F("hConePtSumTrack",
2134                                    Form("Track #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2135                                    nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2136       fhConeSumPtTrack->SetYTitle("#Sigma #it{p}_{T}");
2137       fhConeSumPtTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2138       outputContainer->Add(fhConeSumPtTrack) ;
2139
2140       fhConePtLeadTrack  = new TH2F("hConeLeadPtTrack",
2141                                    Form("Track leading in isolation cone for #it{R} =  %2.2f",r),
2142                                    nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2143       fhConePtLeadTrack->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
2144       fhConePtLeadTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2145       outputContainer->Add(fhConePtLeadTrack) ;
2146
2147       fhPtTrackInCone  = new TH2F("hPtTrackInCone",
2148                                   Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f",r),
2149                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2150       fhPtTrackInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2151       fhPtTrackInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2152       outputContainer->Add(fhPtTrackInCone) ;
2153
2154
2155       if(fFillUEBandSubtractHistograms)
2156       {
2157         fhConeSumPtEtaBandUETrack  = new TH2F("hConePtSumEtaBandUETrack",
2158                                               "#Sigma track #it{p}_{T} in UE Eta Band",
2159                                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2160         fhConeSumPtEtaBandUETrack->SetYTitle("#Sigma #it{p}_{T}");
2161         fhConeSumPtEtaBandUETrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2162         outputContainer->Add(fhConeSumPtEtaBandUETrack) ;
2163
2164         fhConeSumPtPhiBandUETrack  = new TH2F("hConePtSumPhiBandUETrack",
2165                                               "#Sigma track #it{p}_{T} in UE Phi Band",
2166                                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax*8);
2167         fhConeSumPtPhiBandUETrack->SetYTitle("#Sigma #it{p}_{T}");
2168         fhConeSumPtPhiBandUETrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2169         outputContainer->Add(fhConeSumPtPhiBandUETrack) ;
2170
2171
2172         fhConeSumPtEtaBandUETrackTrigEtaPhi  = new TH2F("hConePtSumEtaBandUETrackTrigEtaPhi",
2173                                                         "Trigger #eta vs #phi, #Sigma track #it{p}_{T} in UE Eta Band",
2174                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2175         fhConeSumPtEtaBandUETrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2176         fhConeSumPtEtaBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2177         fhConeSumPtEtaBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2178         outputContainer->Add(fhConeSumPtEtaBandUETrackTrigEtaPhi) ;
2179
2180         fhConeSumPtPhiBandUETrackTrigEtaPhi  = new TH2F("hConePtSumPhiBandUETrackTrigEtaPhi",
2181                                                         "Trigger #eta vs #phi, #Sigma track #it{p}_{T} in UE Phi Band",
2182                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2183         fhConeSumPtPhiBandUETrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2184         fhConeSumPtPhiBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2185         fhConeSumPtPhiBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2186         outputContainer->Add(fhConeSumPtPhiBandUETrackTrigEtaPhi) ;
2187
2188         fhEtaBandTrack  = new TH2F("hEtaBandTrack",
2189                                    Form("#eta vs #phi of tracks in #eta band isolation cone for #it{R} =  %2.2f",r),
2190                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
2191         fhEtaBandTrack->SetXTitle("#eta");
2192         fhEtaBandTrack->SetYTitle("#phi");
2193         outputContainer->Add(fhEtaBandTrack) ;
2194
2195         fhPhiBandTrack  = new TH2F("hPhiBandTrack",
2196                                    Form("#eta vs #phi of tracks in #phi band isolation cone for #it{R} =  %2.2f",r),
2197                                    netabins,-1,1,nphibins,0,TMath::TwoPi());
2198         fhPhiBandTrack->SetXTitle("#eta");
2199         fhPhiBandTrack->SetYTitle("#phi");
2200         outputContainer->Add(fhPhiBandTrack) ;
2201
2202         fhConeSumPtEtaUESubTrack  = new TH2F("hConeSumPtEtaUESubTrack",
2203                                              Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2204                                              nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2205         fhConeSumPtEtaUESubTrack->SetYTitle("#Sigma #it{p}_{T}");
2206         fhConeSumPtEtaUESubTrack->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2207         outputContainer->Add(fhConeSumPtEtaUESubTrack) ;
2208
2209         fhConeSumPtPhiUESubTrack  = new TH2F("hConeSumPtPhiUESubTrack",
2210                                              Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2211                                              nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2212         fhConeSumPtPhiUESubTrack->SetYTitle("#Sigma #it{p}_{T}");
2213         fhConeSumPtPhiUESubTrack->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2214         outputContainer->Add(fhConeSumPtPhiUESubTrack) ;
2215
2216         fhConeSumPtEtaUESubTrackTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrackTrigEtaPhi",
2217                                                        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),
2218                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
2219         fhConeSumPtEtaUESubTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2220         fhConeSumPtEtaUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2221         fhConeSumPtEtaUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2222         outputContainer->Add(fhConeSumPtEtaUESubTrackTrigEtaPhi) ;
2223
2224         fhConeSumPtPhiUESubTrackTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrackTrigEtaPhi",
2225                                                        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),
2226                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
2227         fhConeSumPtPhiUESubTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2228         fhConeSumPtPhiUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2229         fhConeSumPtPhiUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2230         outputContainer->Add(fhConeSumPtPhiUESubTrackTrigEtaPhi) ;
2231
2232         fhFractionTrackOutConeEta  = new TH2F("hFractionTrackOutConeEta",
2233                                               Form("Fraction of the isolation cone #it{R} =  %2.2f, out of tracks #eta acceptance",r),
2234                                               nptbins,ptmin,ptmax,100,0,1);
2235         fhFractionTrackOutConeEta->SetYTitle("#it{fraction}");
2236         fhFractionTrackOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2237         outputContainer->Add(fhFractionTrackOutConeEta) ;
2238
2239         fhFractionTrackOutConeEtaTrigEtaPhi  = new TH2F("hFractionTrackOutConeEtaTrigEtaPhi",
2240                                                         Form("Fraction of the isolation cone #it{R} =  %2.2f, out of tracks #eta acceptance, in trigger #eta-#phi ",r),
2241                                                         netabins,etamin,etamax,nphibins,phimin,phimax);
2242         fhFractionTrackOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2243         fhFractionTrackOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2244         fhFractionTrackOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2245         outputContainer->Add(fhFractionTrackOutConeEtaTrigEtaPhi) ;
2246
2247         fhConeSumPtSubvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubvsConeSumPtTotPhiTrack",
2248                                                         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),
2249                                                         nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2250         fhConeSumPtSubvsConeSumPtTotPhiTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2251         fhConeSumPtSubvsConeSumPtTotPhiTrack->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2252         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiTrack);
2253
2254         fhConeSumPtSubNormvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiTrack",
2255                                                             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),
2256                                                             nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2257         fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2258         fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2259         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiTrack);
2260
2261         fhConeSumPtSubvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubvsConeSumPtTotEtaTrack",
2262                                                         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),
2263                                                         nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2264         fhConeSumPtSubvsConeSumPtTotEtaTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2265         fhConeSumPtSubvsConeSumPtTotEtaTrack->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2266         outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaTrack);
2267
2268         fhConeSumPtSubNormvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaTrack",
2269                                                             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),
2270                                                             nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2271         fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2272         fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2273         outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaTrack);
2274
2275
2276         // UE in perpendicular cone
2277         fhPerpConeSumPt  = new TH2F("hPerpConePtSum",
2278                                     Form("#Sigma #it{p}_{T} in isolation cone at #pm 45 degree phi from trigger particle, #it{R} =  %2.2f",r),
2279                                     nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2280         fhPerpConeSumPt->SetYTitle("#Sigma #it{p}_{T}");
2281         fhPerpConeSumPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2282         outputContainer->Add(fhPerpConeSumPt) ;
2283
2284         fhPtInPerpCone  = new TH2F("hPtInPerpCone",
2285                                    Form("#it{p}_{T} in isolation cone at #pm 45 degree phi from trigger particle, #it{R} =  %2.2f",r),
2286                                    nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2287         fhPtInPerpCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2288         fhPtInPerpCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2289         outputContainer->Add(fhPtInPerpCone) ;
2290
2291         fhEtaPhiTrack= new TH2F("hEtaPhiTrack",
2292                                 Form("#eta vs #phi of all Tracks"),
2293                                 netabins,-1,1,nphibins,0,TMath::TwoPi());
2294         fhEtaPhiTrack->SetXTitle("#eta");
2295         fhEtaPhiTrack->SetYTitle("#phi");
2296         outputContainer->Add(fhEtaPhiTrack) ;
2297
2298         fhEtaPhiInConeTrack= new TH2F("hEtaPhiInConeTrack",
2299                                       Form("#eta vs #phi of Tracks in cone for #it{R} =  %2.2f",r),
2300                                       netabins,-1,1,nphibins,0,TMath::TwoPi());
2301         fhEtaPhiInConeTrack->SetXTitle("#eta");
2302         fhEtaPhiInConeTrack->SetYTitle("#phi");
2303         outputContainer->Add(fhEtaPhiInConeTrack) ;
2304
2305         fhConeSumPtVSUETracksEtaBand  = new TH2F("hConeSumPtVSUETracksEtaBand",
2306                                                  Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in eta band for tracks (before normalization), R=%2.2f",r),
2307                                                  nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
2308         fhConeSumPtVSUETracksEtaBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2309         fhConeSumPtVSUETracksEtaBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2310         outputContainer->Add(fhConeSumPtVSUETracksEtaBand);
2311
2312         fhConeSumPtVSUETracksPhiBand  = new TH2F("hConeSumPtVSUETracksPhiBand",
2313                                                  Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in phi band for tracks (before normalization), R=%2.2f",r),
2314                                                  nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
2315         fhConeSumPtVSUETracksPhiBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2316         fhConeSumPtVSUETracksPhiBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2317         outputContainer->Add(fhConeSumPtVSUETracksPhiBand);
2318       }
2319     }
2320
2321     if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged )
2322     {
2323       fhConeSumPtClustervsTrack   = new TH2F("hConePtSumClustervsTrack",
2324                                              Form("Track vs Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2325                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2326       fhConeSumPtClustervsTrack->SetXTitle("#Sigma #it{p}_{T}^{cluster} (GeV/#it{c})");
2327       fhConeSumPtClustervsTrack->SetYTitle("#Sigma #it{p}_{T}^{track} (GeV/#it{c})");
2328       outputContainer->Add(fhConeSumPtClustervsTrack) ;
2329
2330       fhConeSumPtClusterTrackFrac   = new TH2F("hConePtSumClusterTrackFraction",
2331                                              Form("#Sigma #it{p}_{T}^{cluster}/#Sigma #it{p}_{T}^{track} in isolation cone for #it{R} =  %2.2f",r),
2332                                              nptbins,ptmin,ptmax,200,0,4);
2333       fhConeSumPtClusterTrackFrac->SetYTitle("#Sigma #it{p}^{cluster}_{T} /#Sigma #it{p}_{T}^{track}");
2334       fhConeSumPtClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
2335       outputContainer->Add(fhConeSumPtClusterTrackFrac) ;
2336
2337
2338       fhConePtLeadClustervsTrack   = new TH2F("hConePtLeadClustervsTrack",
2339                                              Form("Track vs Cluster lead #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2340                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2341       fhConePtLeadClustervsTrack->SetXTitle("#it{p}^{leading cluster}_{T} (GeV/#it{c})");
2342       fhConePtLeadClustervsTrack->SetYTitle("#it{p}^{leading track}_{T} (GeV/#it{c})");
2343       outputContainer->Add(fhConePtLeadClustervsTrack) ;
2344
2345       fhConePtLeadClusterTrackFrac   = new TH2F("hConePtLeadClusterTrackFraction",
2346                                                Form(" #it{p}^{leading cluster}_{T}/#it{p}^{leading track}_{T} in isolation cone for #it{R} =  %2.2f",r),
2347                                                nptbins,ptmin,ptmax,200,0,4);
2348       fhConePtLeadClusterTrackFrac->SetYTitle("#it{p}^{leading cluster}_{T}/ #it{p}^{leading track}_{T}");
2349       fhConePtLeadClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
2350       outputContainer->Add(fhConePtLeadClusterTrackFrac) ;
2351
2352
2353       if(fFillCellHistograms)
2354       {
2355         fhConeSumPtCellvsTrack   = new TH2F("hConePtSumCellvsTrack",
2356                                             Form("Track vs cell #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2357                                             nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2358         fhConeSumPtCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2359         fhConeSumPtCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2360         outputContainer->Add(fhConeSumPtCellvsTrack) ;
2361
2362         fhConeSumPtCellTrack = new TH2F("hConePtSumCellTrack",
2363                                         Form("Track and Cell #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2364                                         nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2365         fhConeSumPtCellTrack->SetYTitle("#Sigma #it{p}_{T}");
2366         fhConeSumPtCellTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2367         outputContainer->Add(fhConeSumPtCellTrack) ;
2368
2369         fhConeSumPtCellTrackTrigEtaPhi  = new TH2F("hConePtSumCellTrackTrigEtaPhi",
2370                                                    Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
2371                                                    netabins,etamin,etamax,nphibins,phimin,phimax);
2372         fhConeSumPtCellTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2373         fhConeSumPtCellTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2374         fhConeSumPtCellTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2375         outputContainer->Add(fhConeSumPtCellTrackTrigEtaPhi) ;
2376
2377       }
2378
2379       if(fFillUEBandSubtractHistograms)
2380       {
2381         fhConeSumPtEtaUESub  = new TH2F("hConeSumPtEtaUESub",
2382                                         Form("#Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2383                                         nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2384         fhConeSumPtEtaUESub->SetYTitle("#Sigma #it{p}_{T}");
2385         fhConeSumPtEtaUESub->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2386         outputContainer->Add(fhConeSumPtEtaUESub) ;
2387
2388         fhConeSumPtPhiUESub  = new TH2F("hConeSumPtPhiUESub",
2389                                         Form("#Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2390                                         nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2391         fhConeSumPtPhiUESub->SetYTitle("#Sigma #it{p}_{T}");
2392         fhConeSumPtPhiUESub->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2393         outputContainer->Add(fhConeSumPtPhiUESub) ;
2394
2395         fhConeSumPtEtaUESubTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrigEtaPhi",
2396                                                   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),
2397                                                   netabins,etamin,etamax,nphibins,phimin,phimax);
2398         fhConeSumPtEtaUESubTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2399         fhConeSumPtEtaUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2400         fhConeSumPtEtaUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2401         outputContainer->Add(fhConeSumPtEtaUESubTrigEtaPhi) ;
2402
2403         fhConeSumPtPhiUESubTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrigEtaPhi",
2404                                                   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),
2405                                                   netabins,etamin,etamax,nphibins,phimin,phimax);
2406         fhConeSumPtPhiUESubTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2407         fhConeSumPtPhiUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2408         fhConeSumPtPhiUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2409         outputContainer->Add(fhConeSumPtPhiUESubTrigEtaPhi) ;
2410
2411         fhConeSumPtEtaUESubClustervsTrack   = new TH2F("hConePtSumEtaUESubClustervsTrack",
2412                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} =  %2.2f",r),
2413                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2414         fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2415         fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2416         outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2417
2418         fhConeSumPtPhiUESubClustervsTrack   = new TH2F("hConePhiUESubPtSumClustervsTrack",
2419                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} =  %2.2f",r),
2420                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2421         fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2422         fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2423         outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2424
2425         fhEtaBandClustervsTrack   = new TH2F("hEtaBandClustervsTrack",
2426                                              Form("Track vs Cluster #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2427                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2428         fhEtaBandClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2429         fhEtaBandClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2430         outputContainer->Add(fhEtaBandClustervsTrack) ;
2431
2432         fhPhiBandClustervsTrack   = new TH2F("hPhiBandClustervsTrack",
2433                                              Form("Track vs Cluster #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2434                                              nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2435         fhPhiBandClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2436         fhPhiBandClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2437         outputContainer->Add(fhPhiBandClustervsTrack) ;
2438
2439         fhEtaBandNormClustervsTrack   = new TH2F("hEtaBandNormClustervsTrack",
2440                                                  Form("Track vs Cluster #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2441                                                  nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2442         fhEtaBandNormClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2443         fhEtaBandNormClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2444         outputContainer->Add(fhEtaBandNormClustervsTrack) ;
2445
2446         fhPhiBandNormClustervsTrack   = new TH2F("hPhiBandNormClustervsTrack",
2447                                                  Form("Track vs Cluster #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2448                                                  nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2449         fhPhiBandNormClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2450         fhPhiBandNormClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2451         outputContainer->Add(fhPhiBandNormClustervsTrack) ;
2452
2453         fhConeSumPtEtaUESubClustervsTrack   = new TH2F("hConePtSumEtaUESubClustervsTrack",
2454                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} =  %2.2f",r),
2455                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2456         fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2457         fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2458         outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2459
2460         fhConeSumPtPhiUESubClustervsTrack   = new TH2F("hConePhiUESubPtSumClustervsTrack",
2461                                                        Form("Track vs Cluster #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} =  %2.2f",r),
2462                                                        2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2463         fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2464         fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2465         outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2466
2467         if(fFillCellHistograms)
2468         {
2469
2470           fhConeSumPtEtaUESubCellvsTrack   = new TH2F("hConePtSumEtaUESubCellvsTrack",
2471                                                       Form("Track vs Cell #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} =  %2.2f",r),
2472                                                       2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2473           fhConeSumPtEtaUESubCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2474           fhConeSumPtEtaUESubCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2475           outputContainer->Add(fhConeSumPtEtaUESubCellvsTrack) ;
2476
2477           fhConeSumPtPhiUESubCellvsTrack   = new TH2F("hConePhiUESubPtSumCellvsTrack",
2478                                                       Form("Track vs Cell #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} =  %2.2f",r),
2479                                                       2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2480           fhConeSumPtPhiUESubCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2481           fhConeSumPtPhiUESubCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2482           outputContainer->Add(fhConeSumPtPhiUESubCellvsTrack) ;
2483
2484           fhEtaBandCellvsTrack   = new TH2F("hEtaBandCellvsTrack",
2485                                             Form("Track vs Cell #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2486                                             nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2487           fhEtaBandCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2488           fhEtaBandCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2489           outputContainer->Add(fhEtaBandCellvsTrack) ;
2490
2491           fhPhiBandCellvsTrack   = new TH2F("hPhiBandCellvsTrack",
2492                                             Form("Track vs Cell #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2493                                             nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2494           fhPhiBandCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2495           fhPhiBandCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2496           outputContainer->Add(fhPhiBandCellvsTrack) ;
2497
2498           fhEtaBandNormCellvsTrack   = new TH2F("hEtaBandNormCellvsTrack",
2499                                                 Form("Track vs Cell #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} =  %2.2f",r),
2500                                                 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2501           fhEtaBandNormCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2502           fhEtaBandNormCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2503           outputContainer->Add(fhEtaBandNormCellvsTrack) ;
2504
2505           fhPhiBandNormCellvsTrack   = new TH2F("hPhiBandNormCellvsTrack",
2506                                                 Form("Track vs Cell #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} =  %2.2f",r),
2507                                                 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2508           fhPhiBandNormCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2509           fhPhiBandNormCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2510           outputContainer->Add(fhPhiBandNormCellvsTrack) ;
2511
2512           fhConeSumPtEtaUESubTrackCell  = new TH2F("hConeSumPtEtaUESubTrackCell",
2513                                                    Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} =  %2.2f",r),
2514                                                    nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2515           fhConeSumPtEtaUESubTrackCell->SetYTitle("#Sigma #it{p}_{T}");
2516           fhConeSumPtEtaUESubTrackCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2517           outputContainer->Add(fhConeSumPtEtaUESubTrackCell) ;
2518
2519           fhConeSumPtPhiUESubTrackCell  = new TH2F("hConeSumPtPhiUESubTrackCell",
2520                                                    Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} =  %2.2f",r),
2521                                                    nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2522           fhConeSumPtPhiUESubTrackCell->SetYTitle("#Sigma #it{p}_{T}");
2523           fhConeSumPtPhiUESubTrackCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2524           outputContainer->Add(fhConeSumPtPhiUESubTrackCell) ;
2525
2526           fhConeSumPtEtaUESubTrackCellTrigEtaPhi  = new TH2F("hConeSumPtEtaUESubTrackCellTrigEtaPhi",
2527                                                              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),
2528                                                              netabins,etamin,etamax,nphibins,phimin,phimax);
2529           fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2530           fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2531           fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2532           outputContainer->Add(fhConeSumPtEtaUESubTrackCellTrigEtaPhi) ;
2533
2534           fhConeSumPtPhiUESubTrackCellTrigEtaPhi  = new TH2F("hConeSumPtPhiUESubTrackCellTrigEtaPhi",
2535                                                              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),
2536                                                              netabins,etamin,etamax,nphibins,phimin,phimax);
2537           fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2538           fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2539           fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2540           outputContainer->Add(fhConeSumPtPhiUESubTrackCellTrigEtaPhi) ;
2541         }
2542       }
2543     }
2544
2545     for(Int_t iso = 0; iso < 2; iso++)
2546     {
2547       if(fFillTMHisto)
2548       {
2549         fhTrackMatchedDEta[iso]  = new TH2F
2550         (Form("hTrackMatchedDEta%s",isoName[iso].Data()),
2551          Form("%s - d#eta of cluster-track vs cluster energy, %s",isoTitle[iso].Data(),parTitle.Data()),
2552          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2553         fhTrackMatchedDEta[iso]->SetYTitle("d#eta");
2554         fhTrackMatchedDEta[iso]->SetXTitle("E_{cluster} (GeV)");
2555
2556         fhTrackMatchedDPhi[iso]  = new TH2F
2557         (Form("hTrackMatchedDPhi%s",isoName[iso].Data()),
2558          Form("%s - d#phi of cluster-track vs cluster energy, %s",isoTitle[iso].Data(),parTitle.Data()),
2559          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2560         fhTrackMatchedDPhi[iso]->SetYTitle("d#phi (rad)");
2561         fhTrackMatchedDPhi[iso]->SetXTitle("E_{cluster} (GeV)");
2562
2563         fhTrackMatchedDEtaDPhi[iso]  = new TH2F
2564         (Form("hTrackMatchedDEtaDPhi%s",isoName[iso].Data()),
2565          Form("%s - d#eta vs d#phi of cluster-track, %s",isoTitle[iso].Data(),parTitle.Data()),
2566          nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2567         fhTrackMatchedDEtaDPhi[iso]->SetYTitle("d#phi (rad)");
2568         fhTrackMatchedDEtaDPhi[iso]->SetXTitle("d#eta");
2569
2570         outputContainer->Add(fhTrackMatchedDEta[iso]) ;
2571         outputContainer->Add(fhTrackMatchedDPhi[iso]) ;
2572         outputContainer->Add(fhTrackMatchedDEtaDPhi[iso]) ;
2573
2574         fhdEdx[iso]  = new TH2F
2575         (Form("hdEdx%s",isoName[iso].Data()),
2576          Form("%s - Matched track <d#it{E}/d#it{x}> vs cluster #it{E}, %s",isoTitle[iso].Data(),parTitle.Data()),
2577          nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2578         fhdEdx[iso]->SetXTitle("#it{E} (GeV)");
2579         fhdEdx[iso]->SetYTitle("<d#it{E}/d#it{x}>");
2580         outputContainer->Add(fhdEdx[iso]);
2581
2582         fhEOverP[iso]  = new TH2F
2583         (Form("hEOverP%s",isoName[iso].Data()),
2584          Form("%s - Matched track #it{E}/#it{p} vs cluster, %s",isoTitle[iso].Data(),parTitle.Data()),
2585          nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2586         fhEOverP[iso]->SetXTitle("#it{E} (GeV)");
2587         fhEOverP[iso]->SetYTitle("#it{E}/#it{p}");
2588         outputContainer->Add(fhEOverP[iso]);
2589
2590         if(IsDataMC())
2591         {
2592           fhTrackMatchedMCParticle[iso]  = new TH2F
2593           (Form("hTrackMatchedMCParticle%s",isoName[iso].Data()),
2594            Form("%s - Origin of particle vs cluster #it{E}, %s",isoTitle[iso].Data(),parTitle.Data()),
2595            nptbins,ptmin,ptmax,8,0,8);
2596           fhTrackMatchedMCParticle[iso]->SetXTitle("#it{E} (GeV)");
2597           //fhTrackMatchedMCParticle[iso]->SetYTitle("Particle type");
2598
2599           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(1 ,"Photon");
2600           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(2 ,"Electron");
2601           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
2602           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(4 ,"Rest");
2603           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
2604           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
2605           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
2606           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
2607
2608           outputContainer->Add(fhTrackMatchedMCParticle[iso]);
2609         }
2610       }
2611
2612       if(fFillSSHisto)
2613       {
2614         fhELambda0[iso]  = new TH2F
2615         (Form("hELambda0%s",isoName[iso].Data()),
2616          Form("%s cluster : #it{E} vs #lambda_{0}, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2617         fhELambda0[iso]->SetYTitle("#lambda_{0}^{2}");
2618         fhELambda0[iso]->SetXTitle("#it{E} (GeV)");
2619         outputContainer->Add(fhELambda0[iso]) ;
2620
2621         fhELambda1[iso]  = new TH2F
2622         (Form("hELambda1%s",isoName[iso].Data()),
2623          Form("%s cluster: #it{E} vs #lambda_{1}, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2624         fhELambda1[iso]->SetYTitle("#lambda_{1}^{2}");
2625         fhELambda1[iso]->SetXTitle("#it{E} (GeV)");
2626         outputContainer->Add(fhELambda1[iso]) ;
2627
2628         fhPtLambda0[iso]  = new TH2F
2629         (Form("hPtLambda0%s",isoName[iso].Data()),
2630          Form("%s cluster : #it{p}_{T} vs #lambda_{0}, %s",isoTitle[iso].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2631         fhPtLambda0[iso]->SetYTitle("#lambda_{0}^{2}");
2632         fhPtLambda0[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2633         outputContainer->Add(fhPtLambda0[iso]) ;
2634
2635         if(IsDataMC())
2636         {
2637           for(Int_t imc = 0; imc < 9; imc++)
2638           {
2639             fhPtLambda0MC[imc][iso]  = new TH2F(Form("hPtLambda0%s_MC%s",isoName[iso].Data(),mcPartName[imc].Data()),
2640                                                 Form("%s cluster : #it{p}_{T} vs #lambda_{0}: %s %s",isoTitle[iso].Data(),mcPartType[imc].Data(),parTitle.Data()),
2641                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2642             fhPtLambda0MC[imc][iso]->SetYTitle("#lambda_{0}^{2}");
2643             fhPtLambda0MC[imc][iso]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
2644             outputContainer->Add( fhPtLambda0MC[imc][iso]) ;
2645           }
2646         }
2647
2648         if(fIsoDetector=="EMCAL" &&  GetFirstSMCoveredByTRD() >= 0)
2649         {
2650           fhPtLambda0TRD[iso]  = new TH2F
2651           (Form("hPtLambda0TRD%s",isoName[iso].Data()),
2652            Form("%s cluster: #it{p}_{T} vs #lambda_{0}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2653           fhPtLambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
2654           fhPtLambda0TRD[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2655           outputContainer->Add(fhPtLambda0TRD[iso]) ;
2656
2657           fhELambda0TRD[iso]  = new TH2F
2658           (Form("hELambda0TRD%s",isoName[iso].Data()),
2659            Form("%s cluster: #it{E} vs #lambda_{0}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2660           fhELambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
2661           fhELambda0TRD[iso]->SetXTitle("#it{E} (GeV)");
2662           outputContainer->Add(fhELambda0TRD[iso]) ;
2663
2664           fhELambda1TRD[iso]  = new TH2F
2665           (Form("hELambda1TRD%s",isoName[iso].Data()),
2666            Form("%s cluster: #it{E} vs #lambda_{1}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2667           fhELambda1TRD[iso]->SetYTitle("#lambda_{1}^{2}");
2668           fhELambda1TRD[iso]->SetXTitle("#it{E} (GeV)");
2669           outputContainer->Add(fhELambda1TRD[iso]) ;
2670         }
2671
2672         if(fFillNLMHistograms)
2673         {
2674           fhNLocMax[iso] = new TH2F
2675           (Form("hNLocMax%s",isoName[iso].Data()),
2676            Form("%s - Number of local maxima in cluster, %s",isoTitle[iso].Data(),parTitle.Data()),
2677            nptbins,ptmin,ptmax,10,0,10);
2678           fhNLocMax[iso]->SetYTitle("#it{NLM}");
2679           fhNLocMax[iso]->SetXTitle("#it{E} (GeV)");
2680           outputContainer->Add(fhNLocMax[iso]) ;
2681
2682           fhELambda0LocMax1[iso]  = new TH2F
2683           (Form("hELambda0LocMax1%s",isoName[iso].Data()),
2684            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);
2685           fhELambda0LocMax1[iso]->SetYTitle("#lambda_{0}^{2}");
2686           fhELambda0LocMax1[iso]->SetXTitle("#it{E} (GeV)");
2687           outputContainer->Add(fhELambda0LocMax1[iso]) ;
2688
2689           fhELambda1LocMax1[iso]  = new TH2F
2690           (Form("hELambda1LocMax1%s",isoName[iso].Data()),
2691            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);
2692           fhELambda1LocMax1[iso]->SetYTitle("#lambda_{1}^{2}");
2693           fhELambda1LocMax1[iso]->SetXTitle("#it{E} (GeV)");
2694           outputContainer->Add(fhELambda1LocMax1[iso]) ;
2695
2696           fhELambda0LocMax2[iso]  = new TH2F
2697           (Form("hELambda0LocMax2%s",isoName[iso].Data()),
2698            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);
2699           fhELambda0LocMax2[iso]->SetYTitle("#lambda_{0}^{2}");
2700           fhELambda0LocMax2[iso]->SetXTitle("#it{E} (GeV)");
2701           outputContainer->Add(fhELambda0LocMax2[iso]) ;
2702
2703           fhELambda1LocMax2[iso]  = new TH2F
2704           (Form("hELambda1LocMax2%s",isoName[iso].Data()),
2705            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);
2706           fhELambda1LocMax2[iso]->SetYTitle("#lambda_{1}^{2}");
2707           fhELambda1LocMax2[iso]->SetXTitle("#it{E} (GeV)");
2708           outputContainer->Add(fhELambda1LocMax2[iso]) ;
2709
2710           fhELambda0LocMaxN[iso]  = new TH2F
2711           ( Form("hELambda0LocMaxN%s",isoName[iso].Data()),
2712            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);
2713           fhELambda0LocMaxN[iso]->SetYTitle("#lambda_{0}^{2}");
2714           fhELambda0LocMaxN[iso]->SetXTitle("#it{E} (GeV)");
2715           outputContainer->Add(fhELambda0LocMaxN[iso]) ;
2716
2717           fhELambda1LocMaxN[iso]  = new TH2F
2718           (Form("hELambda1LocMaxN%s",isoName[iso].Data()),
2719            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);
2720           fhELambda1LocMaxN[iso]->SetYTitle("#lambda_{1}^{2}");
2721           fhELambda1LocMaxN[iso]->SetXTitle("#it{E} (GeV)");
2722           outputContainer->Add(fhELambda1LocMaxN[iso]) ;
2723         } // NLM
2724       } // SS histo
2725     } // control histograms for isolated and non isolated objects
2726
2727
2728     if(fFillPileUpHistograms)
2729     {
2730       fhPtTrackInConeOtherBC  = new TH2F("hPtTrackInConeOtherBC",
2731                                          Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC!=0",r),
2732                                          nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2733       fhPtTrackInConeOtherBC->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2734       fhPtTrackInConeOtherBC->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2735       outputContainer->Add(fhPtTrackInConeOtherBC) ;
2736
2737       fhPtTrackInConeOtherBCPileUpSPD  = new TH2F("hPtTrackInConeOtherBCPileUpSPD",
2738                                                   Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC!=0, pile-up from SPD",r),
2739                                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2740       fhPtTrackInConeOtherBCPileUpSPD->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2741       fhPtTrackInConeOtherBCPileUpSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2742       outputContainer->Add(fhPtTrackInConeOtherBCPileUpSPD) ;
2743
2744       fhPtTrackInConeBC0  = new TH2F("hPtTrackInConeBC0",
2745                                      Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC==0",r),
2746                                      nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2747       fhPtTrackInConeBC0->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2748       fhPtTrackInConeBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2749       outputContainer->Add(fhPtTrackInConeBC0) ;
2750
2751       fhPtTrackInConeVtxBC0  = new TH2F("hPtTrackInConeVtxBC0",
2752                                         Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC==0",r),
2753                                         nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2754       fhPtTrackInConeVtxBC0->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2755       fhPtTrackInConeVtxBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2756       outputContainer->Add(fhPtTrackInConeVtxBC0) ;
2757
2758
2759       fhPtTrackInConeBC0PileUpSPD  = new TH2F("hPtTrackInConeBC0PileUpSPD",
2760                                               Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f, TOF from BC==0, pile-up from SPD",r),
2761                                               nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2762       fhPtTrackInConeBC0PileUpSPD->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2763       fhPtTrackInConeBC0PileUpSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2764       outputContainer->Add(fhPtTrackInConeBC0PileUpSPD) ;
2765
2766
2767       for (Int_t i = 0; i < 7 ; i++)
2768       {
2769         fhPtInConePileUp[i]  = new TH2F(Form("hPtInConePileUp%s",pileUpName[i].Data()),
2770                                         Form("#it{p}_{T} in isolation cone for #it{R} =  %2.2f, from pile-up (%s)",r,pileUpName[i].Data()),
2771                                         nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2772         fhPtInConePileUp[i]->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2773         fhPtInConePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2774         outputContainer->Add(fhPtInConePileUp[i]) ;
2775       }
2776     }
2777
2778     if(IsDataMC())
2779     {
2780       // For histograms in arrays, index in the array, corresponding to any particle origin
2781
2782       for(Int_t i = 0; i < 6; i++)
2783       {
2784         fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
2785                                  Form("primary photon  %s : #it{E}, %s",pptype[i].Data(),parTitle.Data()),
2786                                  nptbins,ptmin,ptmax);
2787         fhEPrimMC[i]->SetXTitle("#it{E} (GeV)");
2788         outputContainer->Add(fhEPrimMC[i]) ;
2789
2790         fhPtPrimMCiso[i]  = new TH1F(Form("hPtPrim_MCiso%s",ppname[i].Data()),
2791                                      Form("primary isolated photon %s : #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2792                                      nptbins,ptmin,ptmax);
2793         fhPtPrimMCiso[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2794         outputContainer->Add(fhPtPrimMCiso[i]) ;
2795
2796         fhEtaPrimMC[i]  = new TH2F(Form("hEtaPrim_MC%s",ppname[i].Data()),
2797                                    Form("primary photon %s : #eta vs #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2798                                    nptbins,ptmin,ptmax,200,-2,2);
2799         fhEtaPrimMC[i]->SetYTitle("#eta");
2800         fhEtaPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2801         outputContainer->Add(fhEtaPrimMC[i]) ;
2802
2803         fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
2804                                    Form("primary photon %s : #phi vs #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2805                                    nptbins,ptmin,ptmax,200,0.,TMath::TwoPi());
2806         fhPhiPrimMC[i]->SetYTitle("#phi");
2807         fhPhiPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2808         outputContainer->Add(fhPhiPrimMC[i]) ;
2809       }
2810
2811     }//Histos with MC
2812
2813   }
2814
2815   if(fMakeSeveralIC)
2816   {
2817     const Int_t buffersize = 255;
2818     char name[buffersize];
2819     char title[buffersize];
2820     for(Int_t icone = 0; icone<fNCones; icone++)
2821     {
2822       // sum pt in cone vs. pt leading
2823       snprintf(name, buffersize,"hSumPtLeadingPt_Cone_%d",icone);
2824       snprintf(title, buffersize,"#Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",fConeSizes[icone]);
2825       fhSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2826       fhSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}#it{p}_{T} (GeV/#it{c})");//#Sigma #it{p}_{T}
2827       fhSumPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2828       outputContainer->Add(fhSumPtLeadingPt[icone]) ;
2829
2830       // pt in cone vs. pt leading
2831       snprintf(name, buffersize,"hPtLeadingPt_Cone_%d",icone);
2832       snprintf(title, buffersize,"#it{p}_{T} in isolation cone for #it{R} =  %2.2f",fConeSizes[icone]);
2833       fhPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2834       fhPtLeadingPt[icone] ->SetYTitle("#it{p}_{T}^{cone} (GeV/#it{c})");
2835       fhPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2836       outputContainer->Add(fhPtLeadingPt[icone]) ;
2837
2838       // sum pt in cone vs. pt leading in the forward region (for background subtraction studies)
2839       snprintf(name, buffersize,"hPerpSumPtLeadingPt_Cone_%d",icone);
2840       snprintf(title, buffersize,"#Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",fConeSizes[icone]);
2841       fhPerpSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2842       fhPerpSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}#it{p}_{T} (GeV/#it{c})");//#Sigma #it{p}_{T}
2843       fhPerpSumPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2844       outputContainer->Add(fhPerpSumPtLeadingPt[icone]) ;
2845
2846       // pt in cone vs. pt leading in the forward region (for background subtraction studies)
2847       snprintf(name, buffersize,"hPerpPtLeadingPt_Cone_%d",icone);
2848       snprintf(title, buffersize,"#it{p}_{T} in isolation cone for #it{R} =  %2.2f",fConeSizes[icone]);
2849       fhPerpPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2850       fhPerpPtLeadingPt[icone] ->SetYTitle("#it{p}_{T}^{cone} (GeV/#it{c})");
2851       fhPerpPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2852       outputContainer->Add(fhPerpPtLeadingPt[icone]) ;
2853
2854       if(IsDataMC())
2855       {
2856         for(Int_t imc = 0; imc < 9; imc++)
2857         {
2858           snprintf(name , buffersize,"hSumPtLeadingPt_MC%s_Cone_%d",mcPartName[imc].Data(),icone);
2859           snprintf(title, buffersize,"Candidate %s #it{p}_{T} vs cone #Sigma #it{p}_{T} for #it{R}=%2.2f",mcPartType[imc].Data(),fConeSizes[icone]);
2860           fhSumPtLeadingPtMC[imc][icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2861           fhSumPtLeadingPtMC[imc][icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2862           fhSumPtLeadingPtMC[imc][icone]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
2863           outputContainer->Add(fhSumPtLeadingPtMC[imc][icone]) ;
2864         }
2865       }//Histos with MC
2866
2867       for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++)
2868       {
2869         snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
2870         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]);
2871         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2872         fhPtThresIsolated[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2873         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
2874
2875         snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
2876         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]);
2877         fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2878         fhPtFracIsolated[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2879         outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
2880
2881         snprintf(name, buffersize,"hSumPt_Cone_%d_Pt%d",icone,ipt);
2882         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]);
2883         fhSumPtIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2884         // fhSumPtIsolated[icone][ipt]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2885         fhSumPtIsolated[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2886         outputContainer->Add(fhSumPtIsolated[icone][ipt]) ;
2887
2888         snprintf(name, buffersize,"hPtSumDensity_Cone_%d_Pt%d",icone,ipt);
2889         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]);
2890         fhPtSumDensityIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2891         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2892         fhPtSumDensityIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2893         outputContainer->Add(fhPtSumDensityIso[icone][ipt]) ;
2894
2895         snprintf(name, buffersize,"hPtFracPtSum_Cone_%d_Pt%d",icone,ipt);
2896         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]);
2897         fhPtFracPtSumIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2898         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2899         fhPtFracPtSumIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2900         outputContainer->Add(fhPtFracPtSumIso[icone][ipt]) ;
2901
2902         // eta:phi
2903         snprintf(name, buffersize,"hEtaPhiPtThres_Cone_%d_Pt%d",icone,ipt);
2904         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]);
2905         fhEtaPhiPtThresIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2906         fhEtaPhiPtThresIso[icone][ipt]->SetXTitle("#eta");
2907         fhEtaPhiPtThresIso[icone][ipt]->SetYTitle("#phi");
2908         outputContainer->Add(fhEtaPhiPtThresIso[icone][ipt]) ;
2909
2910         snprintf(name, buffersize,"hEtaPhiPtFrac_Cone_%d_Pt%d",icone,ipt);
2911         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]);
2912         fhEtaPhiPtFracIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2913         fhEtaPhiPtFracIso[icone][ipt]->SetXTitle("#eta");
2914         fhEtaPhiPtFracIso[icone][ipt]->SetYTitle("#phi");
2915         outputContainer->Add(fhEtaPhiPtFracIso[icone][ipt]) ;
2916
2917         snprintf(name, buffersize,"hEtaPhiPtSum_Cone_%d_Pt%d",icone,ipt);
2918         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]);
2919         fhEtaPhiPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2920         fhEtaPhiPtSumIso[icone][ipt]->SetXTitle("#eta");
2921         fhEtaPhiPtSumIso[icone][ipt]->SetYTitle("#phi");
2922         outputContainer->Add(fhEtaPhiPtSumIso[icone][ipt]) ;
2923
2924         snprintf(name, buffersize,"hEtaPhiSumDensity_Cone_%d_Pt%d",icone,ipt);
2925         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]);
2926         fhEtaPhiSumDensityIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2927         fhEtaPhiSumDensityIso[icone][ipt]->SetXTitle("#eta");
2928         fhEtaPhiSumDensityIso[icone][ipt]->SetYTitle("#phi");
2929         outputContainer->Add(fhEtaPhiSumDensityIso[icone][ipt]) ;
2930
2931         snprintf(name, buffersize,"hEtaPhiFracPtSum_Cone_%d_Pt%d",icone,ipt);
2932         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]);
2933         fhEtaPhiFracPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2934         fhEtaPhiFracPtSumIso[icone][ipt]->SetXTitle("#eta");
2935         fhEtaPhiFracPtSumIso[icone][ipt]->SetYTitle("#phi");
2936         outputContainer->Add(fhEtaPhiFracPtSumIso[icone][ipt]) ;
2937
2938         if(fFillTaggedDecayHistograms)
2939         {
2940           // pt decays isolated
2941           snprintf(name, buffersize,"hPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
2942           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]);
2943           fhPtPtThresDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2944           fhPtPtThresDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2945           outputContainer->Add(fhPtPtThresDecayIso[icone][ipt]) ;
2946
2947           snprintf(name, buffersize,"hPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
2948           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]);
2949           fhPtPtFracDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
2950           fhPtPtFracDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2951           outputContainer->Add(fhPtPtFracDecayIso[icone][ipt]) ;
2952
2953           snprintf(name, buffersize,"hPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2954           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]);
2955           fhPtPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2956           //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2957           fhPtPtSumDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2958           outputContainer->Add(fhPtPtSumDecayIso[icone][ipt]) ;
2959
2960           snprintf(name, buffersize,"hPtSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
2961           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]);
2962           fhPtSumDensityDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2963           //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2964           fhPtSumDensityDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2965           outputContainer->Add(fhPtSumDensityDecayIso[icone][ipt]) ;
2966
2967           snprintf(name, buffersize,"hPtFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2968           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]);
2969           fhPtFracPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
2970           //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2971           fhPtFracPtSumDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2972           outputContainer->Add(fhPtFracPtSumDecayIso[icone][ipt]) ;
2973
2974           // eta:phi decays
2975           snprintf(name, buffersize,"hEtaPhiPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
2976           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]);
2977           fhEtaPhiPtThresDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2978           fhEtaPhiPtThresDecayIso[icone][ipt]->SetXTitle("#eta");
2979           fhEtaPhiPtThresDecayIso[icone][ipt]->SetYTitle("#phi");
2980           outputContainer->Add(fhEtaPhiPtThresDecayIso[icone][ipt]) ;
2981
2982           snprintf(name, buffersize,"hEtaPhiPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
2983           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]);
2984           fhEtaPhiPtFracDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2985           fhEtaPhiPtFracDecayIso[icone][ipt]->SetXTitle("#eta");
2986           fhEtaPhiPtFracDecayIso[icone][ipt]->SetYTitle("#phi");
2987           outputContainer->Add(fhEtaPhiPtFracDecayIso[icone][ipt]) ;
2988
2989
2990           snprintf(name, buffersize,"hEtaPhiPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
2991           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]);
2992           fhEtaPhiPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
2993           fhEtaPhiPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
2994           fhEtaPhiPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
2995           outputContainer->Add(fhEtaPhiPtSumDecayIso[icone][ipt]) ;
2996
2997           snprintf(name, buffersize,"hEtaPhiSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
2998           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]);
2999           fhEtaPhiSumDensityDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3000           fhEtaPhiSumDensityDecayIso[icone][ipt]->SetXTitle("#eta");
3001           fhEtaPhiSumDensityDecayIso[icone][ipt]->SetYTitle("#phi");
3002           outputContainer->Add(fhEtaPhiSumDensityDecayIso[icone][ipt]) ;
3003
3004           snprintf(name, buffersize,"hEtaPhiFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
3005           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]);
3006           fhEtaPhiFracPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3007           fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
3008           fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
3009           outputContainer->Add(fhEtaPhiFracPtSumDecayIso[icone][ipt]) ;
3010
3011         }
3012
3013         if(IsDataMC())
3014         {
3015           for(Int_t imc = 0; imc < 9; imc++)
3016           {
3017             snprintf(name , buffersize,"hPtThreshMC%s_Cone_%d_Pt%d",mcPartName[imc].Data(),icone,ipt);
3018             snprintf(title, buffersize,"Isolated %s #it{p}_{T} for #it{R}=%2.2f and #it{p}_{T}^{th}=%2.2f",
3019                      mcPartType[imc].Data(),fConeSizes[icone], fPtThresholds[ipt]);
3020             fhPtThresIsolatedMC[imc][icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3021             fhPtThresIsolatedMC[imc][icone][ipt]->SetYTitle("#it{counts}");
3022             fhPtThresIsolatedMC[imc][icone][ipt]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3023             outputContainer->Add(fhPtThresIsolatedMC[imc][icone][ipt]) ;
3024
3025
3026             snprintf(name , buffersize,"hPtFracMC%s_Cone_%d_Pt%d",mcPartName[imc].Data(),icone,ipt);
3027             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",
3028                      mcPartType[imc].Data(),fConeSizes[icone], fPtFractions[ipt]);
3029             fhPtFracIsolatedMC[imc][icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3030             fhPtFracIsolatedMC[imc][icone][ipt]->SetYTitle("#it{counts}");
3031             fhPtFracIsolatedMC[imc][icone][ipt]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3032             outputContainer->Add(fhPtFracIsolatedMC[imc][icone][ipt]) ;
3033
3034             snprintf(name , buffersize,"hSumPtMC%s_Cone_%d_Pt%d",mcPartName[imc].Data(),icone,ipt);
3035             snprintf(title, buffersize,"Isolated %s #it{p}_{T} for #it{R}=%2.2f and #Sigma #it{p}_{T}^{in cone}=%2.2f",
3036                      mcPartType[imc].Data(),fConeSizes[icone], fSumPtThresholds[ipt]);
3037             fhSumPtIsolatedMC[imc][icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3038             fhSumPtIsolatedMC[imc][icone][ipt]->SetYTitle("#it{counts}");
3039             fhSumPtIsolatedMC[imc][icone][ipt]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3040             outputContainer->Add(fhSumPtIsolatedMC[imc][icone][ipt]) ;
3041           }
3042         }//Histos with MC
3043       }//icone loop
3044     }//ipt loop
3045   }
3046
3047   if(fFillPileUpHistograms)
3048   {
3049     for (Int_t i = 0; i < 7 ; i++)
3050     {
3051       fhEIsoPileUp[i]   = new TH1F(Form("hEPileUp%s",pileUpName[i].Data()),
3052                                    Form("Number of isolated particles vs E, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3053                                    nptbins,ptmin,ptmax);
3054       fhEIsoPileUp[i]->SetYTitle("d#it{N} / d#it{E}");
3055       fhEIsoPileUp[i]->SetXTitle("#it{E} (GeV)");
3056       outputContainer->Add(fhEIsoPileUp[i]) ;
3057
3058       fhPtIsoPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
3059                                    Form("Number of isolated particles vs #it{p}_{T}, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3060                                    nptbins,ptmin,ptmax);
3061       fhPtIsoPileUp[i]->SetYTitle("d#it{N} / #it{p}_{T}");
3062       fhPtIsoPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3063       outputContainer->Add(fhPtIsoPileUp[i]) ;
3064
3065       fhENoIsoPileUp[i]   = new TH1F(Form("hENoIsoPileUp%s",pileUpName[i].Data()),
3066                                      Form("Number of not isolated particles vs E, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3067                                      nptbins,ptmin,ptmax);
3068       fhENoIsoPileUp[i]->SetYTitle("d#it{N} / dE");
3069       fhENoIsoPileUp[i]->SetXTitle("#it{E} (GeV)");
3070       outputContainer->Add(fhENoIsoPileUp[i]) ;
3071
3072       fhPtNoIsoPileUp[i]  = new TH1F(Form("hPtNoIsoPileUp%s",pileUpName[i].Data()),
3073                                      Form("Number of not isolated particles vs #it{p}_{T}, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3074                                      nptbins,ptmin,ptmax);
3075       fhPtNoIsoPileUp[i]->SetYTitle("d#it{N} / #it{p}_{T}");
3076       fhPtNoIsoPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3077       outputContainer->Add(fhPtNoIsoPileUp[i]) ;
3078     }
3079
3080     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3081     fhTimeENoCut->SetXTitle("#it{E} (GeV)");
3082     fhTimeENoCut->SetYTitle("#it{time} (ns)");
3083     outputContainer->Add(fhTimeENoCut);
3084
3085     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3086     fhTimeESPD->SetXTitle("#it{E} (GeV)");
3087     fhTimeESPD->SetYTitle("#it{time} (ns)");
3088     outputContainer->Add(fhTimeESPD);
3089
3090     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3091     fhTimeESPDMulti->SetXTitle("#it{E} (GeV)");
3092     fhTimeESPDMulti->SetYTitle("#it{time} (ns)");
3093     outputContainer->Add(fhTimeESPDMulti);
3094
3095     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3096     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
3097     fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
3098     outputContainer->Add(fhTimeNPileUpVertSPD);
3099
3100     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
3101     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
3102     fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
3103     outputContainer->Add(fhTimeNPileUpVertTrack);
3104
3105     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3106     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
3107     fhTimeNPileUpVertContributors->SetXTitle("#it{time} (ns)");
3108     outputContainer->Add(fhTimeNPileUpVertContributors);
3109
3110     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);
3111     fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{z} (cm) ");
3112     fhTimePileUpMainVertexZDistance->SetXTitle("#it{time} (ns)");
3113     outputContainer->Add(fhTimePileUpMainVertexZDistance);
3114
3115     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
3116     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{z} (cm) ");
3117     fhTimePileUpMainVertexZDiamond->SetXTitle("#it{time} (ns)");
3118     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
3119   }
3120
3121   return outputContainer ;
3122
3123 }
3124
3125 //____________________________________________________
3126 Int_t AliAnaParticleIsolation::GetMCIndex(Int_t mcTag)
3127 {
3128   // Histogram index depending on origin of candidate
3129
3130   if(!IsDataMC()) return -1;
3131
3132   if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))
3133   {
3134     return kmcPrompt;
3135   }
3136   else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation) ||
3137           GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCISR))
3138   {
3139     return kmcFragment;
3140   }
3141   else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))
3142   {
3143     return kmcPi0;
3144   }
3145   else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))
3146   {
3147     return kmcPi0Decay;
3148   }
3149   else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))
3150   {
3151     return kmcEtaDecay;
3152   }
3153   else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))
3154   {
3155     return kmcOtherDecay;
3156   }
3157   else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
3158   {
3159     return kmcElectron;
3160   }
3161   else // anything else
3162   {
3163     // careful can contain also other decays, to be checked.
3164     return kmcHadron;
3165   }
3166 }
3167
3168 //__________________________________
3169 void AliAnaParticleIsolation::Init()
3170 {
3171   // Do some checks and init stuff
3172
3173   // In case of several cone and thresholds analysis, open the cuts for the filling of the
3174   // track and cluster reference arrays in cone when done in the MakeAnalysisFillAOD().
3175   // The different cones, thresholds are tested for this list of tracks, clusters.
3176   if(fMakeSeveralIC)
3177   {
3178     printf("AliAnaParticleIsolation::Init() - Open default isolation cuts for multiple Isolation analysis\n");
3179     GetIsolationCut()->SetPtThreshold(100);
3180     GetIsolationCut()->SetPtFraction(100);
3181     GetIsolationCut()->SetConeSize(1);
3182   }
3183
3184   if(!GetReader()->IsCTSSwitchedOn() && GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyNeutral)
3185     AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
3186
3187 }
3188
3189 //____________________________________________
3190 void AliAnaParticleIsolation::InitParameters()
3191 {
3192
3193   //Initialize the parameters of the analysis.
3194   SetInputAODName("PWG4Particle");
3195   SetAODObjArrayName("IsolationCone");
3196   AddToHistogramsName("AnaIsolation_");
3197
3198   fCalorimeter = "EMCAL" ;
3199   fIsoDetector = "EMCAL" ;
3200
3201   fReMakeIC = kFALSE ;
3202   fMakeSeveralIC = kFALSE ;
3203
3204   fLeadingOnly = kTRUE;
3205   fCheckLeadingWithNeutralClusters = kTRUE;
3206
3207   fNDecayBits = 1;
3208   fDecayBits[0] = AliNeutralMesonSelection::kPi0;
3209   fDecayBits[1] = AliNeutralMesonSelection::kEta;
3210   fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
3211   fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
3212
3213   fNBkgBin = 11;
3214   fBkgBinLimit[ 0] = 00.0; fBkgBinLimit[ 1] = 00.2; fBkgBinLimit[ 2] = 00.3; fBkgBinLimit[ 3] = 00.4; fBkgBinLimit[ 4] = 00.5;
3215   fBkgBinLimit[ 5] = 01.0; fBkgBinLimit[ 6] = 01.5; fBkgBinLimit[ 7] = 02.0; fBkgBinLimit[ 8] = 03.0; fBkgBinLimit[ 9] = 05.0;
3216   fBkgBinLimit[10] = 10.0; fBkgBinLimit[11] = 100.;
3217   for(Int_t ibin = 12; ibin < 20; ibin++) fBkgBinLimit[ibin] = 00.0;
3218
3219   //----------- Several IC-----------------
3220   fNCones             = 5 ;
3221   fNPtThresFrac       = 5 ;
3222   fConeSizes      [0] = 0.1;     fConeSizes      [1] = 0.2;   fConeSizes      [2] = 0.3; fConeSizes      [3] = 0.4;  fConeSizes      [4] = 0.5;
3223   fPtThresholds   [0] = 1.;      fPtThresholds   [1] = 2.;    fPtThresholds   [2] = 3.;  fPtThresholds   [3] = 4.;   fPtThresholds   [4] = 5.;
3224   fPtFractions    [0] = 0.05;    fPtFractions    [1] = 0.075; fPtFractions    [2] = 0.1; fPtFractions    [3] = 1.25; fPtFractions    [4] = 1.5;
3225   fSumPtThresholds[0] = 1.;      fSumPtThresholds[1] = 2.;    fSumPtThresholds[2] = 3.;  fSumPtThresholds[3] = 4.;   fSumPtThresholds[4] = 5.;
3226
3227 }
3228
3229 //_________________________________________________________________________________________
3230 Bool_t AliAnaParticleIsolation::IsTriggerTheNearSideEventLeadingParticle(Int_t & idLeading)
3231 {
3232   // Check if the what of the selected isolation candidates is leading particle in the same hemisphere
3233   // comparing with all the candidates, all the tracks or all the clusters.
3234
3235   Double_t ptTrig      = GetMinPt();
3236   Double_t phiTrig     = 0 ;
3237   Int_t index          =-1 ;
3238   AliAODPWG4ParticleCorrelation* pLeading = 0;
3239
3240   // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
3241
3242   for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3243   {
3244     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3245     particle->SetLeadingParticle(kFALSE); // set it later
3246
3247     // Vertex cut in case of mixing
3248     if(GetMixedEvent())
3249     {
3250       Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3251       if(check ==  0) continue;
3252       if(check == -1) return kFALSE; // not sure if it is correct.
3253     }
3254
3255     //check if it is low pt trigger particle
3256     if((particle->Pt() < GetIsolationCut()->GetPtThreshold() ||
3257         particle->Pt() < GetIsolationCut()->GetSumPtThreshold()) &&
3258        !fMakeSeveralIC)
3259     {
3260       continue ; //trigger should not come from underlying event
3261     }
3262
3263     // find the leading particles with highest momentum
3264     if (particle->Pt() > ptTrig)
3265     {
3266       ptTrig   = particle->Pt() ;
3267       phiTrig  = particle->Phi();
3268       index    = iaod     ;
3269       pLeading = particle ;
3270     }
3271   }// finish search of leading trigger particle on the AOD branch.
3272
3273   if(index < 0) return kFALSE;
3274
3275   idLeading = index;
3276
3277   //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3278
3279   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3280
3281   // Compare if it is the leading of all tracks
3282
3283   TVector3 p3;
3284   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3285   {
3286     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3287
3288     if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
3289        track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3)   ) continue ;
3290
3291     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
3292     p3.SetXYZ(mom[0],mom[1],mom[2]);
3293     Float_t pt   = p3.Pt();
3294     Float_t phi  = p3.Phi() ;
3295     if(phi < 0) phi+=TMath::TwoPi();
3296
3297     //skip this event if near side associated particle pt larger than trigger
3298
3299     Float_t deltaPhi = phiTrig-phi;
3300     //
3301     // Calculate deltaPhi shift so that for the particles on the opposite side
3302     // it is defined between 90 and 270 degrees
3303     // Shift [-360,-90]  to [0, 270]
3304     // and [270,360] to [-90,0]
3305     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3306     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3307
3308     if(pt > ptTrig && deltaPhi < TMath::PiOver2())  return kFALSE;
3309
3310   }// track loop
3311
3312   // Compare if it is leading of all calorimeter clusters
3313
3314   if(fCheckLeadingWithNeutralClusters)
3315   {
3316     // Select the calorimeter cluster list
3317     TObjArray * nePl = 0x0;
3318     if      (pLeading->GetDetector() == "PHOS" )
3319       nePl = GetPHOSClusters();
3320     else
3321       nePl = GetEMCALClusters();
3322
3323     if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3324
3325     TLorentzVector lv;
3326     for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3327     {
3328       AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3329
3330       if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3331
3332       cluster->GetMomentum(lv,GetVertex(0));
3333
3334       Float_t pt   = lv.Pt();
3335       Float_t phi  = lv.Phi() ;
3336       if(phi < 0) phi+=TMath::TwoPi();
3337
3338       if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3339
3340       // skip this event if near side associated particle pt larger than trigger
3341       // not really needed for calorimeter, unless DCal is included
3342
3343       Float_t deltaPhi = phiTrig-phi;
3344       if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3345       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3346
3347       if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3348
3349     }// cluster loop
3350   } // check neutral clusters
3351
3352   idLeading = index ;
3353   pLeading->SetLeadingParticle(kTRUE);
3354
3355   if( GetDebug() > 1 )
3356     printf("AliAnaParticleIsolation::IsTriggerTheNearSideEventLeadingParticle() - Particle AOD with index %d is leading with pT %2.2f\n",
3357            idLeading, pLeading->Pt());
3358
3359   return kTRUE;
3360
3361 }
3362
3363 //__________________________________________________
3364 void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
3365 {
3366   // Do analysis and fill aods
3367   // Search for the isolated photon in fCalorimeter with GetMinPt() < pt < GetMaxPt()
3368   // and if the particle is leading in the near side (if requested)
3369
3370   if(!GetInputAODBranch())
3371     AliFatal(Form("No input particles in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3372
3373   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
3374     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()));
3375
3376   Int_t n = 0, nfrac = 0;
3377   Bool_t  isolated  = kFALSE ;
3378   Float_t coneptsum = 0, coneptlead = 0;
3379   TObjArray * pl    = 0x0; ;
3380
3381   //Select the calorimeter for candidate isolation with neutral particles
3382   if      (fCalorimeter == "PHOS" )
3383     pl = GetPHOSClusters();
3384   else if (fCalorimeter == "EMCAL")
3385     pl = GetEMCALClusters();
3386
3387   //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
3388   TLorentzVector mom ;
3389   Int_t idLeading = -1 ;
3390   Int_t iaod0 = 0;
3391   Int_t naod  = GetInputAODBranch()->GetEntriesFast();
3392
3393   if(GetDebug() > 0)
3394     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
3395
3396   if(IsLeadingOnlyOn())
3397   {
3398     Bool_t leading = IsTriggerTheNearSideEventLeadingParticle(idLeading);
3399     if(!leading)
3400     {
3401       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Not leading; End fill AODs \n");
3402       return;
3403     }
3404     iaod0 = idLeading  ; // first entry in particle loop
3405     naod  = idLeading+1; // last entry in particle loop
3406   }
3407
3408   // Check isolation of list of candidate particles or leading particle
3409
3410   for(Int_t iaod = iaod0; iaod < naod; iaod++)
3411   {
3412     AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3413
3414     if(IsLeadingOnlyOn()) aodinput->SetLeadingParticle(kTRUE);
3415
3416     // Check isolation only of clusters in fiducial region
3417     if(IsFiducialCutOn())
3418     {
3419       Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodinput->Momentum(),aodinput->GetDetector()) ;
3420       if(! in ) return ;
3421     }
3422
3423     //If too small or too large pt, skip
3424     Float_t pt = aodinput->Pt();
3425     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3426
3427     //check if it is low pt trigger particle
3428     if( ( pt < GetIsolationCut()->GetPtThreshold() ||  pt < GetIsolationCut()->GetSumPtThreshold() ) &&
3429        !fMakeSeveralIC)
3430     {
3431       continue ; //trigger should not come from underlying event
3432     }
3433
3434     //After cuts, study isolation
3435     n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0; coneptlead = 0;
3436     GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,
3437                                         GetReader(), GetCaloPID(),
3438                                         kTRUE, aodinput, GetAODObjArrayName(),
3439                                         n,nfrac,coneptsum,coneptlead,isolated);
3440
3441     if(!fMakeSeveralIC) aodinput->SetIsolated(isolated);
3442   } // particle isolationloop
3443
3444   if(GetDebug() > 1)
3445   {
3446     if(isolated) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
3447     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
3448   }
3449
3450 }
3451
3452 //_________________________________________________________
3453 void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
3454 {
3455   // Do analysis and fill histograms
3456
3457   // In case of simulated data, fill acceptance histograms
3458   // Not ready for multiple case analysis.
3459   if(IsDataMC() && !fMakeSeveralIC) FillAcceptanceHistograms();
3460
3461   //Loop on stored AOD
3462   Int_t naod = GetInputAODBranch()->GetEntriesFast();
3463   if(GetDebug() > 0)
3464     printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
3465
3466   for(Int_t iaod = 0; iaod < naod ; iaod++)
3467   {
3468     AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3469
3470     if(IsLeadingOnlyOn() && !aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
3471
3472     // Check isolation only of clusters in fiducial region
3473     if(IsFiducialCutOn())
3474     {
3475       Bool_t in = GetFiducialCut()->IsInFiducialCut(*aod->Momentum(),aod->GetDetector()) ;
3476       if(! in ) continue ;
3477     }
3478
3479     Float_t pt         = aod->Pt();
3480
3481     //If too small or too large pt, skip
3482     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3483
3484     Int_t mcTag        = aod->GetTag() ;
3485     Int_t mcIndex      = GetMCIndex(mcTag);
3486
3487     // --- In case of redoing isolation from delta AOD ----
3488     // Not standard case, not used since its implementation
3489     if(fMakeSeveralIC)
3490     {
3491       //Analysis of multiple IC at same time
3492       MakeSeveralICAnalysis(aod,mcIndex);
3493       continue;
3494     }
3495
3496     // --- In case of redoing isolation multiple cuts ----
3497
3498     if(fReMakeIC)
3499     {
3500       //In case a more strict IC is needed in the produced AOD
3501       Bool_t  isolated = kFALSE;
3502       Int_t   n = 0, nfrac = 0;
3503       Float_t coneptsum = 0, coneptlead = 0;
3504
3505       //Recover reference arrays with clusters and tracks
3506       TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
3507       TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
3508
3509       GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
3510                                           GetReader(), GetCaloPID(),
3511                                           kFALSE, aod, "",
3512                                           n,nfrac,coneptsum,coneptlead,isolated);
3513     }
3514
3515     Bool_t  isolated   = aod->IsIsolated();
3516     Float_t energy     = aod->E();
3517     Float_t phi        = aod->Phi();
3518     Float_t eta        = aod->Eta();
3519
3520     Int_t   decayTag = 0;
3521     if(fFillTaggedDecayHistograms)
3522     {
3523       decayTag = aod->GetBtag(); // temporary
3524       if(decayTag < 0) decayTag = 0; // temporary
3525     }
3526
3527     if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f, Isolated %d\n",
3528                               pt, eta, phi, isolated);
3529
3530     //---------------------------------------------------------------
3531     // Fill pt/sum pT distribution of particles in cone or in UE band
3532     //---------------------------------------------------------------
3533
3534     Float_t coneptLeadCluster= 0;
3535     Float_t coneptLeadTrack  = 0;
3536     Float_t coneptsumCluster = 0;
3537     Float_t coneptsumTrack   = 0;
3538     Float_t coneptsumCell    = 0;
3539     Float_t etaBandptsumClusterNorm = 0;
3540     Float_t etaBandptsumTrackNorm   = 0;
3541
3542     CalculateTrackSignalInCone   (aod,coneptsumTrack  , coneptLeadTrack  );
3543     CalculateCaloSignalInCone    (aod,coneptsumCluster, coneptLeadCluster);
3544     if(fFillCellHistograms)
3545       CalculateCaloCellSignalInCone(aod,coneptsumCell   );
3546
3547     if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged)
3548     {
3549       fhConeSumPtClustervsTrack ->Fill(coneptsumCluster, coneptsumTrack );
3550       fhConePtLeadClustervsTrack->Fill(coneptLeadCluster,coneptLeadTrack);
3551
3552       if(coneptsumTrack  > 0) fhConeSumPtClusterTrackFrac ->Fill(pt, coneptsumCluster /coneptsumTrack );
3553       if(coneptLeadTrack > 0) fhConePtLeadClusterTrackFrac->Fill(pt, coneptLeadCluster/coneptLeadTrack);
3554
3555       if(fFillCellHistograms)
3556       {
3557         fhConeSumPtCellvsTrack        ->Fill(coneptsumCell,   coneptsumTrack);
3558         fhConeSumPtCellTrack          ->Fill(pt,     coneptsumTrack+coneptsumCell);
3559         fhConeSumPtCellTrackTrigEtaPhi->Fill(eta,phi,coneptsumTrack+coneptsumCell);
3560       }
3561     }
3562
3563     fhConeSumPt              ->Fill(pt,     coneptsumTrack+coneptsumCluster);
3564     fhConeSumPtTrigEtaPhi    ->Fill(eta,phi,coneptsumTrack+coneptsumCluster);
3565
3566     Float_t coneptLead = coneptLeadTrack;
3567     if(coneptLeadCluster > coneptLeadTrack) coneptLead = coneptLeadCluster;
3568     fhConePtLead->Fill(pt, coneptLead );
3569
3570     if(GetDebug() > 1)
3571       printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f, Leading pT in cone %2.2f\n",
3572              iaod, coneptsumTrack+coneptsumCluster, coneptLead);
3573
3574     //normalize phi/eta band per area unit
3575     if(fFillUEBandSubtractHistograms)
3576       CalculateNormalizeUEBandPerUnitArea(aod, coneptsumCluster, coneptsumCell, coneptsumTrack, etaBandptsumTrackNorm, etaBandptsumClusterNorm) ;
3577
3578     //  printf("Histograms analysis : cluster pt = %f, etaBandTrack = %f, etaBandCluster = %f, isolation = %d\n",aod->Pt(),etaBandptsumTrackNorm,etaBandptsumClusterNorm,aod->IsIsolated());
3579
3580     //---------------------------------------------------------------
3581     // Fill Shower shape and track matching histograms
3582     //---------------------------------------------------------------
3583
3584     FillTrackMatchingShowerShapeControlHistograms(aod, coneptsumTrack+coneptsumCluster, coneptLead, mcIndex);
3585
3586     //---------------------------------------------------------------
3587     // Isolated/ Non isolated histograms
3588     //---------------------------------------------------------------
3589
3590     if(isolated)
3591     {
3592       if(GetDebug() > 1)
3593         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED: fill histograms\n", iaod);
3594
3595       fhEIso      ->Fill(energy);
3596       fhPtIso     ->Fill(pt);
3597       fhPhiIso    ->Fill(pt,phi);
3598       fhEtaIso    ->Fill(pt,eta);
3599       fhEtaPhiIso ->Fill(eta,phi);
3600
3601       if(IsDataMC())
3602       {
3603         // For histograms in arrays, index in the array, corresponding to any particle origin
3604         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3605         {
3606           fhPtIsoMC [kmcPhoton]->Fill(pt);
3607           fhPhiIsoMC[kmcPhoton]->Fill(pt,phi);
3608           fhEtaIsoMC[kmcPhoton]->Fill(pt,eta);
3609         }
3610
3611         fhPtIsoMC [mcIndex]->Fill(pt);
3612         fhPhiIsoMC[mcIndex]->Fill(pt,phi);
3613         fhEtaIsoMC[mcIndex]->Fill(pt,eta);
3614       }//Histograms with MC
3615
3616       // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
3617       if(fFillTaggedDecayHistograms && decayTag > 0)
3618       {
3619         for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
3620         {
3621           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3622           {
3623             fhPtDecayIso[ibit]    ->Fill(pt);
3624             fhEtaPhiDecayIso[ibit]->Fill(eta,phi);
3625
3626             if(IsDataMC())
3627             {
3628               if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3629                 fhPtDecayIsoMC[ibit][kmcPhoton]->Fill(pt);
3630
3631               fhPtDecayIsoMC[ibit][mcIndex]->Fill(pt);
3632             }
3633           } // bit ok
3634         } // bit loop
3635       } // decay histograms
3636
3637       if(fFillNLMHistograms)
3638         fhPtNLocMaxIso ->Fill(pt,aod->GetFiducialArea()) ; // remember to change method name
3639
3640       if(fFillHighMultHistograms)
3641       {
3642         fhPtCentralityIso ->Fill(pt,GetEventCentrality()) ;
3643         fhPtEventPlaneIso ->Fill(pt,GetEventPlaneAngle() ) ;
3644       }
3645
3646       if(fFillPileUpHistograms)
3647       {
3648         if(GetReader()->IsPileUpFromSPD())               { fhEIsoPileUp[0] ->Fill(energy) ; fhPtIsoPileUp[0]->Fill(pt) ; }
3649         if(GetReader()->IsPileUpFromEMCal())             { fhEIsoPileUp[1] ->Fill(energy) ; fhPtIsoPileUp[1]->Fill(pt) ; }
3650         if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhEIsoPileUp[2] ->Fill(energy) ; fhPtIsoPileUp[2]->Fill(pt) ; }
3651         if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhEIsoPileUp[3] ->Fill(energy) ; fhPtIsoPileUp[3]->Fill(pt) ; }
3652         if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhEIsoPileUp[4] ->Fill(energy) ; fhPtIsoPileUp[4]->Fill(pt) ; }
3653         if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhEIsoPileUp[5] ->Fill(energy) ; fhPtIsoPileUp[5]->Fill(pt) ; }
3654         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhEIsoPileUp[6] ->Fill(energy) ; fhPtIsoPileUp[6]->Fill(pt) ; }
3655
3656         // Fill histograms to undertand pile-up before other cuts applied
3657         // Remember to relax time cuts in the reader
3658         FillPileUpHistograms(aod->GetCaloLabel(0));
3659       }
3660
3661     }//Isolated histograms
3662     else // NON isolated
3663     {
3664       if(GetDebug() > 1)
3665         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d NOT ISOLATED, fill histograms\n", iaod);
3666
3667       fhENoIso        ->Fill(energy);
3668       fhPtNoIso       ->Fill(pt);
3669       fhEtaPhiNoIso   ->Fill(eta,phi);
3670
3671       if(IsDataMC())
3672       {
3673         if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) )
3674           fhPtNoIsoMC[kmcPhoton]->Fill(pt);
3675
3676         fhPtNoIsoMC[mcIndex]->Fill(pt);
3677       }
3678
3679       // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
3680       if(fFillTaggedDecayHistograms && decayTag > 0)
3681       {
3682         for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
3683         {
3684           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3685           {
3686             fhPtDecayNoIso[ibit]    ->Fill(pt);
3687             fhEtaPhiDecayNoIso[ibit]->Fill(eta,phi);
3688
3689             if(IsDataMC())
3690             {
3691               if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3692                 fhPtDecayNoIsoMC[ibit][kmcPhoton]->Fill(pt);
3693
3694               fhPtDecayNoIsoMC[ibit][mcIndex]->Fill(pt);
3695             }
3696           } // bit ok
3697         } // bit loop
3698       } // decay histograms
3699
3700       if(fFillNLMHistograms)
3701         fhPtNLocMaxNoIso ->Fill(pt,aod->GetFiducialArea()); // remember to change method name
3702
3703       if(fFillPileUpHistograms)
3704       {
3705         if(GetReader()->IsPileUpFromSPD())                { fhENoIsoPileUp[0] ->Fill(energy) ; fhPtNoIsoPileUp[0]->Fill(pt) ; }
3706         if(GetReader()->IsPileUpFromEMCal())              { fhENoIsoPileUp[1] ->Fill(energy) ; fhPtNoIsoPileUp[1]->Fill(pt) ; }
3707         if(GetReader()->IsPileUpFromSPDOrEMCal())         { fhENoIsoPileUp[2] ->Fill(energy) ; fhPtNoIsoPileUp[2]->Fill(pt) ; }
3708         if(GetReader()->IsPileUpFromSPDAndEMCal())        { fhENoIsoPileUp[3] ->Fill(energy) ; fhPtNoIsoPileUp[3]->Fill(pt) ; }
3709         if(GetReader()->IsPileUpFromSPDAndNotEMCal())     { fhENoIsoPileUp[4] ->Fill(energy) ; fhPtNoIsoPileUp[4]->Fill(pt) ; }
3710         if(GetReader()->IsPileUpFromEMCalAndNotSPD())     { fhENoIsoPileUp[5] ->Fill(energy) ; fhPtNoIsoPileUp[5]->Fill(pt) ; }
3711         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhENoIsoPileUp[6] ->Fill(energy) ; fhPtNoIsoPileUp[6]->Fill(pt) ; }
3712       }
3713     } // non iso
3714   }// aod loop
3715
3716 }
3717
3718 //______________________________________________________________________
3719 void AliAnaParticleIsolation::FillAcceptanceHistograms()
3720 {
3721   // Fill acceptance histograms if MC data is available
3722   // Only when particle is in the calorimeter. Rethink if CTS is used.
3723
3724   if(GetDebug() > 0) printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - Start \n");
3725
3726   //Double_t photonY   = -100 ;
3727   Double_t photonE   = -1 ;
3728   Double_t photonPt  = -1 ;
3729   Double_t photonPhi =  100 ;
3730   Double_t photonEta = -1 ;
3731
3732   Int_t    pdg       =  0 ;
3733   Int_t    status    =  0 ;
3734   Int_t    tag       =  0 ;
3735   Int_t    mcIndex   =  0 ;
3736   Int_t    nprim     = 0;
3737
3738   TParticle        * primStack = 0;
3739   AliAODMCParticle * primAOD   = 0;
3740   TLorentzVector lv;
3741
3742   // Get the ESD MC particles container
3743   AliStack * stack = 0;
3744   if( GetReader()->ReadStack() )
3745   {
3746     stack = GetMCStack();
3747     if(!stack ) return;
3748     nprim = stack->GetNtrack();
3749   }
3750
3751   // Get the AOD MC particles container
3752   TClonesArray * mcparticles = 0;
3753   if( GetReader()->ReadAODMCParticles() )
3754   {
3755     mcparticles = GetReader()->GetAODMCParticles();
3756     if( !mcparticles ) return;
3757     nprim = mcparticles->GetEntriesFast();
3758   }
3759
3760   for(Int_t i=0 ; i < nprim; i++)
3761   {
3762     if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
3763
3764     if(GetReader()->ReadStack())
3765     {
3766       primStack = stack->Particle(i) ;
3767       if(!primStack)
3768       {
3769         printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - ESD primaries pointer not available!!\n");
3770         continue;
3771       }
3772
3773       pdg    = primStack->GetPdgCode();
3774       status = primStack->GetStatusCode();
3775
3776       if(primStack->Energy() == TMath::Abs(primStack->Pz()))  continue ; //Protection against floating point exception
3777
3778       //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
3779       //       prim->GetName(), prim->GetPdgCode());
3780
3781       //photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
3782
3783       //Photon kinematics
3784       primStack->Momentum(lv);
3785
3786     }
3787     else
3788     {
3789       primAOD = (AliAODMCParticle *) mcparticles->At(i);
3790       if(!primAOD)
3791       {
3792         printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - AOD primaries pointer not available!!\n");
3793         continue;
3794       }
3795
3796       pdg    = primAOD->GetPdgCode();
3797       status = primAOD->GetStatus();
3798
3799       if(primAOD->E() == TMath::Abs(primAOD->Pz()))  continue ; //Protection against floating point exception
3800
3801       //photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
3802
3803       //Photon kinematics
3804       lv.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
3805     }
3806
3807     // Select only photons in the final state
3808     if(pdg != 22 ) continue ;
3809
3810     // Consider only final state particles, but this depends on generator,
3811     // status 1 is the usual one, in case of not being ok, leave the possibility
3812     // to not consider this.
3813     if(status != 1 && GetMCAnalysisUtils()->GetMCGenerator()!="" ) continue ;
3814
3815     // If too small or too large pt, skip, same cut as for data analysis
3816     photonPt  = lv.Pt () ;
3817
3818     if(photonPt < GetMinPt() || photonPt > GetMaxPt() ) continue ;
3819
3820     photonE   = lv.E  () ;
3821     photonEta = lv.Eta() ;
3822     photonPhi = lv.Phi() ;
3823
3824     if(photonPhi < 0) photonPhi+=TMath::TwoPi();
3825
3826     // Check if photons hit the Calorimeter acceptance
3827     if(IsRealCaloAcceptanceOn() && fIsoDetector!="CTS") // defined on base class
3828     {
3829       if(GetReader()->ReadStack()          &&
3830          !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fIsoDetector, primStack)) continue ;
3831       if(GetReader()->ReadAODMCParticles() &&
3832          !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fIsoDetector, primAOD  )) continue ;
3833     }
3834
3835     // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
3836     if(!GetFiducialCut()->IsInFiducialCut(lv,fIsoDetector)) continue ;
3837
3838     // Get tag of this particle photon from fragmentation, decay, prompt ...
3839     // Set the origin of the photon.
3840     tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader());
3841
3842     if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
3843     {
3844       // A conversion photon from a hadron, skip this kind of photon
3845       // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
3846       // GetMCAnalysisUtils()->PrintMCTag(tag);
3847
3848       continue;
3849     }
3850
3851     //
3852     if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) )
3853     {
3854       mcIndex = kmcPrimPrompt;
3855     }
3856     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) )
3857     {
3858       mcIndex = kmcPrimFrag ;
3859     }
3860     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
3861     {
3862       mcIndex = kmcPrimISR;
3863     }
3864     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3865     {
3866       mcIndex = kmcPrimPi0Decay;
3867     }
3868     else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
3869               GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)))
3870     {
3871       mcIndex = kmcPrimOtherDecay;
3872     }
3873     else
3874     {
3875       // Other decay but from non final state particle
3876       mcIndex = kmcPrimOtherDecay;
3877     }//Other origin
3878
3879     // ////////////////////ISO MC/////////////////////////
3880     Double_t sumPtInCone = 0; Double_t dR=0. ;
3881     TParticle        * mcisopStack = 0;
3882     AliAODMCParticle * mcisopAOD   = 0;
3883     Int_t partInConeStatus = -1, partInConeMother = -1;
3884     Double_t partInConePt = 0, partInConeEta = 0, partInConePhi = 0;
3885     Int_t npart = 0;
3886     for(Int_t ip = 0; ip < nprim ; ip++)
3887     {
3888       if(ip==i) continue;
3889
3890       if( GetReader()->ReadStack() )
3891       {
3892         mcisopStack = static_cast<TParticle*>(stack->Particle(ip));
3893         if( !mcisopStack ) continue;
3894         partInConeStatus = mcisopStack->GetStatusCode();
3895         partInConeMother = mcisopStack->GetMother(0);
3896         partInConePt     = mcisopStack->Pt();
3897         partInConeEta    = mcisopStack->Eta();
3898         partInConePhi    = mcisopStack->Phi();
3899       }
3900       else
3901       {
3902         mcisopAOD   = (AliAODMCParticle *) mcparticles->At(ip);
3903         if( !mcisopAOD )   continue;
3904         partInConeStatus = mcisopAOD->GetStatus();
3905         partInConeMother = mcisopAOD->GetMother();
3906         partInConePt     = mcisopAOD->Pt();
3907         partInConeEta    = mcisopAOD->Eta();
3908         partInConePhi    = mcisopAOD->Phi();
3909       }
3910
3911       // Consider only final state particles, but this depends on generator,
3912       // status 1 is the usual one, in case of not being ok, leave the possibility
3913       // to not consider this.
3914       if( partInConeStatus != 1 && GetMCAnalysisUtils()->GetMCGenerator()!="" ) continue ;
3915
3916       if( partInConeMother == i ) continue;
3917
3918       if( partInConePt < GetReader()->GetCTSPtMin() ) continue;
3919       // Careful!!!, cut for TPC tracks and calorimeter clusters in cone can be different
3920
3921       // TVector3 vmcv(mcisop->Px(),mcisop->Py(), mcisop->Pz());
3922       // if(vmcv.Perp()>1)
3923       //   continue;
3924
3925       //
3926       // Add here Acceptance cut???, charged particles CTS fid cut, neutral Calo real acceptance.
3927       //
3928
3929       dR = GetIsolationCut()->Radius(photonEta, photonPhi, partInConeEta, partInConePhi);
3930
3931       if(dR > GetIsolationCut()->GetConeSize())
3932         continue;
3933
3934       sumPtInCone += partInConePt;
3935       if(partInConePt > GetIsolationCut()->GetPtThreshold() &&
3936          partInConePt < GetIsolationCut()->GetPtThresholdMax()) npart++;
3937     }
3938
3939     ///////END ISO MC/////////////////////////
3940
3941     // Fill the histograms, only those in the defined calorimeter acceptance
3942
3943     fhEtaPrimMC[kmcPrimPhoton]->Fill(photonPt , photonEta) ;
3944     fhPhiPrimMC[kmcPrimPhoton]->Fill(photonPt , photonPhi) ;
3945     fhEPrimMC  [kmcPrimPhoton]->Fill(photonE) ;
3946
3947     fhEtaPrimMC[mcIndex]->Fill(photonPt , photonEta) ;
3948     fhPhiPrimMC[mcIndex]->Fill(photonPt , photonPhi) ;
3949     fhEPrimMC  [mcIndex]->Fill(photonE ) ;
3950
3951     // Isolated?
3952     Bool_t isolated = kFALSE;
3953     if(GetIsolationCut()->GetICMethod() == AliIsolationCut::kSumPtIC   &&
3954        (sumPtInCone < GetIsolationCut()->GetSumPtThreshold() ||
3955         sumPtInCone > GetIsolationCut()->GetSumPtThresholdMax()))
3956       isolated = kTRUE;
3957
3958     if(GetIsolationCut()->GetICMethod() == AliIsolationCut::kPtThresIC &&
3959        npart == 0)
3960       isolated = kTRUE;
3961
3962     if(isolated)
3963     {
3964       fhPtPrimMCiso [mcIndex]      ->Fill(photonPt) ;
3965       fhPtPrimMCiso [kmcPrimPhoton]->Fill(photonPt) ;
3966     }
3967
3968   }//loop on primaries
3969
3970   if(GetDebug() > 0) printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - End \n");
3971
3972 }
3973
3974
3975 //_____________________________________________________________________________________
3976 void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph,
3977                                                      Int_t mcIndex)
3978 {
3979
3980   //Isolation Cut Analysis for both methods and different pt cuts and cones
3981   Float_t ptC   = ph->Pt();
3982   Float_t etaC  = ph->Eta();
3983   Float_t phiC  = ph->Phi();
3984   Int_t   tag   = ph->GetTag();
3985
3986   Int_t   decayTag = 0;
3987   if(fFillTaggedDecayHistograms)
3988   {
3989     decayTag = ph->GetBtag(); // temporary
3990     if(decayTag < 0) decayTag = 0; // temporary
3991   }
3992
3993   if(GetDebug() > 0)
3994     printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f, decay tag %d\n",ptC, decayTag);
3995
3996   //Keep original setting used when filling AODs, reset at end of analysis
3997   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
3998   Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();
3999   Float_t ptsumcorg  = GetIsolationCut()->GetSumPtThreshold();
4000   Float_t rorg       = GetIsolationCut()->GetConeSize();
4001
4002   Float_t coneptsum = 0, coneptlead = 0;
4003   Int_t   n    [10][10];//[fNCones][fNPtThresFrac];
4004   Int_t   nfrac[10][10];//[fNCones][fNPtThresFrac];
4005   Bool_t  isolated  = kFALSE;
4006
4007   // Fill hist with all particles before isolation criteria
4008   fhENoIso     ->Fill(ph->E());
4009   fhPtNoIso    ->Fill(ptC);
4010   fhEtaPhiNoIso->Fill(etaC,phiC);
4011
4012   if(IsDataMC())
4013   {
4014     if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4015       fhPtNoIsoMC[kmcPhoton]->Fill(ptC);
4016
4017     fhPtNoIsoMC[mcIndex]->Fill(ptC);
4018   }
4019
4020   // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
4021   if(fFillTaggedDecayHistograms && decayTag > 0)
4022   {
4023     for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
4024     {
4025       if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
4026       {
4027         fhPtDecayNoIso[ibit]    ->Fill(ptC);
4028         fhEtaPhiDecayNoIso[ibit]->Fill(etaC,phiC);
4029
4030         if(IsDataMC())
4031         {
4032           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4033             fhPtDecayNoIsoMC[ibit][kmcPhoton]->Fill(ptC);
4034
4035           fhPtDecayNoIsoMC[ibit][mcIndex]->Fill(ptC);
4036         }
4037       } // bit ok
4038     } // bit loop
4039   } // decay histograms
4040
4041   //Get vertex for photon momentum calculation
4042   Double_t vertex[] = {0,0,0} ; //vertex ;
4043   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
4044     GetReader()->GetVertex(vertex);
4045
4046   //Loop on cone sizes
4047   for(Int_t icone = 0; icone<fNCones; icone++)
4048   {
4049     //Recover reference arrays with clusters and tracks
4050     TObjArray * refclusters = ph->GetObjArray(GetAODObjArrayName()+"Clusters");
4051     TObjArray * reftracks   = ph->GetObjArray(GetAODObjArrayName()+"Tracks");
4052
4053     //If too small or too large pt, skip
4054     if(ptC < GetMinPt() || ptC > GetMaxPt() ) continue ;
4055
4056     //In case a more strict IC is needed in the produced AOD
4057
4058     isolated = kFALSE; coneptsum = 0; coneptlead = 0;
4059
4060     GetIsolationCut()->SetSumPtThreshold(100);
4061     GetIsolationCut()->SetPtThreshold(100);
4062     GetIsolationCut()->SetPtFraction(100);
4063     GetIsolationCut()->SetConeSize(fConeSizes[icone]);
4064
4065     // retreive pt tracks to fill histo vs. pt leading
4066     //Fill pt distribution of particles in cone
4067     //fhPtLeadingPt(),fhPerpSumPtLeadingPt(),fhPerpPtLeadingPt(),
4068
4069     // Tracks in perpendicular cones
4070     Double_t sumptPerp = 0. ;
4071     TObjArray * trackList   = GetCTSTracks() ;
4072     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
4073     {
4074       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
4075       //fill the histograms at forward range
4076       if(!track)
4077       {
4078         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
4079         continue;
4080       }
4081
4082       Double_t dPhi = phiC - track->Phi() + TMath::PiOver2();
4083       Double_t dEta = etaC - track->Eta();
4084       Double_t arg  = dPhi*dPhi + dEta*dEta;
4085       if(TMath::Sqrt(arg) < fConeSizes[icone])
4086       {
4087         fhPerpPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
4088         sumptPerp+=track->Pt();
4089       }
4090
4091       dPhi = phiC - track->Phi() - TMath::PiOver2();
4092       arg  = dPhi*dPhi + dEta*dEta;
4093       if(TMath::Sqrt(arg) < fConeSizes[icone])
4094       {
4095         fhPerpPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
4096         sumptPerp+=track->Pt();
4097       }
4098     }
4099
4100     fhPerpSumPtLeadingPt[icone]->Fill(ptC,sumptPerp);
4101
4102     // Tracks in isolation cone, pT distribution and sum
4103     if(reftracks && GetIsolationCut()->GetParticleTypeInCone()!= AliIsolationCut::kOnlyNeutral)
4104     {
4105       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
4106       {
4107         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
4108
4109         Float_t rad = GetIsolationCut()->Radius(etaC, phiC, track->Eta(), track->Phi());
4110
4111         if(rad > fConeSizes[icone]) continue ;
4112
4113         fhPtLeadingPt[icone]->Fill(ptC, track->Pt());
4114         coneptsum += track->Pt();
4115       }
4116     }
4117
4118     // Clusters in isolation cone, pT distribution and sum
4119     if(refclusters && GetIsolationCut()->GetParticleTypeInCone()!= AliIsolationCut::kOnlyCharged)
4120     {
4121       TLorentzVector mom ;
4122       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
4123       {
4124         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
4125
4126         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
4127
4128         Float_t rad = GetIsolationCut()->Radius(etaC, phiC, mom.Eta(), mom.Phi());
4129
4130         if(rad > fConeSizes[icone]) continue ;
4131
4132         fhPtLeadingPt[icone]->Fill(ptC, mom.Pt());
4133         coneptsum += mom.Pt();
4134       }
4135     }
4136
4137     fhSumPtLeadingPt[icone]->Fill(ptC,coneptsum);
4138
4139     if(IsDataMC())
4140     {
4141       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4142         fhSumPtLeadingPtMC[kmcPhoton][icone]->Fill(ptC,coneptsum) ;
4143
4144       fhSumPtLeadingPtMC[mcIndex][icone]->Fill(ptC,coneptsum) ;
4145     }
4146
4147     ///////////////////
4148
4149     //Loop on pt thresholds
4150     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
4151     {
4152       n    [icone][ipt]=0;
4153       nfrac[icone][ipt]=0;
4154       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
4155       GetIsolationCut()->SetPtFraction(fPtFractions[ipt]) ;
4156       GetIsolationCut()->SetSumPtThreshold(fSumPtThresholds[ipt]);
4157
4158       GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
4159                                           GetReader(), GetCaloPID(),
4160                                           kFALSE, ph, "",
4161                                           n[icone][ipt],nfrac[icone][ipt],
4162                                           coneptsum, coneptlead, isolated);
4163
4164       // Normal pT threshold cut
4165
4166       if(GetDebug() > 0)
4167       {
4168         printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - cone size %1.1f, ptThres  %1.1f, sumptThresh  %1.1f\n",
4169                fConeSizes[icone],fPtThresholds[ipt],fSumPtThresholds[ipt]);
4170         printf("\t n %d, nfrac %d, coneptsum %2.2f\n",
4171                n[icone][ipt],nfrac[icone][ipt],coneptsum);
4172
4173         printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - pt %1.1f, eta %1.1f, phi %1.1f\n",ptC, etaC, phiC);
4174       }
4175
4176       if(n[icone][ipt] == 0)
4177       {
4178         if(GetDebug() > 0)
4179           printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling pt threshold loop\n");
4180
4181         fhPtThresIsolated [icone][ipt]->Fill(ptC);
4182         fhEtaPhiPtThresIso[icone][ipt]->Fill(etaC,phiC);
4183
4184         if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4185         {
4186           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4187           {
4188             fhPtPtThresDecayIso    [icone][ipt]->Fill(ptC);
4189             fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
4190           }
4191         }
4192
4193         if(IsDataMC())
4194         {
4195           if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
4196             fhPtThresIsolatedMC[kmcPhoton][icone][ipt]->Fill(ptC) ;
4197
4198           fhPtThresIsolatedMC[mcIndex][icone][ipt]->Fill(ptC) ;
4199
4200         }
4201       }
4202
4203       // pt in cone fraction
4204       if(nfrac[icone][ipt] == 0)
4205       {
4206         if(GetDebug() > 0)
4207           printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling frac loop\n");
4208
4209         fhPtFracIsolated [icone][ipt]->Fill(ptC);
4210         fhEtaPhiPtFracIso[icone][ipt]->Fill(etaC,phiC);
4211
4212         if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4213         {
4214           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4215           {
4216             fhPtPtFracDecayIso    [icone][ipt]->Fill(ptC);
4217             fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
4218           }
4219         }
4220
4221         if(IsDataMC())
4222         {
4223           if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4224             fhPtFracIsolatedMC[kmcPhoton][icone][ipt]->Fill(ptC) ;
4225
4226           fhPtFracIsolatedMC[mcIndex][icone][ipt]->Fill(ptC) ;
4227         }
4228       }
4229
4230       if(GetDebug()>0)
4231         printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - checking IC method : %i\n",GetIsolationCut()->GetICMethod());
4232
4233       //Pt threshold on pt cand/ sum in cone histograms
4234       if(coneptsum < fSumPtThresholds[ipt])
4235       {
4236         if(GetDebug() > 0 )
4237           printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling sum loop\n");
4238
4239         fhSumPtIsolated [icone][ipt]->Fill(ptC) ;
4240         fhEtaPhiPtSumIso[icone][ipt]->Fill(etaC, phiC) ;
4241
4242         if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4243         {
4244           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4245           {
4246             fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
4247             fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
4248           }
4249         }
4250
4251         if(IsDataMC())
4252         {
4253           if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4254             fhSumPtIsolatedMC[kmcPhoton][icone][ipt]->Fill(ptC) ;
4255
4256           fhSumPtIsolatedMC[mcIndex][icone][ipt]->Fill(ptC) ;
4257         }
4258       }
4259
4260       // pt sum pt frac method
4261       //    if( ((fPtFractions[ipt]*ptC < fSumPtThresholds[ipt]) && (coneptsum < fSumPtThresholds[ipt])) || ((fPtFractions[ipt]*ptC > fSumPtThresholds[ipt]) && (coneptsum < fPtFractions[ipt]*ptC)) )
4262
4263       if(coneptsum < fPtFractions[ipt]*ptC)
4264       {
4265         if(GetDebug() > 0)
4266           printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling PtFrac PtSum loop\n");
4267
4268         fhPtFracPtSumIso    [icone][ipt]->Fill(ptC) ;
4269         fhEtaPhiFracPtSumIso[icone][ipt]->Fill(etaC,phiC) ;
4270
4271         if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4272         {
4273           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4274           {
4275             fhPtFracPtSumDecayIso    [icone][ipt]->Fill(ptC);
4276             fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
4277           }
4278         }
4279       }
4280
4281       // density method
4282       Float_t cellDensity = GetIsolationCut()->GetCellDensity( ph, GetReader());
4283       if(coneptsum < fSumPtThresholds[ipt]*cellDensity)
4284       {
4285         if(GetDebug() > 0)
4286           printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling density loop\n");
4287
4288         fhPtSumDensityIso    [icone][ipt]->Fill(ptC) ;
4289         fhEtaPhiSumDensityIso[icone][ipt]->Fill(etaC,phiC) ;
4290
4291         if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4292         {
4293           if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4294           {
4295             fhPtSumDensityDecayIso    [icone][ipt]->Fill(ptC);
4296             fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
4297           }
4298         }
4299       }
4300     }//pt thresh loop
4301
4302
4303   }//cone size loop
4304
4305   //Reset original parameters for AOD analysis
4306   GetIsolationCut()->SetPtThreshold(ptthresorg);
4307   GetIsolationCut()->SetPtFraction(ptfracorg);
4308   GetIsolationCut()->SetSumPtThreshold(ptsumcorg);
4309   GetIsolationCut()->SetConeSize(rorg);
4310
4311 }
4312
4313 //_____________________________________________________________
4314 void AliAnaParticleIsolation::Print(const Option_t * opt) const
4315 {
4316
4317   //Print some relevant parameters set for the analysis
4318   if(! opt)
4319     return;
4320
4321   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
4322   AliAnaCaloTrackCorrBaseClass::Print(" ");
4323
4324   printf("ReMake Isolation          = %d \n",  fReMakeIC) ;
4325   printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;
4326   printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;
4327   printf("Detector for candidate isolation = %s \n", fIsoDetector.Data()) ;
4328
4329   if(fMakeSeveralIC)
4330   {
4331     printf("N Cone Sizes       =     %d\n", fNCones) ;
4332     printf("Cone Sizes          =    \n") ;
4333     for(Int_t i = 0; i < fNCones; i++)
4334       printf("  %1.2f;",  fConeSizes[i]) ;
4335     printf("    \n") ;
4336
4337     printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
4338     printf(" pT thresholds         =    \n") ;
4339     for(Int_t i = 0; i < fNPtThresFrac; i++)
4340       printf("   %2.2f;",  fPtThresholds[i]) ;
4341
4342     printf("    \n") ;
4343
4344     printf(" pT fractions         =    \n") ;
4345     for(Int_t i = 0; i < fNPtThresFrac; i++)
4346       printf("   %2.2f;",  fPtFractions[i]) ;
4347
4348     printf("    \n") ;
4349
4350     printf("sum pT thresholds         =    \n") ;
4351     for(Int_t i = 0; i < fNPtThresFrac; i++)
4352       printf("   %2.2f;",  fSumPtThresholds[i]) ;
4353
4354
4355   }
4356
4357   printf("    \n") ;
4358
4359 }
4360