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