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