PWGGA/CaloTrackCorrelations/AliAnaParticleIsolation.cxx

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