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
  36 // --- Analysis system ---
  37 #include "AliAnaParticleIsolation.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliIsolationCut.h"
  40 #include "AliNeutralMesonSelection.h"
  41 #include "AliAODPWG4ParticleCorrelation.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "AliVTrack.h"
  44 #include "AliVCluster.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47
  48 ClassImp(AliAnaParticleIsolation)
  49
  50 //______________________________________________________________________________
  51 AliAnaParticleIsolation::AliAnaParticleIsolation() :
  52 AliAnaCaloTrackCorrBaseClass(),   fCalorimeter(""),
  53 fReMakeIC(0),                     fMakeSeveralIC(0),
  54 fFillPileUpHistograms(0),
  55 fFillTMHisto(0),                  fFillSSHisto(0),
  56 // Several IC
  57 fNCones(0),                       fNPtThresFrac(0),
  58 fConeSizes(),                     fPtThresholds(),
  59 fPtFractions(),                   fSumPtThresholds(),
  60 // Histograms
  61 fhEIso(0),                        fhPtIso(0),
  62 fhPtCentralityIso(0),             fhPtEventPlaneIso(0),
  63 fhPtNLocMaxIso(0),
  64 fhPhiIso(0),                      fhEtaIso(0),                     fhEtaPhiIso(0),
  65 fhEtaPhiNoIso(0),
  66 fhENoIso(0),                      fhPtNoIso(0),                    fhPtNLocMaxNoIso(0),
  67 fhPtDecayIso(0),                  fhPtDecayNoIso(0),
  68 fhEtaPhiDecayIso(0),              fhEtaPhiDecayNoIso(0),
  69 fhConeSumPt(0),                   fhPtInCone(0),
  70 fhPtTrackInCone(0),
  71 fhPtTrackInConeOtherBC(0),        fhPtTrackInConeOtherBCPileUpSPD(0),
  72 fhPtTrackInConeBC0(0),            fhPtTrackInConeVtxBC0(0),
  73 fhPtTrackInConeBC0PileUpSPD(0),
  74 fhPtInConePileUp(),               fhPtInConeCent(0),
  75 fhFRConeSumPt(0),                 fhPtInFRCone(0),                 fhPhiUEConeSumPt(0),
  76 fhEtaUEConeSumPt(0),              fhEtaBand(0),                    fhPhiBand(0),
  77 fhConeSumPtEtaUESub(0),           fhConeSumPtPhiUESub(0),
  78 // MC histograms
  79 fhPtIsoPrompt(0),                 fhPhiIsoPrompt(0),               fhEtaIsoPrompt(0),
  80 fhPtThresIsolatedPrompt(),        fhPtFracIsolatedPrompt(),        fhPtSumIsolatedPrompt(),
  81 fhPtIsoFragmentation(0),          fhPhiIsoFragmentation(0),        fhEtaIsoFragmentation(0),
  82 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
  83 fhPtIsoPi0(0),                    fhPhiIsoPi0(0),                  fhEtaIsoPi0(0),
  84 fhPtThresIsolatedPi0(),           fhPtFracIsolatedPi0(),           fhPtSumIsolatedPi0(),
  85 fhPtIsoPi0Decay(0),               fhPhiIsoPi0Decay(0),             fhEtaIsoPi0Decay(0),
  86 fhPtThresIsolatedPi0Decay(),      fhPtFracIsolatedPi0Decay(),      fhPtSumIsolatedPi0Decay(),
  87 fhPtIsoEtaDecay(0),               fhPhiIsoEtaDecay(0),             fhEtaIsoEtaDecay(0),
  88 fhPtThresIsolatedEtaDecay(),      fhPtFracIsolatedEtaDecay(),      fhPtSumIsolatedEtaDecay(),
  89 fhPtIsoOtherDecay(0),             fhPhiIsoOtherDecay(0),           fhEtaIsoOtherDecay(0),
  90 fhPtThresIsolatedOtherDecay(),    fhPtFracIsolatedOtherDecay(),    fhPtSumIsolatedOtherDecay(),
  91 //fhPtIsoConversion(0),             fhPhiIsoConversion(0),           fhEtaIsoConversion(0),
  92 //fhPtThresIsolatedConversion(),    fhPtFracIsolatedConversion(),    fhPtSumIsolatedConversion(),
  93 fhPtIsoHadron(0),                 fhPhiIsoHadron(0),               fhEtaIsoHadron(0),
  94 fhPtThresIsolatedHadron(),        fhPtFracIsolatedHadron(),        fhPtSumIsolatedHadron(),
  95 fhPtNoIsoPi0(0),                  fhPtNoIsoPi0Decay(0),
  96 fhPtNoIsoEtaDecay(0),             fhPtNoIsoOtherDecay(0),
  97 fhPtNoIsoPrompt(0),               fhPtIsoMCPhoton(0),              fhPtNoIsoMCPhoton(0),
  98 //fhPtNoIsoConversion(0),
  99 fhPtNoIsoFragmentation(0),        fhPtNoIsoHadron(0),
 100 // Hist several IC
 101 fhSumPtLeadingPt(),               fhPtLeadingPt(),
 102 fhFRSumPtLeadingPt(),             fhFRPtLeadingPt(),
 103 fhPtThresIsolated(),              fhPtFracIsolated(),              fhPtSumIsolated(),
 104 fhEtaPhiPtThresIso(),             fhEtaPhiPtThresDecayIso(),       fhPtPtThresDecayIso(),
 105 fhEtaPhiPtFracIso(),              fhEtaPhiPtFracDecayIso(),        fhPtPtFracDecayIso(),
 106 fhPtPtSumDecayIso(),              fhEtaPhiSumDensityIso(),         fhEtaPhiSumDensityDecayIso(),
 107 fhPtSumDensityIso(),              fhPtSumDensityDecayIso(),
 108 fhPtFracPtSumIso(),               fhPtFracPtSumDecayIso(),
 109 fhEtaPhiFracPtSumIso(),           fhEtaPhiFracPtSumDecayIso(),
 110 // Cluster control histograms
 111 fhTrackMatchedDEta(),             fhTrackMatchedDPhi(),           fhTrackMatchedDEtaDPhi(),
 112 fhdEdx(),                         fhEOverP(),                     fhTrackMatchedMCParticle(),
 113 fhELambda0() ,                    fhELambda1(),                   fhELambda0SSBkg(),
 114 fhELambda0TRD(),                  fhELambda1TRD(),
 115 fhELambda0MCPhoton(),             fhELambda0MCPi0(),              fhELambda0MCPi0Decay(),
 116 fhELambda0MCEtaDecay(),           fhELambda0MCOtherDecay(),       fhELambda0MCHadron(),
 117 // Number of local maxima in cluster
 118 fhNLocMax(),
 119 fhELambda0LocMax1(),              fhELambda1LocMax1(),
 120 fhELambda0LocMax2(),              fhELambda1LocMax2(),
 121 fhELambda0LocMaxN(),              fhELambda1LocMaxN(),
 122 // PileUp
 123 fhEIsoPileUp(),                   fhPtIsoPileUp(),
 124 fhENoIsoPileUp(),                 fhPtNoIsoPileUp(),
 125 fhTimeENoCut(0),                  fhTimeESPD(0),                  fhTimeESPDMulti(0),
 126 fhTimeNPileUpVertSPD(0),          fhTimeNPileUpVertTrack(0),
 127 fhTimeNPileUpVertContributors(0),
 128 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
 129 // Histograms settings
 130 fHistoNPtSumBins(0),              fHistoPtSumMax(0.),              fHistoPtSumMin(0.),
 131 fHistoNPtInConeBins(0),           fHistoPtInConeMax(0.),           fHistoPtInConeMin(0.)
 132 {
 133   //default ctor
 134
 135   //Initialize parameters
 136   InitParameters();
 137
 138   for(Int_t i = 0; i < 5 ; i++)
 139   {
 140     fConeSizes[i]      = 0 ;
 141
 142     fhPtSumIsolatedPrompt       [i] = 0 ;
 143     fhPtSumIsolatedFragmentation[i] = 0 ;
 144     fhPtSumIsolatedPi0Decay     [i] = 0 ;
 145     fhPtSumIsolatedPi0          [i] = 0 ;
 146     fhPtSumIsolatedEtaDecay     [i] = 0 ;
 147     fhPtSumIsolatedOtherDecay   [i] = 0 ;
 148 //    fhPtSumIsolatedConversion   [i] = 0 ;
 149     fhPtSumIsolatedHadron      [i] = 0 ;
 150
 151     for(Int_t j = 0; j < 5 ; j++)
 152     {
 153       fhPtThresIsolated      [i][j] = 0 ;
 154       fhPtFracIsolated       [i][j] = 0 ;
 155       fhPtSumIsolated        [i][j] = 0 ;
 156
 157       fhEtaPhiPtThresIso     [i][j] = 0 ;
 158       fhEtaPhiPtThresDecayIso[i][j] = 0 ;
 159       fhPtPtThresDecayIso    [i][j] = 0 ;
 160
 161       fhEtaPhiPtFracIso      [i][j] = 0 ;
 162       fhEtaPhiPtFracDecayIso [i][j] = 0 ;
 163       fhPtPtFracDecayIso     [i][j] = 0 ;
 164       fhPtPtSumDecayIso      [i][j] = 0 ;
 165       fhPtSumDensityIso      [i][j] = 0 ;
 166       fhPtSumDensityDecayIso [i][j] = 0 ;
 167       fhEtaPhiSumDensityIso      [i][j] = 0 ;
 168       fhEtaPhiSumDensityDecayIso [i][j] = 0 ;
 169       fhPtFracPtSumIso           [i][j] = 0 ;
 170       fhPtFracPtSumDecayIso      [i][j] = 0 ;
 171       fhEtaPhiFracPtSumIso       [i][j] = 0 ;
 172       fhEtaPhiFracPtSumDecayIso  [i][j] = 0 ;
 173
 174       fhPtThresIsolatedPrompt       [i][j] = 0 ;
 175       fhPtThresIsolatedFragmentation[i][j] = 0 ;
 176       fhPtThresIsolatedPi0Decay     [i][j] = 0 ;
 177       fhPtThresIsolatedPi0          [i][j] = 0 ;
 178       fhPtThresIsolatedEtaDecay     [i][j] = 0 ;
 179       fhPtThresIsolatedOtherDecay   [i][j] = 0 ;
 180 //      fhPtThresIsolatedConversion   [i][j] = 0 ;
 181       fhPtThresIsolatedHadron      [i][j] = 0 ;
 182
 183       fhPtFracIsolatedPrompt        [i][j] = 0 ;
 184       fhPtFracIsolatedFragmentation [i][j] = 0 ;
 185       fhPtFracIsolatedPi0           [i][j] = 0 ;
 186       fhPtFracIsolatedPi0Decay      [i][j] = 0 ;
 187       fhPtFracIsolatedEtaDecay      [i][j] = 0 ;
 188       fhPtFracIsolatedOtherDecay    [i][j] = 0 ;
 189 //      fhPtFracIsolatedConversion    [i][j] = 0 ;
 190       fhPtFracIsolatedHadron       [i][j] = 0 ;
 191
 192     }
 193   }
 194
 195   for(Int_t i = 0; i < 5 ; i++)
 196   {
 197     fPtFractions    [i] = 0 ;
 198     fPtThresholds   [i] = 0 ;
 199     fSumPtThresholds[i] = 0 ;
 200   }
 201
 202
 203   for(Int_t i = 0; i < 2 ; i++)
 204   {
 205     fhTrackMatchedDEta[i] = 0 ;             fhTrackMatchedDPhi[i] = 0 ;           fhTrackMatchedDEtaDPhi  [i] = 0 ;
 206     fhdEdx            [i] = 0 ;             fhEOverP          [i] = 0 ;           fhTrackMatchedMCParticle[i] = 0 ;
 207     fhELambda0        [i] = 0 ;             fhELambda1        [i] = 0 ;
 208     fhELambda0TRD     [i] = 0 ;             fhELambda1TRD     [i] = 0 ;
 209
 210     fhELambda0MCPhoton  [i] = 0 ;           fhELambda0MCPi0       [i] = 0 ;       fhELambda0MCPi0Decay[i] = 0 ;
 211     fhELambda0MCEtaDecay[i] = 0 ;           fhELambda0MCOtherDecay[i] = 0 ;       fhELambda0MCHadron  [i] = 0 ;
 212
 213
 214     // Number of local maxima in cluster
 215     fhNLocMax        [i] = 0 ;
 216     fhELambda0LocMax1[i] = 0 ;              fhELambda1LocMax1[i] = 0 ;
 217     fhELambda0LocMax2[i] = 0 ;              fhELambda1LocMax2[i] = 0 ;
 218     fhELambda0LocMaxN[i] = 0 ;              fhELambda1LocMaxN[i] = 0 ;
 219
 220   }
 221
 222   // Pile-Up
 223
 224   for(Int_t i = 0 ; i < 7 ; i++)
 225   {
 226     fhPtInConePileUp[i] = 0 ;
 227     fhEIsoPileUp    [i] = 0 ;
 228     fhPtIsoPileUp   [i] = 0 ;
 229     fhENoIsoPileUp  [i] = 0 ;
 230     fhPtNoIsoPileUp [i] = 0 ;
 231   }
 232
 233 }
 234
 235 //_________________________________________________________________
 236 void AliAnaParticleIsolation::FillPileUpHistograms(Int_t clusterID)
 237 {
 238   // Fill some histograms to understand pile-up
 239   if(!fFillPileUpHistograms) return;
 240
 241   if(clusterID < 0 )
 242   {
 243     printf("AliAnaParticleIsolation::FillPileUpHistograms(), ID of cluster = %d, not possible! ", clusterID);
 244     return;
 245   }
 246
 247   Int_t iclus = -1;
 248   TObjArray* clusters = 0x0;
 249   if     (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
 250   else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
 251
 252   Float_t energy = 0;
 253   Float_t time   = -1000;
 254
 255   if(clusters)
 256   {
 257     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
 258     energy = cluster->E();
 259     time   = cluster->GetTOF()*1e9;
 260   }
 261
 262   //printf("E %f, time %f\n",energy,time);
 263   AliVEvent * event = GetReader()->GetInputEvent();
 264
 265   fhTimeENoCut->Fill(energy,time);
 266   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
 267   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
 268
 269   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
 270
 271   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
 272   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
 273
 274   // N pile up vertices
 275   Int_t nVerticesSPD    = -1;
 276   Int_t nVerticesTracks = -1;
 277
 278   if      (esdEv)
 279   {
 280     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
 281     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
 282
 283   }//ESD
 284   else if (aodEv)
 285   {
 286     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
 287     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
 288   }//AOD
 289
 290   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
 291   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
 292
 293   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
 294   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
 295
 296   Int_t ncont = -1;
 297   Float_t z1 = -1, z2 = -1;
 298   Float_t diamZ = -1;
 299   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
 300   {
 301     if      (esdEv)
 302     {
 303       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
 304       ncont=pv->GetNContributors();
 305       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
 306       z2 = pv->GetZ();
 307       diamZ = esdEv->GetDiamondZ();
 308     }//ESD
 309     else if (aodEv)
 310     {
 311       AliAODVertex *pv=aodEv->GetVertex(iVert);
 312       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
 313       ncont=pv->GetNContributors();
 314       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
 315       z2=pv->GetZ();
 316       diamZ = aodEv->GetDiamondZ();
 317     }// AOD
 318
 319     Double_t distZ  = TMath::Abs(z2-z1);
 320     diamZ  = TMath::Abs(z2-diamZ);
 321
 322     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
 323     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
 324     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
 325
 326   }// loop
 327 }
 328
 329 //________________________________________________________________________________________________
 330 void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(const Bool_t isolated,
 331                                                                             const Int_t  clusterID,
 332                                                                             const Int_t  nMaxima,
 333                                                                             const Int_t  mcTag,
 334                                                                             const TObjArray * plCTS,
 335                                                                             const TObjArray * plNe,
 336                                                                             AliAODPWG4ParticleCorrelation  *pCandidate,
 337                                                                             const AliCaloTrackReader * reader,
 338                                                                             const AliCaloPID * pid)
 339 {
 340   // Fill Track matching and Shower Shape control histograms
 341   if(!fFillTMHisto &&  !fFillSSHisto) return;
 342
 343   if(clusterID < 0 )
 344   {
 345     printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(), ID of cluster = %d, not possible! \n", clusterID);
 346     return;
 347   }
 348
 349   Int_t iclus = -1;
 350   TObjArray* clusters = 0x0;
 351   if     (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
 352   else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
 353
 354   if(clusters)
 355   {
 356
 357     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
 358     Float_t energy = cluster->E();
 359
 360     if(fFillSSHisto)
 361     {
 362       fhELambda0[isolated]->Fill(energy, cluster->GetM02() );
 363       fhELambda1[isolated]->Fill(energy, cluster->GetM20() );
 364
 365       if(IsDataMC())
 366       {
 367         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||
 368                 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) fhELambda0MCPhoton    [isolated]->Fill(energy, cluster->GetM02());
 369         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))           fhELambda0MCPi0       [isolated]->Fill(energy, cluster->GetM02());
 370         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      fhELambda0MCPi0Decay  [isolated]->Fill(energy, cluster->GetM02());
 371         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      fhELambda0MCEtaDecay  [isolated]->Fill(energy, cluster->GetM02());
 372         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    fhELambda0MCOtherDecay[isolated]->Fill(energy, cluster->GetM02());
 373
 374         //        else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(energy, cluster->GetM02());
 375         else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))     fhELambda0MCHadron    [isolated]->Fill(energy, cluster->GetM02());
 376
 377       }
 378
 379       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5) // TO DO: CHANGE FOR 2012
 380       {
 381         fhELambda0TRD[isolated]->Fill(energy, cluster->GetM02() );
 382         fhELambda1TRD[isolated]->Fill(energy, cluster->GetM20() );
 383       }
 384
 385       fhNLocMax[isolated]->Fill(energy,nMaxima);
 386       if     (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax1[isolated]->Fill(energy,cluster->GetM20()); }
 387       else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax2[isolated]->Fill(energy,cluster->GetM20()); }
 388       else                { fhELambda0LocMaxN[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMaxN[isolated]->Fill(energy,cluster->GetM20()); }
 389
 390       if(isolated==0)
 391         {
 392           //Analyse non-isolated events
 393           Int_t   n = 0;
 394           Int_t nfrac = 0;
 395           Bool_t  iso  = kFALSE ;
 396           Float_t coneptsum = 0 ;
 397           GetIsolationCut()->SetPtThresholdMax(1.);
 398           GetIsolationCut()->MakeIsolationCut(plCTS,   plNe,
 399                                               reader, pid,
 400                                               kFALSE, pCandidate, "",
 401                                               n,nfrac,coneptsum, iso);
 402           if (!iso) fhELambda0SSBkg->Fill(energy, cluster->GetM02());
 403
 404
 405           if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
 406         }
 407       GetIsolationCut()->SetPtThresholdMax(10000.);
 408
 409     } // SS histo fill
 410
 411
 412     if(fFillTMHisto)
 413     {
 414       Float_t dZ  = cluster->GetTrackDz();
 415       Float_t dR  = cluster->GetTrackDx();
 416
 417       if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 418       {
 419         dR = 2000., dZ = 2000.;
 420         GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 421       }
 422
 423       //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
 424       if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
 425       {
 426         fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
 427         fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
 428         if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
 429       }
 430
 431       // Check dEdx and E/p of matched clusters
 432
 433       if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 434       {
 435
 436         AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
 437
 438         if(track)
 439         {
 440           Float_t dEdx = track->GetTPCsignal();
 441           fhdEdx[isolated]->Fill(cluster->E(), dEdx);
 442
 443           Float_t eOverp = cluster->E()/track->P();
 444           fhEOverP[isolated]->Fill(cluster->E(),  eOverp);
 445         }
 446         //else
 447         //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
 448
 449
 450         if(IsDataMC())
 451         {
 452           if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
 453           {
 454             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
 455                        GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
 456             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
 457             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
 458             else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
 459
 460           }
 461           else
 462           {
 463             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
 464                        GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
 465             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
 466             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
 467             else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
 468           }
 469
 470         }  // MC
 471
 472       } // match window
 473
 474     }// TM histos fill
 475
 476   } // clusters array available
 477
 478 }
 479
 480 //______________________________________________________
 481 TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
 482 {
 483   //Save parameters used for analysis
 484   TString parList ; //this will be list of parameters used for this analysis.
 485   const Int_t buffersize = 255;
 486   char onePar[buffersize] ;
 487
 488   snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
 489   parList+=onePar ;
 490   snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 491   parList+=onePar ;
 492   snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
 493   parList+=onePar ;
 494   snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
 495   parList+=onePar ;
 496   snprintf(onePar, buffersize,"fFillTMHisto=%d (Flag for track matching histograms) \n",fFillTMHisto) ;
 497   parList+=onePar ;
 498   snprintf(onePar, buffersize,"fFillSSHisto=%d (Flag for shower shape histograms) \n",fFillSSHisto) ;
 499   parList+=onePar ;
 500
 501   if(fMakeSeveralIC)
 502   {
 503     snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
 504     parList+=onePar ;
 505     snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
 506     parList+=onePar ;
 507
 508     for(Int_t icone = 0; icone < fNCones ; icone++)
 509     {
 510       snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
 511       parList+=onePar ;
 512     }
 513     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
 514     {
 515       snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
 516       parList+=onePar ;
 517     }
 518     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
 519     {
 520       snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
 521       parList+=onePar ;
 522     }
 523     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
 524     {
 525       snprintf(onePar, buffersize,"fSumPtThresholds[%d]=%1.2f (isolation sum pt threshold) \n",ipt, fSumPtThresholds[ipt]) ;
 526       parList+=onePar ;
 527     }
 528   }
 529
 530   //Get parameters set in base class.
 531   parList += GetBaseParametersList() ;
 532
 533   //Get parameters set in IC class.
 534   if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
 535
 536   return new TObjString(parList) ;
 537
 538 }
 539
 540 //________________________________________________________
 541 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
 542 {
 543   // Create histograms to be saved in output file and
 544   // store them in outputContainer
 545   TList * outputContainer = new TList() ;
 546   outputContainer->SetName("IsolatedParticleHistos") ;
 547
 548   Int_t   nptbins  = GetHistogramRanges()->GetHistoPtBins();
 549   Int_t   nphibins = GetHistogramRanges()->GetHistoPhiBins();
 550   Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
 551   Float_t ptmax    = GetHistogramRanges()->GetHistoPtMax();
 552   Float_t phimax   = GetHistogramRanges()->GetHistoPhiMax();
 553   Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
 554   Float_t ptmin    = GetHistogramRanges()->GetHistoPtMin();
 555   Float_t phimin   = GetHistogramRanges()->GetHistoPhiMin();
 556   Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
 557   Int_t   ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();
 558   Float_t ssmax    = GetHistogramRanges()->GetHistoShowerShapeMax();
 559   Float_t ssmin    = GetHistogramRanges()->GetHistoShowerShapeMin();
 560   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
 561   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
 562   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
 563
 564   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 565   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 566   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 567   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 568   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 569   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 570
 571   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
 572   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
 573   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
 574   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
 575   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
 576   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
 577
 578   Int_t   nptsumbins    = fHistoNPtSumBins;
 579   Float_t ptsummax      = fHistoPtSumMax;
 580   Float_t ptsummin      = fHistoPtSumMin;
 581   Int_t   nptinconebins = fHistoNPtInConeBins;
 582   Float_t ptinconemax   = fHistoPtInConeMax;
 583   Float_t ptinconemin   = fHistoPtInConeMin;
 584
 585   Float_t ptthre = GetIsolationCut()->GetPtThreshold();
 586   Float_t ptfrac = GetIsolationCut()->GetPtFraction();
 587   Float_t r      = GetIsolationCut()->GetConeSize();
 588
 589   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
 590
 591   if(!fMakeSeveralIC)
 592   {
 593     TString hName [] = {"NoIso",""};
 594     TString hTitle[] = {"Not isolated"  ,"isolated"};
 595
 596     if(fFillSSHisto)
 597     {
 598         fhELambda0SSBkg  = new TH2F
 599       ("hELambda0SSBkg","Non isolated clusters : E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 600       fhELambda0SSBkg->SetYTitle("#lambda_{0}^{2}");
 601       fhELambda0SSBkg->SetXTitle("E (GeV)");
 602       outputContainer->Add(fhELambda0SSBkg) ;
 603     }
 604
 605     for(Int_t iso = 0; iso < 2; iso++)
 606     {
 607       if(fFillTMHisto)
 608       {
 609         fhTrackMatchedDEta[iso]  = new TH2F
 610         (Form("hTrackMatchedDEta%s",hName[iso].Data()),
 611          Form("%s - d#eta of cluster-track vs cluster energy for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
 612          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 613         fhTrackMatchedDEta[iso]->SetYTitle("d#eta");
 614         fhTrackMatchedDEta[iso]->SetXTitle("E_{cluster} (GeV)");
 615
 616         fhTrackMatchedDPhi[iso]  = new TH2F
 617         (Form("hTrackMatchedDPhi%s",hName[iso].Data()),
 618          Form("%s - d#phi of cluster-track vs cluster energy for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
 619          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 620         fhTrackMatchedDPhi[iso]->SetYTitle("d#phi (rad)");
 621         fhTrackMatchedDPhi[iso]->SetXTitle("E_{cluster} (GeV)");
 622
 623         fhTrackMatchedDEtaDPhi[iso]  = new TH2F
 624         (Form("hTrackMatchedDEtaDPhi%s",hName[iso].Data()),
 625          Form("%s - d#eta vs d#phi of cluster-track for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
 626          nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
 627         fhTrackMatchedDEtaDPhi[iso]->SetYTitle("d#phi (rad)");
 628         fhTrackMatchedDEtaDPhi[iso]->SetXTitle("d#eta");
 629
 630         outputContainer->Add(fhTrackMatchedDEta[iso]) ;
 631         outputContainer->Add(fhTrackMatchedDPhi[iso]) ;
 632         outputContainer->Add(fhTrackMatchedDEtaDPhi[iso]) ;
 633
 634         fhdEdx[iso]  = new TH2F
 635         (Form("hdEdx%s",hName[iso].Data()),
 636          Form("%s - Matched track <dE/dx> vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
 637          nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 638         fhdEdx[iso]->SetXTitle("E (GeV)");
 639         fhdEdx[iso]->SetYTitle("<dE/dx>");
 640         outputContainer->Add(fhdEdx[iso]);
 641
 642         fhEOverP[iso]  = new TH2F
 643         (Form("hEOverP%s",hName[iso].Data()),
 644          Form("%s - Matched track E/p vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
 645          nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 646         fhEOverP[iso]->SetXTitle("E (GeV)");
 647         fhEOverP[iso]->SetYTitle("E/p");
 648         outputContainer->Add(fhEOverP[iso]);
 649
 650         if(IsDataMC())
 651         {
 652           fhTrackMatchedMCParticle[iso]  = new TH2F
 653           (Form("hTrackMatchedMCParticle%s",hName[iso].Data()),
 654            Form("%s - Origin of particle vs energy vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",hTitle[iso].Data(),r,ptthre,ptfrac),
 655            nptbins,ptmin,ptmax,8,0,8);
 656           fhTrackMatchedMCParticle[iso]->SetXTitle("E (GeV)");
 657           //fhTrackMatchedMCParticle[iso]->SetYTitle("Particle type");
 658
 659           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(1 ,"Photon");
 660           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(2 ,"Electron");
 661           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
 662           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(4 ,"Rest");
 663           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
 664           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
 665           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
 666           fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
 667
 668           outputContainer->Add(fhTrackMatchedMCParticle[iso]);
 669         }
 670       }
 671
 672       if(fFillSSHisto)
 673       {
 674         fhELambda0[iso]  = new TH2F
 675         (Form("hELambda0%s",hName[iso].Data()),
 676          Form("%s cluster : E vs #lambda_{0}",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 677         fhELambda0[iso]->SetYTitle("#lambda_{0}^{2}");
 678         fhELambda0[iso]->SetXTitle("E (GeV)");
 679         outputContainer->Add(fhELambda0[iso]) ;
 680
 681         if(IsDataMC())
 682         {
 683           fhELambda0MCPhoton[iso]  = new TH2F
 684           (Form("hELambda0%s_MCPhoton",hName[iso].Data()),
 685            Form("%s cluster : E vs #lambda_{0}: Origin is final state photon",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 686           fhELambda0MCPhoton[iso]->SetYTitle("#lambda_{0}^{2}");
 687           fhELambda0MCPhoton[iso]->SetXTitle("E (GeV)");
 688           outputContainer->Add(fhELambda0MCPhoton[iso]) ;
 689
 690           fhELambda0MCPi0[iso]  = new TH2F
 691           (Form("hELambda0%s_MCPi0",hName[iso].Data()),
 692            Form("%s cluster : E vs #lambda_{0}: Origin is pi0 (2 #gamma)",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 693           fhELambda0MCPi0[iso]->SetYTitle("#lambda_{0}^{2}");
 694           fhELambda0MCPi0[iso]->SetXTitle("E (GeV)");
 695           outputContainer->Add(fhELambda0MCPi0[iso]) ;
 696
 697           fhELambda0MCPi0Decay[iso]  = new TH2F
 698           (Form("hELambda0%s_MCPi0Decay",hName[iso].Data()),
 699            Form("%s cluster : E vs #lambda_{0}: Origin is pi0 decay",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 700           fhELambda0MCPi0Decay[iso]->SetYTitle("#lambda_{0}^{2}");
 701           fhELambda0MCPi0Decay[iso]->SetXTitle("E (GeV)");
 702           outputContainer->Add(fhELambda0MCPi0Decay[iso]) ;
 703
 704           fhELambda0MCEtaDecay[iso]  = new TH2F
 705           (Form("hELambda0%s_MCEtaDecay",hName[iso].Data()),
 706            Form("%s cluster : E vs #lambda_{0}: Origin is eta decay",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 707           fhELambda0MCEtaDecay[iso]->SetYTitle("#lambda_{0}^{2}");
 708           fhELambda0MCEtaDecay[iso]->SetXTitle("E (GeV)");
 709           outputContainer->Add(fhELambda0MCEtaDecay[iso]) ;
 710
 711           fhELambda0MCOtherDecay[iso]  = new TH2F
 712           (Form("hELambda0%s_MCOtherDecay",hName[iso].Data()),
 713            Form("%s cluster : E vs #lambda_{0}: Origin is other decay",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 714           fhELambda0MCOtherDecay[iso]->SetYTitle("#lambda_{0}^{2}");
 715           fhELambda0MCOtherDecay[iso]->SetXTitle("E (GeV)");
 716           outputContainer->Add(fhELambda0MCOtherDecay[iso]) ;
 717
 718           fhELambda0MCHadron[iso]  = new TH2F
 719           (Form("hELambda0%s_MCHadron",hName[iso].Data()),
 720            Form("%s cluster : E vs #lambda_{0}: Origin is hadron",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 721           fhELambda0MCHadron[iso]->SetYTitle("#lambda_{0}^{2}");
 722           fhELambda0MCHadron[iso]->SetXTitle("E (GeV)");
 723           outputContainer->Add(fhELambda0MCHadron[iso]) ;
 724         }
 725
 726         fhELambda1[iso]  = new TH2F
 727         (Form("hELambda1%s",hName[iso].Data()),
 728          Form("%s cluster: E vs #lambda_{1}",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 729         fhELambda1[iso]->SetYTitle("#lambda_{1}^{2}");
 730         fhELambda1[iso]->SetXTitle("E (GeV)");
 731         outputContainer->Add(fhELambda1[iso]) ;
 732
 733         if(fCalorimeter=="EMCAL")
 734         {
 735           fhELambda0TRD[iso]  = new TH2F
 736           (Form("hELambda0TRD%s",hName[iso].Data()),
 737            Form("%s cluster: E vs #lambda_{0}, SM behind TRD",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 738           fhELambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
 739           fhELambda0TRD[iso]->SetXTitle("E (GeV)");
 740           outputContainer->Add(fhELambda0TRD[iso]) ;
 741
 742           fhELambda1TRD[iso]  = new TH2F
 743           (Form("hELambda1TRD%s",hName[iso].Data()),
 744            Form("%s cluster: E vs #lambda_{1}, SM behind TRD",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 745           fhELambda1TRD[iso]->SetYTitle("#lambda_{1}^{2}");
 746           fhELambda1TRD[iso]->SetXTitle("E (GeV)");
 747           outputContainer->Add(fhELambda1TRD[iso]) ;
 748         }
 749
 750         fhNLocMax[iso] = new TH2F
 751         (Form("hNLocMax%s",hName[iso].Data()),
 752          Form("%s - Number of local maxima in cluster",hTitle[iso].Data()),
 753          nptbins,ptmin,ptmax,10,0,10);
 754         fhNLocMax[iso]->SetYTitle("N maxima");
 755         fhNLocMax[iso]->SetXTitle("E (GeV)");
 756         outputContainer->Add(fhNLocMax[iso]) ;
 757
 758         fhELambda0LocMax1[iso]  = new TH2F
 759         (Form("hELambda0LocMax1%s",hName[iso].Data()),
 760          Form("%s cluster (#eta) pairs: E vs #lambda_{0}, 1 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 761         fhELambda0LocMax1[iso]->SetYTitle("#lambda_{0}^{2}");
 762         fhELambda0LocMax1[iso]->SetXTitle("E (GeV)");
 763         outputContainer->Add(fhELambda0LocMax1[iso]) ;
 764
 765         fhELambda1LocMax1[iso]  = new TH2F
 766         (Form("hELambda1LocMax1%s",hName[iso].Data()),
 767          Form("%s cluster (#eta) pairs: E vs #lambda_{1}, 1 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 768         fhELambda1LocMax1[iso]->SetYTitle("#lambda_{1}^{2}");
 769         fhELambda1LocMax1[iso]->SetXTitle("E (GeV)");
 770         outputContainer->Add(fhELambda1LocMax1[iso]) ;
 771
 772         fhELambda0LocMax2[iso]  = new TH2F
 773         (Form("hELambda0LocMax2%s",hName[iso].Data()),
 774          Form("%s cluster (#eta) pairs: E vs #lambda_{0}, 2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 775         fhELambda0LocMax2[iso]->SetYTitle("#lambda_{0}^{2}");
 776         fhELambda0LocMax2[iso]->SetXTitle("E (GeV)");
 777         outputContainer->Add(fhELambda0LocMax2[iso]) ;
 778
 779         fhELambda1LocMax2[iso]  = new TH2F
 780         (Form("hELambda1LocMax2%s",hName[iso].Data()),
 781          Form("%s cluster (#eta) pairs: E vs #lambda_{1}, 2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 782         fhELambda1LocMax2[iso]->SetYTitle("#lambda_{1}^{2}");
 783         fhELambda1LocMax2[iso]->SetXTitle("E (GeV)");
 784         outputContainer->Add(fhELambda1LocMax2[iso]) ;
 785
 786         fhELambda0LocMaxN[iso]  = new TH2F
 787         ( Form("hELambda0LocMaxN%s",hName[iso].Data()),
 788          Form("%s cluster (#eta) pairs: E vs #lambda_{0}, N>2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 789         fhELambda0LocMaxN[iso]->SetYTitle("#lambda_{0}^{2}");
 790         fhELambda0LocMaxN[iso]->SetXTitle("E (GeV)");
 791         outputContainer->Add(fhELambda0LocMaxN[iso]) ;
 792
 793         fhELambda1LocMaxN[iso]  = new TH2F
 794         (Form("hELambda1LocMaxN%s",hName[iso].Data()),
 795          Form("%s cluster (#eta) pairs: E vs #lambda_{1}, N>2 Local maxima",hTitle[iso].Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 796         fhELambda1LocMaxN[iso]->SetYTitle("#lambda_{1}^{2}");
 797         fhELambda1LocMaxN[iso]->SetXTitle("E (GeV)");
 798         outputContainer->Add(fhELambda1LocMaxN[iso]) ;
 799
 800       }
 801     } // control histograms for isolated and non isolated objects
 802
 803     fhConeSumPt  = new TH2F("hConePtSum",
 804                             Form("#Sigma p_{T} in isolation cone for R = %2.2f",r),
 805                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 806     fhConeSumPt->SetYTitle("#Sigma p_{T}");
 807     fhConeSumPt->SetXTitle("p_{T} (GeV/c)");
 808     outputContainer->Add(fhConeSumPt) ;
 809
 810     fhPtInCone  = new TH2F("hPtInCone",
 811                            Form("p_{T} of clusters and tracks in isolation cone for R = %2.2f",r),
 812                            nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 813     fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
 814     fhPtInCone->SetXTitle("p_{T} (GeV/c)");
 815     outputContainer->Add(fhPtInCone) ;
 816
 817     fhPtTrackInCone  = new TH2F("hPtTrackInCone",
 818                            Form("p_{T} of tracks in isolation cone for R = %2.2f",r),
 819                            nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 820     fhPtTrackInCone->SetYTitle("p_{T in cone} (GeV/c)");
 821     fhPtTrackInCone->SetXTitle("p_{T} (GeV/c)");
 822     outputContainer->Add(fhPtTrackInCone) ;
 823
 824     if(fFillPileUpHistograms)
 825     {
 826       fhPtTrackInConeOtherBC  = new TH2F("hPtTrackInConeOtherBC",
 827                                          Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC!=0",r),
 828                                          nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 829       fhPtTrackInConeOtherBC->SetYTitle("p_{T in cone} (GeV/c)");
 830       fhPtTrackInConeOtherBC->SetXTitle("p_{T} (GeV/c)");
 831       outputContainer->Add(fhPtTrackInConeOtherBC) ;
 832
 833       fhPtTrackInConeOtherBCPileUpSPD  = new TH2F("hPtTrackInConeOtherBCPileUpSPD",
 834                                          Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC!=0, pile-up from SPD",r),
 835                                          nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 836       fhPtTrackInConeOtherBCPileUpSPD->SetYTitle("p_{T in cone} (GeV/c)");
 837       fhPtTrackInConeOtherBCPileUpSPD->SetXTitle("p_{T} (GeV/c)");
 838       outputContainer->Add(fhPtTrackInConeOtherBCPileUpSPD) ;
 839
 840       fhPtTrackInConeBC0  = new TH2F("hPtTrackInConeBC0",
 841                                          Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC==0",r),
 842                                          nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 843       fhPtTrackInConeBC0->SetYTitle("p_{T in cone} (GeV/c)");
 844       fhPtTrackInConeBC0->SetXTitle("p_{T} (GeV/c)");
 845       outputContainer->Add(fhPtTrackInConeBC0) ;
 846
 847       fhPtTrackInConeVtxBC0  = new TH2F("hPtTrackInConeVtxBC0",
 848                                      Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC==0",r),
 849                                      nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 850       fhPtTrackInConeVtxBC0->SetYTitle("p_{T in cone} (GeV/c)");
 851       fhPtTrackInConeVtxBC0->SetXTitle("p_{T} (GeV/c)");
 852       outputContainer->Add(fhPtTrackInConeVtxBC0) ;
 853
 854
 855       fhPtTrackInConeBC0PileUpSPD  = new TH2F("hPtTrackInConeBC0PileUpSPD",
 856                                                   Form("p_{T} of tracks in isolation cone for R = %2.2f, TOF from BC==0, pile-up from SPD",r),
 857                                                   nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 858       fhPtTrackInConeBC0PileUpSPD->SetYTitle("p_{T in cone} (GeV/c)");
 859       fhPtTrackInConeBC0PileUpSPD->SetXTitle("p_{T} (GeV/c)");
 860       outputContainer->Add(fhPtTrackInConeBC0PileUpSPD) ;
 861
 862
 863       for (Int_t i = 0; i < 7 ; i++)
 864       {
 865         fhPtInConePileUp[i]  = new TH2F(Form("hPtInConePileUp%s",pileUpName[i].Data()),
 866                                         Form("p_{T} in isolation cone for R = %2.2f, from pile-up (%s)",r,pileUpName[i].Data()),
 867                                         nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 868         fhPtInConePileUp[i]->SetYTitle("p_{T in cone} (GeV/c)");
 869         fhPtInConePileUp[i]->SetXTitle("p_{T} (GeV/c)");
 870         outputContainer->Add(fhPtInConePileUp[i]) ;
 871       }
 872     }
 873
 874     fhPtInConeCent  = new TH2F("hPtInConeCent",
 875                                Form("p_{T} in isolation cone for R = %2.2f",r),
 876                                100,0,100,nptinconebins,ptinconemin,ptinconemax);
 877     fhPtInConeCent->SetYTitle("p_{T in cone} (GeV/c)");
 878     fhPtInConeCent->SetXTitle("centrality");
 879     outputContainer->Add(fhPtInConeCent) ;
 880
 881     fhFRConeSumPt  = new TH2F("hFRConePtSum",
 882                               Form("#Sigma p_{T} in the forward region isolation cone for R = %2.2f",r),
 883                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 884     fhFRConeSumPt->SetYTitle("#Sigma p_{T}");
 885     fhFRConeSumPt->SetXTitle("p_{T} (GeV/c)");
 886     outputContainer->Add(fhFRConeSumPt) ;
 887
 888     fhPtInFRCone  = new TH2F("hPtInFRCone",
 889                              Form("p_{T} in forward region isolation cone for R = %2.2f",r),
 890                              nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 891     fhPtInFRCone->SetYTitle("p_{T in cone} (GeV/c)");
 892     fhPtInFRCone->SetXTitle("p_{T} (GeV/c)");
 893     outputContainer->Add(fhPtInFRCone) ;
 894
 895     fhPhiUEConeSumPt  = new TH2F("hPhiUEConeSumPt",
 896                                  Form("p_{T} in phi band around isolation cone for R = %2.2f",r),
 897                                  nptbins,ptmin,ptmax,4*nptinconebins,ptinconemin,4*ptinconemax);
 898     fhPhiUEConeSumPt->SetYTitle("p_{T in band} (GeV/c)");
 899     fhPhiUEConeSumPt->SetXTitle("p_{T} (GeV/c)");
 900     outputContainer->Add(fhPhiUEConeSumPt) ;
 901
 902     fhEtaUEConeSumPt  = new TH2F("hEtaUEConeSumPt",
 903                                  Form("p_{T} in eta band around isolation cone for R = %2.2f",r),
 904                                  nptbins,ptmin,ptmax,4*nptinconebins,ptinconemin,4*ptinconemax);
 905     fhEtaUEConeSumPt->SetYTitle("p_{T in band} (GeV/c)");
 906     fhEtaUEConeSumPt->SetXTitle("p_{T} (GeV/c)");
 907     outputContainer->Add(fhEtaUEConeSumPt) ;
 908
 909     fhEtaBand  = new TH2F("fhEtaBand",
 910                           Form("Eta/Phi of particle in Eta band isolation cone for R = %2.2f",r),
 911                           netabins,etamin,etamax,nphibins,phimin,phimax);
 912     fhEtaBand->SetXTitle("#eta");
 913     fhEtaBand->SetYTitle("#phi");
 914     outputContainer->Add(fhEtaBand) ;
 915
 916     fhPhiBand  = new TH2F("fhPhiBand",
 917                           Form("Eta/Phi of particle in Phi band isolation cone for R = %2.2f",r),
 918                           netabins,etamin,etamax,nphibins,phimin,phimax);
 919     fhPhiBand->SetXTitle("#eta");
 920     fhPhiBand->SetYTitle("#phi");
 921     outputContainer->Add(fhPhiBand) ;
 922
 923     fhConeSumPtEtaUESub  = new TH2F("hConeSumPtEtaUESub",
 924                                     Form("#Sigma p_{T} after bkg subtraction from eta band in the isolation cone for R = %2.2f",r),
 925                                     nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
 926     fhConeSumPtEtaUESub->SetYTitle("#Sigma p_{T}");
 927     fhConeSumPtEtaUESub->SetXTitle("p_{T} (GeV/c)");
 928     outputContainer->Add(fhConeSumPtEtaUESub) ;
 929
 930     fhConeSumPtPhiUESub  = new TH2F("hConeSumPtPhiUESub",
 931                                     Form("#Sigma p_{T} after bkg subtraction from phi band in the isolation cone for R = %2.2f",r),
 932                                     nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
 933     fhConeSumPtPhiUESub->SetYTitle("#Sigma p_{T}");
 934     fhConeSumPtPhiUESub->SetXTitle("p_{T} (GeV/c)");
 935     outputContainer->Add(fhConeSumPtPhiUESub) ;
 936
 937     fhEIso   = new TH1F("hE",
 938                         Form("Number of isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 939                         nptbins,ptmin,ptmax);
 940     fhEIso->SetYTitle("dN / dE");
 941     fhEIso->SetXTitle("E (GeV/c)");
 942     outputContainer->Add(fhEIso) ;
 943
 944     fhPtIso  = new TH1F("hPt",
 945                         Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 946                         nptbins,ptmin,ptmax);
 947     fhPtIso->SetYTitle("dN / p_{T}");
 948     fhPtIso->SetXTitle("p_{T} (GeV/c)");
 949     outputContainer->Add(fhPtIso) ;
 950
 951     fhPtCentralityIso  = new TH2F("hPtCentrality","centrality vs p_{T} for isolated particles",nptbins,ptmin,ptmax, 100,0,100);
 952     fhPtCentralityIso->SetYTitle("centrality");
 953     fhPtCentralityIso->SetXTitle("p_{T}(GeV/c)");
 954     outputContainer->Add(fhPtCentralityIso) ;
 955
 956     fhPtEventPlaneIso  = new TH2F("hPtEventPlane","event plane angle vs p_{T} for isolated particles",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
 957     fhPtEventPlaneIso->SetYTitle("Event plane angle (rad)");
 958     fhPtEventPlaneIso->SetXTitle("p_{T} (GeV/c)");
 959     outputContainer->Add(fhPtEventPlaneIso) ;
 960
 961
 962     fhPtNLocMaxIso  = new TH2F("hPtNLocMax",
 963                                Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f vs NLM, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 964                                nptbins,ptmin,ptmax,10,0,10);
 965     fhPtNLocMaxIso->SetYTitle("NLM");
 966     fhPtNLocMaxIso->SetXTitle("p_{T} (GeV/c)");
 967     outputContainer->Add(fhPtNLocMaxIso) ;
 968
 969     fhPhiIso  = new TH2F("hPhi",
 970                          Form("Number of isolated particles vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 971                          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 972     fhPhiIso->SetYTitle("#phi");
 973     fhPhiIso->SetXTitle("p_{T} (GeV/c)");
 974     outputContainer->Add(fhPhiIso) ;
 975
 976     fhEtaIso  = new TH2F("hEta",
 977                          Form("Number of isolated particles vs #eta for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 978                          nptbins,ptmin,ptmax,netabins,etamin,etamax);
 979     fhEtaIso->SetYTitle("#eta");
 980     fhEtaIso->SetXTitle("p_{T} (GeV/c)");
 981     outputContainer->Add(fhEtaIso) ;
 982
 983     fhEtaPhiIso  = new TH2F("hEtaPhiIso",
 984                             Form("Number of isolated particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 985                             netabins,etamin,etamax,nphibins,phimin,phimax);
 986     fhEtaPhiIso->SetXTitle("#eta");
 987     fhEtaPhiIso->SetYTitle("#phi");
 988     outputContainer->Add(fhEtaPhiIso) ;
 989
 990     fhPtDecayIso  = new TH1F("hPtDecayIso",
 991                              Form("Number of isolated #pi^{0} decay particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 992                              nptbins,ptmin,ptmax);
 993     fhPtDecayIso->SetYTitle("N");
 994     fhPtDecayIso->SetXTitle("p_{T}(GeV/c)");
 995     outputContainer->Add(fhPtDecayIso) ;
 996
 997     fhEtaPhiDecayIso  = new TH2F("hEtaPhiDecayIso",
 998                                  Form("Number of isolated Pi0 decay particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 999                                  netabins,etamin,etamax,nphibins,phimin,phimax);
1000     fhEtaPhiDecayIso->SetXTitle("#eta");
1001     fhEtaPhiDecayIso->SetYTitle("#phi");
1002     outputContainer->Add(fhEtaPhiDecayIso) ;
1003
1004     if(IsDataMC())
1005     {
1006       fhPtIsoPrompt  = new TH1F("hPtMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax);
1007       fhPtIsoPrompt->SetYTitle("N");
1008       fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
1009       outputContainer->Add(fhPtIsoPrompt) ;
1010
1011       fhPhiIsoPrompt  = new TH2F
1012       ("hPhiMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1013       fhPhiIsoPrompt->SetYTitle("#phi");
1014       fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
1015       outputContainer->Add(fhPhiIsoPrompt) ;
1016
1017       fhEtaIsoPrompt  = new TH2F
1018       ("hEtaMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1019       fhEtaIsoPrompt->SetYTitle("#eta");
1020       fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
1021       outputContainer->Add(fhEtaIsoPrompt) ;
1022
1023       fhPtIsoFragmentation  = new TH1F("hPtMCFragmentation","Number of isolated #gamma",nptbins,ptmin,ptmax);
1024       fhPtIsoFragmentation->SetYTitle("N");
1025       fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
1026       outputContainer->Add(fhPtIsoFragmentation) ;
1027
1028       fhPhiIsoFragmentation  = new TH2F
1029       ("hPhiMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1030       fhPhiIsoFragmentation->SetYTitle("#phi");
1031       fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
1032       outputContainer->Add(fhPhiIsoFragmentation) ;
1033
1034       fhEtaIsoFragmentation  = new TH2F
1035       ("hEtaMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1036       fhEtaIsoFragmentation->SetYTitle("#eta");
1037       fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
1038       outputContainer->Add(fhEtaIsoFragmentation) ;
1039
1040       fhPtIsoPi0  = new TH1F("hPtMCPi0","Number of isolated #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax);
1041       fhPtIsoPi0->SetYTitle("N");
1042       fhPtIsoPi0->SetXTitle("p_{T #gamma}(GeV/c)");
1043       outputContainer->Add(fhPtIsoPi0) ;
1044
1045       fhPhiIsoPi0  = new TH2F
1046       ("hPhiMCPi0","Number of isolated #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1047       fhPhiIsoPi0->SetYTitle("#phi");
1048       fhPhiIsoPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1049       outputContainer->Add(fhPhiIsoPi0) ;
1050
1051       fhEtaIsoPi0  = new TH2F
1052       ("hEtaMCPi0","Number of isolated #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1053       fhEtaIsoPi0->SetYTitle("#eta");
1054       fhEtaIsoPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1055       outputContainer->Add(fhEtaIsoPi0) ;
1056
1057       fhPtIsoPi0Decay  = new TH1F("hPtMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
1058       fhPtIsoPi0Decay->SetYTitle("N");
1059       fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
1060       outputContainer->Add(fhPtIsoPi0Decay) ;
1061
1062       fhPhiIsoPi0Decay  = new TH2F
1063       ("hPhiMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1064       fhPhiIsoPi0Decay->SetYTitle("#phi");
1065       fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
1066       outputContainer->Add(fhPhiIsoPi0Decay) ;
1067
1068       fhEtaIsoPi0Decay  = new TH2F
1069       ("hEtaMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1070       fhEtaIsoPi0Decay->SetYTitle("#eta");
1071       fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
1072       outputContainer->Add(fhEtaIsoPi0Decay) ;
1073
1074       fhPtIsoEtaDecay  = new TH1F("hPtMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax);
1075       fhPtIsoEtaDecay->SetYTitle("N");
1076       fhPtIsoEtaDecay->SetXTitle("p_{T #gamma}(GeV/c)");
1077       outputContainer->Add(fhPtIsoEtaDecay) ;
1078
1079       fhPhiIsoEtaDecay  = new TH2F
1080       ("hPhiMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1081       fhPhiIsoEtaDecay->SetYTitle("#phi");
1082       fhPhiIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
1083       outputContainer->Add(fhPhiIsoEtaDecay) ;
1084
1085       fhEtaIsoEtaDecay  = new TH2F
1086       ("hEtaMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1087       fhEtaIsoEtaDecay->SetYTitle("#eta");
1088       fhEtaIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
1089       outputContainer->Add(fhEtaIsoEtaDecay) ;
1090
1091       fhPtIsoOtherDecay  = new TH1F("hPtMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax);
1092       fhPtIsoOtherDecay->SetYTitle("N");
1093       fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
1094       outputContainer->Add(fhPtIsoOtherDecay) ;
1095
1096       fhPhiIsoOtherDecay  = new TH2F
1097       ("hPhiMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1098       fhPhiIsoOtherDecay->SetYTitle("#phi");
1099       fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
1100       outputContainer->Add(fhPhiIsoOtherDecay) ;
1101
1102       fhEtaIsoOtherDecay  = new TH2F
1103       ("hEtaMCOtherDecay","Number of isolated #gamma non-#pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1104       fhEtaIsoOtherDecay->SetYTitle("#eta");
1105       fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
1106       outputContainer->Add(fhEtaIsoOtherDecay) ;
1107
1108       //      fhPtIsoConversion  = new TH1F("hPtMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax);
1109       //      fhPtIsoConversion->SetYTitle("N");
1110       //      fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
1111       //      outputContainer->Add(fhPtIsoConversion) ;
1112       //
1113       //      fhPhiIsoConversion  = new TH2F
1114       //      ("hPhiMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1115       //      fhPhiIsoConversion->SetYTitle("#phi");
1116       //      fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
1117       //      outputContainer->Add(fhPhiIsoConversion) ;
1118       //
1119       //      fhEtaIsoConversion  = new TH2F
1120       //      ("hEtaMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1121       //      fhEtaIsoConversion->SetYTitle("#eta");
1122       //      fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
1123       //      outputContainer->Add(fhEtaIsoConversion) ;
1124
1125       fhPtIsoHadron  = new TH1F("hPtMCHadron","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax);
1126       fhPtIsoHadron->SetYTitle("N");
1127       fhPtIsoHadron->SetXTitle("p_{T}(GeV/c)");
1128       outputContainer->Add(fhPtIsoHadron) ;
1129
1130
1131       fhPhiIsoHadron  = new TH2F
1132       ("hPhiMCHadron","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1133       fhPhiIsoHadron->SetYTitle("#phi");
1134       fhPhiIsoHadron->SetXTitle("p_{T} (GeV/c)");
1135       outputContainer->Add(fhPhiIsoHadron) ;
1136
1137       fhEtaIsoHadron  = new TH2F
1138       ("hEtaMCHadron","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1139       fhEtaIsoHadron->SetYTitle("#eta");
1140       fhEtaIsoHadron->SetXTitle("p_{T} (GeV/c)");
1141       outputContainer->Add(fhEtaIsoHadron) ;
1142
1143     }//Histos with MC
1144
1145   }
1146
1147   // Not Isolated histograms, reference histograms
1148
1149   fhENoIso  = new TH1F("hENoIso",
1150                         Form("Number of not isolated leading particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1151                         nptbins,ptmin,ptmax);
1152   fhENoIso->SetYTitle("N");
1153   fhENoIso->SetXTitle("p_{T}(GeV/c)");
1154   outputContainer->Add(fhENoIso) ;
1155
1156   fhPtNoIso  = new TH1F("hPtNoIso",
1157                         Form("Number of not isolated leading particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1158                         nptbins,ptmin,ptmax);
1159   fhPtNoIso->SetYTitle("N");
1160   fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
1161   outputContainer->Add(fhPtNoIso) ;
1162
1163   fhPtNLocMaxNoIso  = new TH2F("hPtNLocMaxNoIso",
1164                                Form("Number of not isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f vs NLM, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1165                                nptbins,ptmin,ptmax,10,0,10);
1166   fhPtNLocMaxNoIso->SetYTitle("NLM");
1167   fhPtNLocMaxNoIso->SetXTitle("p_{T} (GeV/c)");
1168   outputContainer->Add(fhPtNLocMaxNoIso) ;
1169
1170   fhEtaPhiNoIso  = new TH2F("hEtaPhiNoIso",
1171                             Form("Number of not isolated leading particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1172                             netabins,etamin,etamax,nphibins,phimin,phimax);
1173   fhEtaPhiNoIso->SetXTitle("#eta");
1174   fhEtaPhiNoIso->SetYTitle("#phi");
1175   outputContainer->Add(fhEtaPhiNoIso) ;
1176
1177   fhPtDecayNoIso  = new TH1F("hPtDecayNoIso",
1178                              Form("Number of not isolated leading pi0 decay particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1179                              nptbins,ptmin,ptmax);
1180   fhPtDecayNoIso->SetYTitle("N");
1181   fhPtDecayNoIso->SetXTitle("p_{T}(GeV/c)");
1182   outputContainer->Add(fhPtDecayNoIso) ;
1183
1184   fhEtaPhiDecayNoIso  = new TH2F("hEtaPhiDecayNoIso",
1185                                  Form("Number of not isolated leading Pi0 decay particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
1186                                  netabins,etamin,etamax,nphibins,phimin,phimax);
1187   fhEtaPhiDecayNoIso->SetXTitle("#eta");
1188   fhEtaPhiDecayNoIso->SetYTitle("#phi");
1189   outputContainer->Add(fhEtaPhiDecayNoIso) ;
1190
1191
1192
1193   if(IsDataMC())
1194   {
1195     fhPtNoIsoPi0  = new TH1F
1196     ("hPtNoIsoPi0","Number of not isolated leading #gamma from #pi^{0} (2 #gamma)",nptbins,ptmin,ptmax);
1197     fhPtNoIsoPi0->SetYTitle("N");
1198     fhPtNoIsoPi0->SetXTitle("p_{T} (GeV/c)");
1199     outputContainer->Add(fhPtNoIsoPi0) ;
1200
1201     fhPtNoIsoPi0Decay  = new TH1F
1202     ("hPtNoIsoPi0Decay","Number of not isolated leading #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
1203     fhPtNoIsoPi0Decay->SetYTitle("N");
1204     fhPtNoIsoPi0Decay->SetXTitle("p_{T} (GeV/c)");
1205     outputContainer->Add(fhPtNoIsoPi0Decay) ;
1206
1207     fhPtNoIsoEtaDecay  = new TH1F
1208     ("hPtNoIsoEtaDecay","Number of not isolated leading #gamma from eta decay",nptbins,ptmin,ptmax);
1209     fhPtNoIsoEtaDecay->SetYTitle("N");
1210     fhPtNoIsoEtaDecay->SetXTitle("p_{T} (GeV/c)");
1211     outputContainer->Add(fhPtNoIsoEtaDecay) ;
1212
1213     fhPtNoIsoOtherDecay  = new TH1F
1214     ("hPtNoIsoOtherDecay","Number of not isolated leading #gamma from other decay",nptbins,ptmin,ptmax);
1215     fhPtNoIsoOtherDecay->SetYTitle("N");
1216     fhPtNoIsoOtherDecay->SetXTitle("p_{T} (GeV/c)");
1217     outputContainer->Add(fhPtNoIsoOtherDecay) ;
1218
1219     fhPtNoIsoPrompt  = new TH1F
1220     ("hPtNoIsoPrompt","Number of not isolated leading prompt #gamma",nptbins,ptmin,ptmax);
1221     fhPtNoIsoPrompt->SetYTitle("N");
1222     fhPtNoIsoPrompt->SetXTitle("p_{T} (GeV/c)");
1223     outputContainer->Add(fhPtNoIsoPrompt) ;
1224
1225     fhPtIsoMCPhoton  = new TH1F
1226     ("hPtIsoMCPhoton","Number of isolated leading  #gamma",nptbins,ptmin,ptmax);
1227     fhPtIsoMCPhoton->SetYTitle("N");
1228     fhPtIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
1229     outputContainer->Add(fhPtIsoMCPhoton) ;
1230
1231     fhPtNoIsoMCPhoton  = new TH1F
1232     ("hPtNoIsoMCPhoton","Number of not isolated leading #gamma",nptbins,ptmin,ptmax);
1233     fhPtNoIsoMCPhoton->SetYTitle("N");
1234     fhPtNoIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
1235     outputContainer->Add(fhPtNoIsoMCPhoton) ;
1236
1237 //    fhPtNoIsoConversion  = new TH1F
1238 //    ("hPtNoIsoConversion","Number of not isolated leading conversion #gamma",nptbins,ptmin,ptmax);
1239 //    fhPtNoIsoConversion->SetYTitle("N");
1240 //    fhPtNoIsoConversion->SetXTitle("p_{T} (GeV/c)");
1241 //    outputContainer->Add(fhPtNoIsoConversion) ;
1242
1243     fhPtNoIsoFragmentation  = new TH1F
1244     ("hPtNoIsoFragmentation","Number of not isolated leading fragmentation #gamma",nptbins,ptmin,ptmax);
1245     fhPtNoIsoFragmentation->SetYTitle("N");
1246     fhPtNoIsoFragmentation->SetXTitle("p_{T} (GeV/c)");
1247     outputContainer->Add(fhPtNoIsoFragmentation) ;
1248
1249     fhPtNoIsoHadron  = new TH1F
1250     ("hPtNoIsoHadron","Number of not isolated leading hadrons",nptbins,ptmin,ptmax);
1251     fhPtNoIsoHadron->SetYTitle("N");
1252     fhPtNoIsoHadron->SetXTitle("p_{T} (GeV/c)");
1253     outputContainer->Add(fhPtNoIsoHadron) ;
1254
1255   }//Histos with MC
1256
1257
1258   if(fMakeSeveralIC)
1259   {
1260     const Int_t buffersize = 255;
1261     char name[buffersize];
1262     char title[buffersize];
1263     for(Int_t icone = 0; icone<fNCones; icone++)
1264     {
1265       // sum pt in cone vs. pt leading
1266       snprintf(name, buffersize,"hSumPtLeadingPt_Cone_%d",icone);
1267       snprintf(title, buffersize,"#Sigma p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
1268       fhSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1269       fhSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}p_{T} (GeV/c)");//#Sigma p_{T}
1270       fhSumPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
1271       outputContainer->Add(fhSumPtLeadingPt[icone]) ;
1272
1273       // pt in cone vs. pt leading
1274       snprintf(name, buffersize,"hPtLeadingPt_Cone_%d",icone);
1275       snprintf(title, buffersize,"p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
1276       fhPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1277       fhPtLeadingPt[icone] ->SetYTitle("p_{T}^{cone} (GeV/c)");
1278       fhPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
1279       outputContainer->Add(fhPtLeadingPt[icone]) ;
1280
1281        // sum pt in cone vs. pt leading in the forward region (for background subtraction studies)
1282         snprintf(name, buffersize,"hFRSumPtLeadingPt_Cone_%d",icone);
1283       snprintf(title, buffersize,"#Sigma p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
1284       fhFRSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1285       fhFRSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}p_{T} (GeV/c)");//#Sigma p_{T}
1286       fhFRSumPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
1287       outputContainer->Add(fhFRSumPtLeadingPt[icone]) ;
1288
1289       // pt in cone vs. pt leading in the forward region (for background subtraction studies)
1290         snprintf(name, buffersize,"hFRPtLeadingPt_Cone_%d",icone);
1291       snprintf(title, buffersize,"p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
1292       fhFRPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1293       fhFRPtLeadingPt[icone] ->SetYTitle("p_{T}^{cone} (GeV/c)");
1294       fhFRPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
1295       outputContainer->Add(fhFRPtLeadingPt[icone]) ;
1296
1297
1298       if(IsDataMC())
1299       {
1300         snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);
1301         snprintf(title, buffersize,"Candidate Prompt cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1302         fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1303         fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1304         fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
1305         outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
1306
1307         snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);
1308         snprintf(title, buffersize,"Candidate Fragmentation cone sum p_{T} for R = %2.2fvs candidate p_{T}",fConeSizes[icone]);
1309         fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1310         fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1311         fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
1312         outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
1313
1314         snprintf(name, buffersize,"hPtSumPi0_Cone_%d",icone);
1315         snprintf(title, buffersize,"Candidate Pi0 cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1316         fhPtSumIsolatedPi0[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1317         fhPtSumIsolatedPi0[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1318         fhPtSumIsolatedPi0[icone]->SetXTitle("p_{T} (GeV/c)");
1319         outputContainer->Add(fhPtSumIsolatedPi0[icone]) ;
1320
1321         snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);
1322         snprintf(title, buffersize,"Candidate Pi0Decay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1323         fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1324         fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1325         fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
1326         outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
1327
1328         snprintf(name, buffersize,"hPtSumEtaDecay_Cone_%d",icone);
1329         snprintf(title, buffersize,"Candidate EtaDecay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1330         fhPtSumIsolatedEtaDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1331         fhPtSumIsolatedEtaDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1332         fhPtSumIsolatedEtaDecay[icone]->SetXTitle("p_{T} (GeV/c)");
1333         outputContainer->Add(fhPtSumIsolatedEtaDecay[icone]) ;
1334
1335         snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);
1336         snprintf(title, buffersize,"Candidate OtherDecay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1337         fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1338         fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1339         fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
1340         outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
1341
1342 //        snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);
1343 //        snprintf(title, buffersize,"Candidate Conversion cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1344 //        fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1345 //        fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1346 //        fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
1347 //        outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
1348
1349         snprintf(name, buffersize,"hPtSumHadron_Cone_%d",icone);
1350         snprintf(title, buffersize,"Candidate Hadron cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
1351         fhPtSumIsolatedHadron[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1352         fhPtSumIsolatedHadron[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1353         fhPtSumIsolatedHadron[icone]->SetXTitle("p_{T} (GeV/c)");
1354         outputContainer->Add(fhPtSumIsolatedHadron[icone]) ;
1355
1356       }//Histos with MC
1357
1358       for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++)
1359       {
1360
1361         snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
1362         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
1363         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1364         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1365         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
1366
1367         snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
1368         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1369         fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1370         fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1371         outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
1372
1373
1374         snprintf(name, buffersize,"hPtSum_Cone_%d_Pt%d",icone,ipt);
1375         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1376         fhPtSumIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1377         // fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
1378         fhPtSumIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1379         outputContainer->Add(fhPtSumIsolated[icone][ipt]) ;
1380
1381         snprintf(name, buffersize,"hPtSumDensity_Cone_%d_Pt%d",icone,ipt);
1382         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for density in R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1383         fhPtSumDensityIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
1384         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
1385         fhPtSumDensityIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1386         outputContainer->Add(fhPtSumDensityIso[icone][ipt]) ;
1387
1388         snprintf(name, buffersize,"hPtFracPtSum_Cone_%d_Pt%d",icone,ipt);
1389         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for PtFracPtSum in R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1390         fhPtFracPtSumIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
1391         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
1392         fhPtFracPtSumIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1393         outputContainer->Add(fhPtFracPtSumIso[icone][ipt]) ;
1394
1395         // pt decays isolated
1396         snprintf(name, buffersize,"hPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
1397         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
1398         fhPtPtThresDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1399         fhPtPtThresDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1400         outputContainer->Add(fhPtPtThresDecayIso[icone][ipt]) ;
1401
1402         snprintf(name, buffersize,"hPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
1403         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1404         fhPtPtFracDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1405         fhPtPtFracDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1406         outputContainer->Add(fhPtPtFracDecayIso[icone][ipt]) ;
1407
1408         snprintf(name, buffersize,"hPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
1409         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1410         fhPtPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
1411         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1412         fhPtPtSumDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1413         outputContainer->Add(fhPtPtSumDecayIso[icone][ipt]) ;
1414
1415         snprintf(name, buffersize,"hPtSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
1416         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for density in R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1417         fhPtSumDensityDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
1418         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1419         fhPtSumDensityDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1420         outputContainer->Add(fhPtSumDensityDecayIso[icone][ipt]) ;
1421
1422         snprintf(name, buffersize,"hPtFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
1423         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for PtFracPtSum in R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1424         fhPtFracPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
1425         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
1426         fhPtFracPtSumDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1427         outputContainer->Add(fhPtFracPtSumDecayIso[icone][ipt]) ;
1428
1429
1430         // eta:phi
1431         snprintf(name, buffersize,"hEtaPhiPtThres_Cone_%d_Pt%d",icone,ipt);
1432         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
1433         fhEtaPhiPtThresIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1434         fhEtaPhiPtThresIso[icone][ipt]->SetXTitle("#eta");
1435         fhEtaPhiPtThresIso[icone][ipt]->SetYTitle("#phi");
1436         outputContainer->Add(fhEtaPhiPtThresIso[icone][ipt]) ;
1437
1438         snprintf(name, buffersize,"hEtaPhiPtFrac_Cone_%d_Pt%d",icone,ipt);
1439         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1440         fhEtaPhiPtFracIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1441         fhEtaPhiPtFracIso[icone][ipt]->SetXTitle("#eta");
1442         fhEtaPhiPtFracIso[icone][ipt]->SetYTitle("#phi");
1443         outputContainer->Add(fhEtaPhiPtFracIso[icone][ipt]) ;
1444
1445         snprintf(name, buffersize,"hEtaPhiPtSum_Cone_%d_Pt%d",icone,ipt);
1446         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1447         fhEtaPhiPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1448         fhEtaPhiPtSumIso[icone][ipt]->SetXTitle("#eta");
1449         fhEtaPhiPtSumIso[icone][ipt]->SetYTitle("#phi");
1450         outputContainer->Add(fhEtaPhiPtSumIso[icone][ipt]) ;
1451
1452         snprintf(name, buffersize,"hEtaPhiSumDensity_Cone_%d_Pt%d",icone,ipt);
1453         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for density R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1454         fhEtaPhiSumDensityIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1455         fhEtaPhiSumDensityIso[icone][ipt]->SetXTitle("#eta");
1456         fhEtaPhiSumDensityIso[icone][ipt]->SetYTitle("#phi");
1457         outputContainer->Add(fhEtaPhiSumDensityIso[icone][ipt]) ;
1458
1459         snprintf(name, buffersize,"hEtaPhiFracPtSum_Cone_%d_Pt%d",icone,ipt);
1460         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for FracPtSum R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1461         fhEtaPhiFracPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1462         fhEtaPhiFracPtSumIso[icone][ipt]->SetXTitle("#eta");
1463         fhEtaPhiFracPtSumIso[icone][ipt]->SetYTitle("#phi");
1464         outputContainer->Add(fhEtaPhiFracPtSumIso[icone][ipt]) ;
1465
1466         // eta:phi decays
1467         snprintf(name, buffersize,"hEtaPhiPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
1468         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
1469         fhEtaPhiPtThresDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1470         fhEtaPhiPtThresDecayIso[icone][ipt]->SetXTitle("#eta");
1471         fhEtaPhiPtThresDecayIso[icone][ipt]->SetYTitle("#phi");
1472         outputContainer->Add(fhEtaPhiPtThresDecayIso[icone][ipt]) ;
1473
1474         snprintf(name, buffersize,"hEtaPhiPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
1475         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1476         fhEtaPhiPtFracDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1477         fhEtaPhiPtFracDecayIso[icone][ipt]->SetXTitle("#eta");
1478         fhEtaPhiPtFracDecayIso[icone][ipt]->SetYTitle("#phi");
1479         outputContainer->Add(fhEtaPhiPtFracDecayIso[icone][ipt]) ;
1480
1481
1482         snprintf(name, buffersize,"hEtaPhiPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
1483         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1484         fhEtaPhiPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1485         fhEtaPhiPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
1486         fhEtaPhiPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
1487         outputContainer->Add(fhEtaPhiPtSumDecayIso[icone][ipt]) ;
1488
1489         snprintf(name, buffersize,"hEtaPhiSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
1490         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for density R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1491         fhEtaPhiSumDensityDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1492         fhEtaPhiSumDensityDecayIso[icone][ipt]->SetXTitle("#eta");
1493         fhEtaPhiSumDensityDecayIso[icone][ipt]->SetYTitle("#phi");
1494         outputContainer->Add(fhEtaPhiSumDensityDecayIso[icone][ipt]) ;
1495
1496         snprintf(name, buffersize,"hEtaPhiFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
1497         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for FracPtSum R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1498         fhEtaPhiFracPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1499         fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
1500         fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
1501         outputContainer->Add(fhEtaPhiFracPtSumDecayIso[icone][ipt]) ;
1502
1503
1504         if(IsDataMC())
1505         {
1506           snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);
1507           snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1508           fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1509           fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1510           outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
1511
1512           snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);
1513           snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1514           fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1515           fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1516           outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
1517
1518           snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);
1519           snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1520           fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1521           fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1522           outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
1523
1524           snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);
1525           snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1526           fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1527           fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1528           outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
1529
1530           snprintf(name, buffersize,"hPtThresMCPi0_Cone_%d_Pt%d",icone,ipt);
1531           snprintf(title, buffersize,"Isolated candidate Pi0 p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1532           fhPtThresIsolatedPi0[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1533           fhPtThresIsolatedPi0[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1534           outputContainer->Add(fhPtThresIsolatedPi0[icone][ipt]) ;
1535
1536           snprintf(name, buffersize,"hPtFracMCPi0_Cone_%d_Pt%d",icone,ipt);
1537           snprintf(title, buffersize,"Isolated candidate Pi0 p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1538           fhPtFracIsolatedPi0[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1539           fhPtFracIsolatedPi0[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1540           outputContainer->Add(fhPtFracIsolatedPi0[icone][ipt]) ;
1541
1542           snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
1543           snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1544           fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1545           fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1546           outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
1547
1548           snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
1549           snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1550           fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1551           fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1552           outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
1553
1554           snprintf(name, buffersize,"hPtThresMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
1555           snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1556           fhPtThresIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1557           fhPtThresIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1558           outputContainer->Add(fhPtThresIsolatedEtaDecay[icone][ipt]) ;
1559
1560           snprintf(name, buffersize,"hPtFracMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
1561           snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1562           fhPtFracIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1563           fhPtFracIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1564           outputContainer->Add(fhPtFracIsolatedEtaDecay[icone][ipt]) ;
1565
1566
1567           snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
1568           snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1569           fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1570           fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1571           outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
1572
1573           snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
1574           snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1575           fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1576           fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1577           outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
1578
1579 //          snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);
1580 //          snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1581 //          fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1582 //          fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1583 //          outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
1584
1585 //          snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);
1586 //          snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1587 //          fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1588 //          fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1589 //          outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
1590
1591           snprintf(name, buffersize,"hPtThresMCHadron_Cone_%d_Pt%d",icone,ipt);
1592           snprintf(title, buffersize,"Isolated candidate Hadron p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1593           fhPtThresIsolatedHadron[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1594           fhPtThresIsolatedHadron[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1595           outputContainer->Add(fhPtThresIsolatedHadron[icone][ipt]) ;
1596
1597           snprintf(name, buffersize,"hPtFracMCHadron_Cone_%d_Pt%d",icone,ipt);
1598           snprintf(title, buffersize,"Isolated candidate Hadron p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1599           fhPtFracIsolatedHadron[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1600           fhPtFracIsolatedHadron[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1601           outputContainer->Add(fhPtFracIsolatedHadron[icone][ipt]) ;
1602
1603         }//Histos with MC
1604       }//icone loop
1605     }//ipt loop
1606   }
1607
1608   if(fFillPileUpHistograms)
1609   {
1610     for (Int_t i = 0; i < 7 ; i++)
1611     {
1612       fhEIsoPileUp[i]   = new TH1F(Form("hEPileUp%s",pileUpName[i].Data()),
1613                                 Form("Number of isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f, pile-up event by %s",r,ptthre,ptfrac,pileUpName[i].Data()),
1614                                 nptbins,ptmin,ptmax);
1615       fhEIsoPileUp[i]->SetYTitle("dN / dE");
1616       fhEIsoPileUp[i]->SetXTitle("E (GeV/c)");
1617       outputContainer->Add(fhEIsoPileUp[i]) ;
1618
1619       fhPtIsoPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
1620                                 Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr}= %2.2f, pile-up event by %s ",r,ptthre,ptfrac,pileUpName[i].Data()),
1621                                 nptbins,ptmin,ptmax);
1622       fhPtIsoPileUp[i]->SetYTitle("dN / p_{T}");
1623       fhPtIsoPileUp[i]->SetXTitle("p_{T} (GeV/c)");
1624       outputContainer->Add(fhPtIsoPileUp[i]) ;
1625
1626       fhENoIsoPileUp[i]   = new TH1F(Form("hENoIsoPileUp%s",pileUpName[i].Data()),
1627                                   Form("Number of not isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f, pile-up event by %s",r,ptthre,ptfrac,pileUpName[i].Data()),
1628                                   nptbins,ptmin,ptmax);
1629       fhENoIsoPileUp[i]->SetYTitle("dN / dE");
1630       fhENoIsoPileUp[i]->SetXTitle("E (GeV/c)");
1631       outputContainer->Add(fhENoIsoPileUp[i]) ;
1632
1633       fhPtNoIsoPileUp[i]  = new TH1F(Form("hPtNoIsoPileUp%s",pileUpName[i].Data()),
1634                                   Form("Number of not isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr}= %2.2f, pile-up event by %s ",r,ptthre,ptfrac,pileUpName[i].Data()),
1635                                   nptbins,ptmin,ptmax);
1636       fhPtNoIsoPileUp[i]->SetYTitle("dN / p_{T}");
1637       fhPtNoIsoPileUp[i]->SetXTitle("p_{T} (GeV/c)");
1638       outputContainer->Add(fhPtNoIsoPileUp[i]) ;
1639     }
1640
1641     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1642     fhTimeENoCut->SetXTitle("E (GeV)");
1643     fhTimeENoCut->SetYTitle("time (ns)");
1644     outputContainer->Add(fhTimeENoCut);
1645
1646     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1647     fhTimeESPD->SetXTitle("E (GeV)");
1648     fhTimeESPD->SetYTitle("time (ns)");
1649     outputContainer->Add(fhTimeESPD);
1650
1651     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1652     fhTimeESPDMulti->SetXTitle("E (GeV)");
1653     fhTimeESPDMulti->SetYTitle("time (ns)");
1654     outputContainer->Add(fhTimeESPDMulti);
1655
1656     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1657     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1658     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1659     outputContainer->Add(fhTimeNPileUpVertSPD);
1660
1661     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
1662     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1663     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1664     outputContainer->Add(fhTimeNPileUpVertTrack);
1665
1666     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1667     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1668     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1669     outputContainer->Add(fhTimeNPileUpVertContributors);
1670
1671     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);
1672     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1673     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1674     outputContainer->Add(fhTimePileUpMainVertexZDistance);
1675
1676     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
1677     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1678     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1679     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
1680   }
1681
1682   return outputContainer ;
1683
1684 }
1685
1686 //__________________________________
1687 void AliAnaParticleIsolation::Init()
1688 {
1689   // Do some checks and init stuff
1690
1691   // In case of several cone and thresholds analysis, open the cuts for the filling of the
1692   // track and cluster reference arrays in cone when done in the MakeAnalysisFillAOD().
1693   // The different cones, thresholds are tested for this list of tracks, clusters.
1694   if(fMakeSeveralIC)
1695   {
1696     printf("AliAnaParticleIsolation::Init() - Open default isolation cuts for multiple Isolation analysis\n");
1697     GetIsolationCut()->SetPtThreshold(100);
1698     GetIsolationCut()->SetPtFraction(100);
1699     GetIsolationCut()->SetConeSize(1);
1700   }
1701 }
1702
1703 //____________________________________________
1704 void AliAnaParticleIsolation::InitParameters()
1705 {
1706
1707   //Initialize the parameters of the analysis.
1708   SetInputAODName("PWG4Particle");
1709   SetAODObjArrayName("IsolationCone");
1710   AddToHistogramsName("AnaIsolation_");
1711
1712   fCalorimeter = "PHOS" ;
1713   fReMakeIC = kFALSE ;
1714   fMakeSeveralIC = kFALSE ;
1715
1716   //----------- Several IC-----------------
1717   fNCones             = 5 ;
1718   fNPtThresFrac       = 5 ;
1719   fConeSizes      [0] = 0.1;     fConeSizes      [1] = 0.2;   fConeSizes      [2] = 0.3; fConeSizes      [3] = 0.4;  fConeSizes      [4] = 0.5;
1720   fPtThresholds   [0] = 1.;      fPtThresholds   [1] = 2.;    fPtThresholds   [2] = 3.;  fPtThresholds   [3] = 4.;   fPtThresholds   [4] = 5.;
1721   fPtFractions    [0] = 0.05;    fPtFractions    [1] = 0.075; fPtFractions    [2] = 0.1; fPtFractions    [3] = 1.25; fPtFractions    [4] = 1.5;
1722   fSumPtThresholds[0] = 1.;      fSumPtThresholds[1] = 2.;    fSumPtThresholds[2] = 3.;  fSumPtThresholds[3] = 4.;   fSumPtThresholds[4] = 5.;
1723
1724   //------------- Histograms settings -------
1725   fHistoNPtSumBins = 100 ;
1726   fHistoPtSumMax   = 50 ;
1727   fHistoPtSumMin   = 0.  ;
1728
1729   fHistoNPtInConeBins = 100 ;
1730   fHistoPtInConeMax   = 50 ;
1731   fHistoPtInConeMin   = 0.  ;
1732
1733 }
1734
1735 //__________________________________________________
1736 void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
1737 {
1738   //Do analysis and fill aods
1739   //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
1740
1741   if(!GetInputAODBranch())
1742   {
1743     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1744     abort();
1745   }
1746
1747   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1748   {
1749     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
1750     abort();
1751   }
1752
1753   Int_t n = 0, nfrac = 0;
1754   Bool_t  isolated  = kFALSE ;
1755   Float_t coneptsum = 0 ;
1756   TObjArray * pl    = 0x0; ;
1757
1758   //Select the calorimeter for candidate isolation with neutral particles
1759   if      (fCalorimeter == "PHOS" )
1760     pl = GetPHOSClusters();
1761   else if (fCalorimeter == "EMCAL")
1762     pl = GetEMCALClusters();
1763
1764   //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
1765   Double_t ptLeading = 0. ;
1766   Int_t    idLeading = -1 ;
1767   TLorentzVector mom ;
1768   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1769   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
1770
1771   for(Int_t iaod = 0; iaod < naod; iaod++)
1772   {
1773     AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1774
1775     //If too small or too large pt, skip
1776     if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ;
1777
1778     //check if it is low pt trigger particle
1779     if((aodinput->Pt() < GetIsolationCut()->GetPtThreshold() ||
1780         aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold()) &&
1781        !fMakeSeveralIC)
1782     {
1783       continue ; //trigger should not come from underlying event
1784     }
1785
1786     //vertex cut in case of mixing
1787     Int_t check = CheckMixedEventVertex(aodinput->GetCaloLabel(0), aodinput->GetTrackLabel(0));
1788     if(check ==  0) continue;
1789     if(check == -1) return;
1790
1791     //find the leading particles with highest momentum
1792     if ( aodinput->Pt() > ptLeading )
1793     {
1794       ptLeading = aodinput->Pt() ;
1795       idLeading = iaod ;
1796     }
1797
1798     aodinput->SetLeadingParticle(kFALSE);
1799
1800   }//finish searching for leading trigger particle
1801
1802   // Check isolation of leading particle
1803   if(idLeading < 0) return;
1804
1805   AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(idLeading));
1806   aodinput->SetLeadingParticle(kTRUE);
1807
1808   // Check isolation only of clusters in fiducial region
1809   if(IsFiducialCutOn())
1810   {
1811     Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodinput->Momentum(),aodinput->GetDetector()) ;
1812     if(! in ) return ;
1813   }
1814
1815   //After cuts, study isolation
1816   n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
1817   GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,
1818                                       GetReader(), GetCaloPID(),
1819                                       kTRUE, aodinput, GetAODObjArrayName(),
1820                                       n,nfrac,coneptsum, isolated);
1821
1822   if(!fMakeSeveralIC) aodinput->SetIsolated(isolated);
1823
1824   if(GetDebug() > 1)
1825   {
1826     if(isolated)printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
1827     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
1828   }
1829
1830 }
1831
1832 //_________________________________________________________
1833 void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
1834 {
1835   //Do analysis and fill histograms
1836
1837   Int_t   n = 0, nfrac = 0;
1838   Bool_t  isolated  = kFALSE ;
1839   Float_t coneptsum = 0 ;
1840   Float_t etaUEptsum = 0 ;
1841   Float_t phiUEptsum = 0 ;
1842
1843   Double_t bz = GetReader()->GetInputEvent()->GetMagneticField();
1844
1845   //Loop on stored AOD
1846   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1847   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
1848
1849   //Get vertex for photon momentum calculation
1850   Double_t vertex[] = {0,0,0} ; //vertex ;
1851   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1852   {
1853     GetReader()->GetVertex(vertex);
1854   }
1855
1856   Float_t cen = GetEventCentrality();
1857   Float_t ep  = GetEventPlaneAngle();
1858
1859   for(Int_t iaod = 0; iaod < naod ; iaod++)
1860   {
1861     AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1862
1863     if(!aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
1864
1865     // Check isolation only of clusters in fiducial region
1866     if(IsFiducialCutOn())
1867     {
1868       Bool_t in = GetFiducialCut()->IsInFiducialCut(*aod->Momentum(),aod->GetDetector()) ;
1869       if(! in ) continue ;
1870     }
1871
1872     Bool_t  isolation  = aod->IsIsolated();
1873     Bool_t  decay      = aod->IsTagged();
1874     Float_t energy     = aod->E();
1875     Float_t pt         = aod->Pt();
1876     Float_t phi        = aod->Phi();
1877     Float_t eta        = aod->Eta();
1878     Float_t conesize   = GetIsolationCut()->GetConeSize();
1879
1880     //Recover reference arrays with clusters and tracks
1881     TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
1882     TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
1883
1884     //If too small or too large pt, skip
1885     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
1886
1887     // --- In case of redoing isolation from delta AOD ----
1888
1889     if(fMakeSeveralIC)
1890     {
1891       //Analysis of multiple IC at same time
1892       MakeSeveralICAnalysis(aod);
1893       continue;
1894     }
1895     else if(fReMakeIC)
1896     {
1897       //In case a more strict IC is needed in the produced AOD
1898       n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
1899       GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
1900                                           GetReader(), GetCaloPID(),
1901                                           kFALSE, aod, "",
1902                                           n,nfrac,coneptsum, isolated);
1903       fhConeSumPt->Fill(pt,coneptsum);
1904       if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
1905     }
1906
1907     // ---------------------------------------------------
1908
1909     //Fill pt distribution of particles in cone
1910     //Tracks
1911     coneptsum = 0;
1912     Double_t sumptFR = 0. ;
1913     TObjArray * trackList   = GetCTSTracks() ;
1914     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
1915     {
1916       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
1917
1918       if(!track)
1919       {
1920         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
1921         continue;
1922       }
1923
1924       //fill histogram for UE in phi band
1925       if(track->Eta() > (eta-conesize) && track->Eta()  < (eta+conesize))
1926       {
1927         phiUEptsum+=track->Pt();
1928         fhPhiBand->Fill(track->Eta(),track->Phi());
1929       }
1930       //fill histogram for UE in eta band in EMCal acceptance
1931       if(track->Phi() > (phi-conesize) && track->Phi() < (phi+conesize) && track->Eta() > -0.6 && track->Eta() < 0.6)
1932       {
1933         etaUEptsum+=track->Pt();
1934         fhEtaBand->Fill(track->Eta(),track->Phi());
1935       }
1936
1937       //fill the histograms at forward range
1938       Double_t dPhi = phi - track->Phi() + TMath::PiOver2();
1939       Double_t dEta = eta - track->Eta();
1940       Double_t arg  = dPhi*dPhi + dEta*dEta;
1941       if(TMath::Sqrt(arg) < conesize)
1942       {
1943         fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1944         sumptFR+=track->Pt();
1945       }
1946
1947       dPhi = phi - track->Phi() - TMath::PiOver2();
1948       arg  = dPhi*dPhi + dEta*dEta;
1949       if(TMath::Sqrt(arg) < conesize)
1950       {
1951         fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1952         sumptFR+=track->Pt();
1953       }
1954     }
1955
1956     fhFRConeSumPt->Fill(pt,sumptFR);
1957     if(reftracks)
1958     {
1959       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
1960       {
1961         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
1962         Float_t pTtrack = track->Pt();
1963         fhPtInCone     ->Fill(pt,pTtrack);
1964         fhPtTrackInCone->Fill(pt,pTtrack);
1965         if(fFillPileUpHistograms)
1966         {
1967           ULong_t status = track->GetStatus();
1968           Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1969           //Double32_t tof = track->GetTOFsignal()*1e-3;
1970           Int_t trackBC = track->GetTOFBunchCrossing(bz);
1971
1972           if     ( okTOF && trackBC!=0 )fhPtTrackInConeOtherBC->Fill(pt,pTtrack);
1973           else if( okTOF && trackBC==0 )fhPtTrackInConeBC0    ->Fill(pt,pTtrack);
1974
1975           Int_t vtxBC = GetReader()->GetVertexBC();
1976           if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA) fhPtTrackInConeVtxBC0->Fill(pt,pTtrack);
1977
1978           if(GetReader()->IsPileUpFromSPD())             { fhPtInConePileUp[0]->Fill(pt,pTtrack);
1979             if(okTOF && trackBC!=0 )           fhPtTrackInConeOtherBCPileUpSPD->Fill(pt,pTtrack);
1980             if(okTOF && trackBC==0 )           fhPtTrackInConeBC0PileUpSPD    ->Fill(pt,pTtrack); }
1981           if(GetReader()->IsPileUpFromEMCal())             fhPtInConePileUp[1]->Fill(pt,pTtrack);
1982           if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtInConePileUp[2]->Fill(pt,pTtrack);
1983           if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtInConePileUp[3]->Fill(pt,pTtrack);
1984           if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtInConePileUp[4]->Fill(pt,pTtrack);
1985           if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtInConePileUp[5]->Fill(pt,pTtrack);
1986           if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(pt,pTtrack);
1987         }
1988
1989         if (GetEventCentrality()) fhPtInConeCent->Fill(GetEventCentrality(),pTtrack);
1990         coneptsum+=pTtrack;
1991       }
1992     }
1993
1994     //CaloClusters
1995     if(refclusters)
1996     {
1997       TLorentzVector mom ;
1998       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
1999       {
2000         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
2001         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
2002
2003         //fill histogram for UE in phi band
2004         if(mom.Eta() > (eta-conesize) && mom.Eta()  < (eta+conesize))
2005         {
2006           phiUEptsum+=mom.Pt();
2007           fhPhiBand->Fill(mom.Eta(),mom.Phi());
2008         }
2009         //fill histogram for UE in eta band in EMCal acceptance
2010         if(mom.Phi() > (phi-conesize) && mom.Phi() < (phi+conesize))
2011         {
2012           etaUEptsum+=mom.Pt();
2013           fhEtaBand->Fill(mom.Eta(),mom.Phi());
2014         }
2015
2016         fhPtInCone->Fill(pt, mom.Pt());
2017         if(fFillPileUpHistograms)
2018         {
2019           if(GetReader()->IsPileUpFromSPD())               fhPtInConePileUp[0]->Fill(pt,mom.Pt());
2020           if(GetReader()->IsPileUpFromEMCal())             fhPtInConePileUp[1]->Fill(pt,mom.Pt());
2021           if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtInConePileUp[2]->Fill(pt,mom.Pt());
2022           if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtInConePileUp[3]->Fill(pt,mom.Pt());
2023           if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtInConePileUp[4]->Fill(pt,mom.Pt());
2024           if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtInConePileUp[5]->Fill(pt,mom.Pt());
2025           if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(pt,mom.Pt());
2026         }
2027
2028         if (GetEventCentrality()) fhPtInConeCent->Fill(GetEventCentrality(),mom.Pt());
2029         coneptsum+=mom.Pt();
2030       }
2031     }
2032
2033     //normalize phi/eta band per area unit
2034     fhPhiUEConeSumPt->Fill(pt, phiUEptsum*(TMath::Pi()*conesize*conesize)/(2*conesize*2*TMath::Pi()));
2035     fhEtaUEConeSumPt->Fill(pt, etaUEptsum*(TMath::Pi()*conesize*conesize)/(2*conesize*2*0.6));
2036
2037     Double_t sumPhiUESub = coneptsum-(phiUEptsum*(TMath::Pi()*conesize*conesize)/(2*conesize*2*TMath::Pi()));
2038     Double_t sumEtaUESub = coneptsum-(etaUEptsum*(TMath::Pi()*conesize*conesize)/(2*conesize*2*0.6));
2039
2040     fhConeSumPtPhiUESub->Fill(pt,sumPhiUESub);
2041     fhConeSumPtEtaUESub->Fill(pt,sumEtaUESub);
2042
2043
2044     if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);
2045
2046     if(!fReMakeIC) fhConeSumPt->Fill(pt,coneptsum);
2047
2048     Int_t mcTag = aod->GetTag() ;
2049     Int_t clID  = aod->GetCaloLabel(0) ;
2050     Int_t nlm = aod->GetFiducialArea();
2051     if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f\n",pt, eta, phi);
2052
2053     FillTrackMatchingShowerShapeControlHistograms(isolation, clID,nlm,mcTag,reftracks,refclusters,aod,GetReader(), GetCaloPID());
2054
2055     if(isolation)
2056     {
2057       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED: fill histograms\n", iaod);
2058
2059       // Fill histograms to undertand pile-up before other cuts applied
2060       // Remember to relax time cuts in the reader
2061       FillPileUpHistograms(clID);
2062
2063       fhEIso      ->Fill(energy);
2064       fhPtIso     ->Fill(pt);
2065       fhPhiIso    ->Fill(pt,phi);
2066       fhEtaIso    ->Fill(pt,eta);
2067       fhEtaPhiIso ->Fill(eta,phi);
2068
2069       fhPtNLocMaxIso    ->Fill(pt,nlm) ;
2070       fhPtCentralityIso ->Fill(pt,cen) ;
2071       fhPtEventPlaneIso ->Fill(pt,ep ) ;
2072
2073       if(decay)
2074       {
2075         fhPtDecayIso->Fill(pt);
2076         fhEtaPhiDecayIso->Fill(eta,phi);
2077       }
2078
2079       if(fFillPileUpHistograms)
2080       {
2081         if(GetReader()->IsPileUpFromSPD())               { fhEIsoPileUp[0] ->Fill(energy) ; fhPtIsoPileUp[0]->Fill(pt) ; }
2082         if(GetReader()->IsPileUpFromEMCal())             { fhEIsoPileUp[1] ->Fill(energy) ; fhPtIsoPileUp[1]->Fill(pt) ; }
2083         if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhEIsoPileUp[2] ->Fill(energy) ; fhPtIsoPileUp[2]->Fill(pt) ; }
2084         if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhEIsoPileUp[3] ->Fill(energy) ; fhPtIsoPileUp[3]->Fill(pt) ; }
2085         if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhEIsoPileUp[4] ->Fill(energy) ; fhPtIsoPileUp[4]->Fill(pt) ; }
2086         if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhEIsoPileUp[5] ->Fill(energy) ; fhPtIsoPileUp[5]->Fill(pt) ; }
2087         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhEIsoPileUp[6] ->Fill(energy) ; fhPtIsoPileUp[6]->Fill(pt) ; }
2088       }
2089
2090       if(IsDataMC())
2091       {
2092
2093         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
2094         {
2095           fhPtIsoMCPhoton  ->Fill(pt);
2096         }
2097
2098         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))
2099         {
2100           fhPtIsoPrompt  ->Fill(pt);
2101           fhPhiIsoPrompt ->Fill(pt,phi);
2102           fhEtaIsoPrompt ->Fill(pt,eta);
2103         }
2104         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation))
2105         {
2106           fhPtIsoFragmentation  ->Fill(pt);
2107           fhPhiIsoFragmentation ->Fill(pt,phi);
2108           fhEtaIsoFragmentation ->Fill(pt,eta);
2109         }
2110         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))
2111         {
2112           fhPtIsoPi0  ->Fill(pt);
2113           fhPhiIsoPi0 ->Fill(pt,phi);
2114           fhEtaIsoPi0 ->Fill(pt,eta);
2115         }
2116         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))
2117         {
2118           fhPtIsoPi0Decay  ->Fill(pt);
2119           fhPhiIsoPi0Decay ->Fill(pt,phi);
2120           fhEtaIsoPi0Decay ->Fill(pt,eta);
2121         }
2122         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))
2123         {
2124           fhPtIsoEtaDecay  ->Fill(pt);
2125           fhPhiIsoEtaDecay ->Fill(pt,phi);
2126           fhEtaIsoEtaDecay ->Fill(pt,eta);
2127         }
2128         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))
2129         {
2130           fhPtIsoOtherDecay  ->Fill(pt);
2131           fhPhiIsoOtherDecay ->Fill(pt,phi);
2132           fhEtaIsoOtherDecay ->Fill(pt,eta);
2133         }
2134         //        else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))
2135         //        {
2136         //          fhPtIsoConversion  ->Fill(pt);
2137         //          fhPhiIsoConversion ->Fill(pt,phi);
2138         //          fhEtaIsoConversion ->Fill(pt,eta);
2139         //        }
2140         else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))// anything else but electrons
2141         {
2142           fhPtIsoHadron  ->Fill(pt);
2143           fhPhiIsoHadron ->Fill(pt,phi);
2144           fhEtaIsoHadron ->Fill(pt,eta);
2145         }
2146       }//Histograms with MC
2147
2148     }//Isolated histograms
2149     else // NON isolated
2150     {
2151       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d NOT ISOLATED, fill histograms\n", iaod);
2152
2153       fhENoIso        ->Fill(energy);
2154       fhPtNoIso       ->Fill(pt);
2155       fhEtaPhiNoIso   ->Fill(eta,phi);
2156       fhPtNLocMaxNoIso->Fill(pt,nlm);
2157
2158       if(fFillPileUpHistograms)
2159       {
2160         if(GetReader()->IsPileUpFromSPD())                { fhENoIsoPileUp[0] ->Fill(energy) ; fhPtNoIsoPileUp[0]->Fill(pt) ; }
2161         if(GetReader()->IsPileUpFromEMCal())              { fhENoIsoPileUp[1] ->Fill(energy) ; fhPtNoIsoPileUp[1]->Fill(pt) ; }
2162         if(GetReader()->IsPileUpFromSPDOrEMCal())         { fhENoIsoPileUp[2] ->Fill(energy) ; fhPtNoIsoPileUp[2]->Fill(pt) ; }
2163         if(GetReader()->IsPileUpFromSPDAndEMCal())        { fhENoIsoPileUp[3] ->Fill(energy) ; fhPtNoIsoPileUp[3]->Fill(pt) ; }
2164         if(GetReader()->IsPileUpFromSPDAndNotEMCal())     { fhENoIsoPileUp[4] ->Fill(energy) ; fhPtNoIsoPileUp[4]->Fill(pt) ; }
2165         if(GetReader()->IsPileUpFromEMCalAndNotSPD())     { fhENoIsoPileUp[5] ->Fill(energy) ; fhPtNoIsoPileUp[5]->Fill(pt) ; }
2166         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhENoIsoPileUp[6] ->Fill(energy) ; fhPtNoIsoPileUp[6]->Fill(pt) ; }
2167       }
2168
2169       if(decay)
2170       {
2171         fhPtDecayNoIso    ->Fill(pt);
2172         fhEtaPhiDecayNoIso->Fill(eta,phi);
2173       }
2174
2175       if(IsDataMC())
2176       {
2177         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))        fhPtNoIsoMCPhoton     ->Fill(pt);
2178         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))           fhPtNoIsoPi0          ->Fill(pt);
2179         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtNoIsoPi0Decay     ->Fill(pt);
2180         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtNoIsoEtaDecay     ->Fill(pt);
2181         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtNoIsoOtherDecay   ->Fill(pt);
2182         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))        fhPtNoIsoPrompt       ->Fill(pt);
2183         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) fhPtNoIsoFragmentation->Fill(pt);
2184         //        else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(pt);
2185         else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))     fhPtNoIsoHadron      ->Fill(pt);
2186       }
2187     }
2188   }// aod loop
2189
2190 }
2191
2192
2193 //_____________________________________________________________________________________
2194 void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph)
2195 {
2196
2197   //Isolation Cut Analysis for both methods and different pt cuts and cones
2198   Float_t ptC   = ph->Pt();
2199   Float_t etaC  = ph->Eta();
2200   Float_t phiC  = ph->Phi();
2201   Int_t   tag   = ph->GetTag();
2202   Bool_t  decay = ph->IsTagged();
2203
2204   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f\n",ptC);
2205
2206   //Keep original setting used when filling AODs, reset at end of analysis
2207   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
2208   Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();
2209   Float_t rorg       = GetIsolationCut()->GetConeSize();
2210
2211   Float_t coneptsum = 0 ;
2212   Int_t   n    [10][10];//[fNCones][fNPtThresFrac];
2213   Int_t   nfrac[10][10];//[fNCones][fNPtThresFrac];
2214   Bool_t  isolated  = kFALSE;
2215   Int_t   nCone     = 0;
2216   Int_t   nFracCone = 0;
2217
2218   // fill hist with all particles before isolation criteria
2219   fhPtNoIso    ->Fill(ptC);
2220   fhEtaPhiNoIso->Fill(etaC,phiC);
2221
2222   if(IsDataMC())
2223   {
2224     if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))        fhPtNoIsoMCPhoton     ->Fill(ptC);
2225     if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtNoIsoPi0          ->Fill(ptC);
2226     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtNoIsoPi0Decay     ->Fill(ptC);
2227     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtNoIsoEtaDecay     ->Fill(ptC);
2228     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtNoIsoOtherDecay   ->Fill(ptC);
2229     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtNoIsoPrompt       ->Fill(ptC);
2230     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtNoIsoFragmentation->Fill(ptC);
2231 //    else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(ptC);
2232     else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))    fhPtNoIsoHadron      ->Fill(ptC);
2233   }
2234
2235   if(decay)
2236   {
2237     fhPtDecayNoIso    ->Fill(ptC);
2238     fhEtaPhiDecayNoIso->Fill(etaC,phiC);
2239   }
2240   //Get vertex for photon momentum calculation
2241   Double_t vertex[] = {0,0,0} ; //vertex ;
2242   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2243     {
2244       GetReader()->GetVertex(vertex);
2245     }
2246
2247   //Loop on cone sizes
2248   for(Int_t icone = 0; icone<fNCones; icone++)
2249   {
2250     //Recover reference arrays with clusters and tracks
2251     TObjArray * refclusters = ph->GetObjArray(GetAODObjArrayName()+"Clusters");
2252     TObjArray * reftracks   = ph->GetObjArray(GetAODObjArrayName()+"Tracks");
2253
2254     //If too small or too large pt, skip
2255     if(ptC < GetMinPt() || ptC > GetMaxPt() ) continue ;
2256
2257    //In case a more strict IC is needed in the produced AOD
2258
2259     nCone=0; nFracCone = 0; isolated = kFALSE; coneptsum = 0;
2260
2261   GetIsolationCut()->SetSumPtThreshold(100);
2262   GetIsolationCut()->SetPtThreshold(100);
2263   GetIsolationCut()->SetPtFraction(100);
2264   GetIsolationCut()->SetConeSize(fConeSizes[icone]);
2265   GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
2266                                           GetReader(), GetCaloPID(),
2267                                           kFALSE, ph, "",
2268                                       nCone,nFracCone,coneptsum, isolated);
2269
2270
2271     fhSumPtLeadingPt[icone]->Fill(ptC,coneptsum);
2272
2273     // retreive pt tracks to fill histo vs. pt leading
2274     //Fill pt distribution of particles in cone
2275     //fhPtLeadingPt(),fhFRSumPtLeadingPt(),fhFRPtLeadingPt(),
2276
2277     //Tracks
2278     coneptsum = 0;
2279     Double_t sumptFR = 0. ;
2280     TObjArray * trackList   = GetCTSTracks() ;
2281     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
2282     {
2283       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
2284       //fill the histograms at forward range
2285       if(!track)
2286       {
2287         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
2288         continue;
2289       }
2290
2291       Double_t dPhi = phiC - track->Phi() + TMath::PiOver2();
2292       Double_t dEta = etaC - track->Eta();
2293       Double_t arg  = dPhi*dPhi + dEta*dEta;
2294       if(TMath::Sqrt(arg) < fConeSizes[icone])
2295       {
2296         fhFRPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
2297         sumptFR+=track->Pt();
2298       }
2299
2300       dPhi = phiC - track->Phi() - TMath::PiOver2();
2301       arg  = dPhi*dPhi + dEta*dEta;
2302       if(TMath::Sqrt(arg) < fConeSizes[icone])
2303       {
2304         fhFRPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
2305         sumptFR+=track->Pt();
2306       }
2307     }
2308     fhFRSumPtLeadingPt[icone]->Fill(ptC,sumptFR);
2309     if(reftracks)
2310     {
2311       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
2312       {
2313         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
2314         fhPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
2315         coneptsum+=track->Pt();
2316       }
2317     }
2318     //CaloClusters
2319     if(refclusters)
2320     {
2321       TLorentzVector mom ;
2322       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
2323       {
2324         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
2325         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
2326
2327         fhPtLeadingPt[icone]->Fill(ptC, mom.Pt());
2328         coneptsum+=mom.Pt();
2329       }
2330     }
2331     ///////////////////
2332
2333
2334     //Loop on ptthresholds
2335     for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++)
2336     {
2337       n    [icone][ipt]=0;
2338       nfrac[icone][ipt]=0;
2339       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
2340       GetIsolationCut()->SetPtFraction(fPtFractions[ipt]) ;
2341       GetIsolationCut()->SetSumPtThreshold(fSumPtThresholds[ipt]);
2342
2343       GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
2344                                           GetReader(), GetCaloPID(),
2345                                           kFALSE, ph, "",
2346                                           n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
2347
2348       if(!isolated) continue;
2349       //Normal ptThreshold cut
2350
2351       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - cone size %1.1f, ptThres  %1.1f, sumptThresh  %1.1f, n %d, nfrac %d, coneptsum %2.2f, isolated %d\n",
2352                                 fConeSizes[icone],fPtThresholds[ipt],fSumPtThresholds[ipt],n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
2353       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - pt %1.1f, eta %1.1f, phi %1.1f\n",ptC, etaC, phiC);
2354
2355       if(n[icone][ipt] == 0)
2356       {
2357         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling pt threshold loop\n");
2358         fhPtThresIsolated[icone][ipt]->Fill(ptC);
2359         fhEtaPhiPtThresIso[icone][ipt]->Fill(etaC,phiC);
2360
2361         if(decay)
2362         {
2363           fhPtPtThresDecayIso[icone][ipt]->Fill(ptC);
2364           //      fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
2365         }
2366
2367         if(IsDataMC())
2368         {
2369           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtThresIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
2370 //          else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtThresIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
2371           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
2372           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtThresIsolatedPi0[icone][ipt]           ->Fill(ptC) ;
2373           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtThresIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
2374           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtThresIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
2375           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtThresIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
2376           else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))      fhPtThresIsolatedHadron[icone][ipt]      ->Fill(ptC) ;
2377         }
2378       }
2379
2380       // pt in cone fraction
2381       if(nfrac[icone][ipt] == 0)
2382       {
2383         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling frac loop\n");
2384         fhPtFracIsolated[icone][ipt]->Fill(ptC);
2385         fhEtaPhiPtFracIso[icone][ipt]->Fill(etaC,phiC);
2386
2387         if(decay)
2388         {
2389           fhPtPtFracDecayIso[icone][ipt]->Fill(ptC);
2390           fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
2391         }
2392
2393         if(IsDataMC())
2394         {
2395           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtFracIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
2396 //          else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtFracIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
2397           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
2398           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtFracIsolatedPi0[icone][ipt]          ->Fill(ptC) ;
2399           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtFracIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
2400           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtFracIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
2401           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtFracIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
2402           else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))      fhPtFracIsolatedHadron[icone][ipt]->Fill(ptC) ;
2403         }
2404       }
2405
2406       if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - checking IC method : %i\n",GetIsolationCut()->GetICMethod());
2407
2408       //Pt threshold on pt cand/ sum in cone histograms
2409       if(coneptsum<fSumPtThresholds[ipt])
2410       {//      if((GetIsolationCut()->GetICMethod())==1){//kSumPtIC){
2411         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling sum loop\n");
2412         fhPtSumIsolated[icone][ipt]->Fill(ptC) ;
2413         fhEtaPhiPtSumIso[icone][ipt]->Fill(etaC, phiC) ;
2414         if(decay)
2415         {
2416           fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
2417           fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
2418         }
2419       }
2420
2421     // pt sum pt frac method
2422 //    if( ((fPtFractions[ipt]*ptC < fSumPtThresholds[ipt]) && (coneptsum < fSumPtThresholds[ipt])) || ((fPtFractions[ipt]*ptC > fSumPtThresholds[ipt]) && (coneptsum < fPtFractions[ipt]*ptC)) )
2423
2424       if(coneptsum < fPtFractions[ipt]*ptC)
2425        {
2426         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling PtFrac PtSum loop\n");
2427         fhPtFracPtSumIso[icone][ipt]->Fill(ptC) ;
2428         fhEtaPhiFracPtSumIso[icone][ipt]->Fill(etaC,phiC) ;
2429
2430         if(decay)
2431         {
2432           fhPtFracPtSumDecayIso[icone][ipt]->Fill(ptC);
2433           fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
2434         }
2435       }
2436
2437       // density method
2438       Float_t cellDensity = GetIsolationCut()->GetCellDensity( ph, GetReader());
2439       if(coneptsum<fSumPtThresholds[ipt]*cellDensity)
2440       {//(GetIsolationCut()->GetICMethod())==4){//kSumDensityIC) {
2441         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling density loop\n");
2442         fhPtSumDensityIso[icone][ipt]->Fill(ptC) ;
2443         fhEtaPhiSumDensityIso[icone][ipt]->Fill(etaC,phiC) ;
2444
2445         if(decay)
2446         {
2447           fhPtSumDensityDecayIso[icone][ipt]->Fill(ptC);
2448           fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
2449         }
2450
2451       }
2452     }//pt thresh loop
2453
2454     if(IsDataMC())
2455     {
2456       if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtSumIsolatedPrompt[icone]       ->Fill(ptC,coneptsum) ;
2457 //      else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtSumIsolatedConversion[icone]   ->Fill(ptC,coneptsum) ;
2458       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
2459       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))           fhPtSumIsolatedPi0[icone]          ->Fill(ptC,coneptsum) ;
2460       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtSumIsolatedPi0Decay[icone]     ->Fill(ptC,coneptsum) ;
2461       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtSumIsolatedEtaDecay[icone]     ->Fill(ptC,coneptsum) ;
2462       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtSumIsolatedOtherDecay[icone]   ->Fill(ptC,coneptsum) ;
2463       else if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))      fhPtSumIsolatedHadron[icone]->Fill(ptC,coneptsum) ;
2464     }
2465
2466   }//cone size loop
2467
2468   //Reset original parameters for AOD analysis
2469   GetIsolationCut()->SetPtThreshold(ptthresorg);
2470   GetIsolationCut()->SetPtFraction(ptfracorg);
2471   GetIsolationCut()->SetConeSize(rorg);
2472
2473 }
2474
2475 //_____________________________________________________________
2476 void AliAnaParticleIsolation::Print(const Option_t * opt) const
2477 {
2478
2479   //Print some relevant parameters set for the analysis
2480   if(! opt)
2481     return;
2482
2483   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2484   AliAnaCaloTrackCorrBaseClass::Print(" ");
2485
2486   printf("ReMake Isolation          = %d \n",  fReMakeIC) ;
2487   printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;
2488   printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;
2489
2490   if(fMakeSeveralIC)
2491   {
2492     printf("N Cone Sizes       =     %d\n", fNCones) ;
2493     printf("Cone Sizes          =    \n") ;
2494     for(Int_t i = 0; i < fNCones; i++)
2495       printf("  %1.2f;",  fConeSizes[i]) ;
2496     printf("    \n") ;
2497
2498     printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
2499     printf(" pT thresholds         =    \n") ;
2500     for(Int_t i = 0; i < fNPtThresFrac; i++)
2501       printf("   %2.2f;",  fPtThresholds[i]) ;
2502
2503     printf("    \n") ;
2504
2505     printf(" pT fractions         =    \n") ;
2506     for(Int_t i = 0; i < fNPtThresFrac; i++)
2507       printf("   %2.2f;",  fPtFractions[i]) ;
2508
2509     printf("    \n") ;
2510
2511     printf("sum pT thresholds         =    \n") ;
2512     for(Int_t i = 0; i < fNPtThresFrac; i++)
2513       printf("   %2.2f;",  fSumPtThresholds[i]) ;
2514
2515
2516   }
2517
2518   printf("Histograms: %3.1f < pT sum < %3.1f,  Nbin = %d\n",    fHistoPtSumMin,    fHistoPtSumMax,    fHistoNPtSumBins   );
2519   printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);
2520
2521   printf("    \n") ;
2522
2523 }
2524