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