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