PWGGA/CaloTrackCorrelations/AliAnaParticleIsolation.cxx

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