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