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