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