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