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