PWG4/PartCorrDep/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 /* $Id: AliAnaParticleIsolation.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 // Class for analysis of particle isolation
  19 // Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)
  20 //
  21 //  Class created from old AliPHOSGammaJet
  22 //  (see AliRoot versions previous Release 4-09)
  23 //
  24 // -- Author: Gustavo Conesa (LNF-INFN)
  25
  26 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
  27 //////////////////////////////////////////////////////////////////////////////
  28
  29
  30 // --- ROOT system ---
  31 #include <TClonesArray.h>
  32 #include <TList.h>
  33 #include <TObjString.h>
  34 #include <TH2F.h>
  35 #include <TClass.h>
  36
  37 // --- Analysis system ---
  38 #include "AliAnaParticleIsolation.h"
  39 #include "AliCaloTrackReader.h"
  40 #include "AliIsolationCut.h"
  41 #include "AliNeutralMesonSelection.h"
  42 #include "AliAODPWG4ParticleCorrelation.h"
  43 #include "AliMCAnalysisUtils.h"
  44 #include "AliVTrack.h"
  45 #include "AliVCluster.h"
  46
  47 ClassImp(AliAnaParticleIsolation)
  48
  49 //______________________________________________________________________________
  50   AliAnaParticleIsolation::AliAnaParticleIsolation() :
  51     AliAnaPartCorrBaseClass(),        fCalorimeter(""),
  52     fReMakeIC(0),                     fMakeSeveralIC(0),
  53     // Several IC
  54     fNCones(0),                       fNPtThresFrac(0),
  55     fConeSizes(),                     fPtThresholds(),                 fPtFractions(),
  56     // Histograms
  57     fhPtIso(0),                       fhPhiIso(0),                     fhEtaIso(0),
  58     fhPtNoIso(0),                     fhPtDecayIso(0),                 fhPtDecayNoIso(0),
  59     fhConeSumPt(0),                   fhPtInCone(0),
  60     fhFRConeSumPt(0),                 fhPtInFRCone(0),
  61     // MC histograms
  62     fhPtIsoPrompt(0),                 fhPhiIsoPrompt(0),               fhEtaIsoPrompt(0),
  63     fhPtThresIsolatedPrompt(),        fhPtFracIsolatedPrompt(),        fhPtSumIsolatedPrompt(),
  64     fhPtIsoFragmentation(0),          fhPhiIsoFragmentation(0),        fhEtaIsoFragmentation(0),
  65     fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
  66     fhPtIsoPi0Decay(0),               fhPhiIsoPi0Decay(0),             fhEtaIsoPi0Decay(0),
  67     fhPtThresIsolatedPi0Decay(),      fhPtFracIsolatedPi0Decay(),      fhPtSumIsolatedPi0Decay(),
  68     fhPtIsoEtaDecay(0),               fhPhiIsoEtaDecay(0),             fhEtaIsoEtaDecay(0),
  69     fhPtThresIsolatedEtaDecay(),      fhPtFracIsolatedEtaDecay(),      fhPtSumIsolatedEtaDecay(),
  70     fhPtIsoOtherDecay(0),             fhPhiIsoOtherDecay(0),           fhEtaIsoOtherDecay(0),
  71     fhPtThresIsolatedOtherDecay(),    fhPtFracIsolatedOtherDecay(),    fhPtSumIsolatedOtherDecay(),
  72     fhPtIsoConversion(0),             fhPhiIsoConversion(0),           fhEtaIsoConversion(0),
  73     fhPtThresIsolatedConversion(),    fhPtFracIsolatedConversion(),    fhPtSumIsolatedConversion(),
  74     fhPtIsoUnknown(0),                fhPhiIsoUnknown(0),              fhEtaIsoUnknown(0),
  75     fhPtThresIsolatedUnknown(),       fhPtFracIsolatedUnknown(),       fhPtSumIsolatedUnknown(),
  76     fhPtNoIsoPi0Decay(0),             fhPtNoIsoEtaDecay(0),            fhPtNoIsoOtherDecay(0),
  77     fhPtNoIsoPrompt(0),               fhPtIsoMCPhoton(0),              fhPtNoIsoMCPhoton(0),
  78     //Histograms settings
  79     fHistoNPtSumBins(0),              fHistoPtSumMax(0.),              fHistoPtSumMin(0.),
  80     fHistoNPtInConeBins(0),           fHistoPtInConeMax(0.),           fHistoPtInConeMin(0.)
  81 {
  82   //default ctor
  83
  84   //Initialize parameters
  85   InitParameters();
  86
  87   for(Int_t i = 0; i < 5 ; i++){
  88     fConeSizes[i]      = 0 ;
  89     fhPtSumIsolated[i] = 0 ;
  90
  91     fhPtSumIsolatedPrompt[i]        = 0 ;
  92     fhPtSumIsolatedFragmentation[i] = 0 ;
  93     fhPtSumIsolatedPi0Decay[i]      = 0 ;
  94     fhPtSumIsolatedEtaDecay[i]      = 0 ;
  95     fhPtSumIsolatedOtherDecay[i]    = 0 ;
  96     fhPtSumIsolatedConversion[i]    = 0 ;
  97     fhPtSumIsolatedUnknown[i]       = 0 ;
  98
  99     for(Int_t j = 0; j < 5 ; j++){
 100       fhPtThresIsolated[i][j] = 0 ;
 101       fhPtFracIsolated[i][j]  = 0 ;
 102
 103       fhPtThresIsolatedPrompt[i][j]       = 0 ;
 104       fhPtThresIsolatedFragmentation[i][j]= 0 ;
 105       fhPtThresIsolatedPi0Decay[i][j]     = 0 ;
 106       fhPtThresIsolatedEtaDecay[i][j]     = 0 ;
 107       fhPtThresIsolatedOtherDecay[i][j]   = 0 ;
 108       fhPtThresIsolatedConversion[i][j]   = 0 ;
 109       fhPtThresIsolatedUnknown[i][j]      = 0 ;
 110
 111       fhPtFracIsolatedPrompt[i][j]        = 0 ;
 112       fhPtFracIsolatedFragmentation[i][j] = 0 ;
 113       fhPtFracIsolatedPi0Decay[i][j]      = 0 ;
 114       fhPtFracIsolatedEtaDecay[i][j]      = 0 ;
 115       fhPtFracIsolatedOtherDecay[i][j]    = 0 ;
 116       fhPtFracIsolatedConversion[i][j]    = 0 ;
 117       fhPtFracIsolatedUnknown[i][j]       = 0 ;
 118
 119     }
 120   }
 121
 122   for(Int_t i = 0; i < 5 ; i++){
 123     fPtFractions[i] = 0 ;
 124     fPtThresholds[i]= 0 ;
 125   }
 126
 127 }
 128
 129 //______________________________________________________
 130 TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
 131 {
 132         //Save parameters used for analysis
 133   TString parList ; //this will be list of parameters used for this analysis.
 134   const Int_t buffersize = 255;
 135   char onePar[buffersize] ;
 136
 137   snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
 138   parList+=onePar ;
 139   snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 140   parList+=onePar ;
 141   snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
 142   parList+=onePar ;
 143   snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
 144   parList+=onePar ;
 145
 146   if(fMakeSeveralIC){
 147     snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
 148     parList+=onePar ;
 149     snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
 150     parList+=onePar ;
 151
 152     for(Int_t icone = 0; icone < fNCones ; icone++){
 153       snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
 154       parList+=onePar ;
 155     }
 156     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
 157       snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
 158       parList+=onePar ;
 159     }
 160     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
 161       snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
 162       parList+=onePar ;
 163     }
 164   }
 165
 166   //Get parameters set in base class.
 167   parList += GetBaseParametersList() ;
 168
 169   //Get parameters set in IC class.
 170   if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
 171
 172   return new TObjString(parList) ;
 173
 174 }
 175
 176 //________________________________________________________
 177 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
 178 {
 179   // Create histograms to be saved in output file and
 180   // store them in outputContainer
 181   TList * outputContainer = new TList() ;
 182   outputContainer->SetName("IsolatedParticleHistos") ;
 183
 184   Int_t   nptbins  = GetHistoPtBins();
 185   Int_t   nphibins = GetHistoPhiBins();
 186   Int_t   netabins = GetHistoEtaBins();
 187   Float_t ptmax    = GetHistoPtMax();
 188   Float_t phimax   = GetHistoPhiMax();
 189   Float_t etamax   = GetHistoEtaMax();
 190   Float_t ptmin    = GetHistoPtMin();
 191   Float_t phimin   = GetHistoPhiMin();
 192   Float_t etamin   = GetHistoEtaMin();
 193
 194   Int_t   nptsumbins    = fHistoNPtSumBins;
 195   Float_t ptsummax      = fHistoPtSumMax;
 196   Float_t ptsummin      = fHistoPtSumMin;
 197   Int_t   nptinconebins = fHistoNPtInConeBins;
 198   Float_t ptinconemax   = fHistoPtInConeMax;
 199   Float_t ptinconemin   = fHistoPtInConeMin;
 200
 201   if(!fMakeSeveralIC){
 202
 203     fhConeSumPt  = new TH2F
 204     ("hConePtSum","#Sigma p_{T} in isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 205     fhConeSumPt->SetYTitle("#Sigma p_{T}");
 206     fhConeSumPt->SetXTitle("p_{T} (GeV/c)");
 207     outputContainer->Add(fhConeSumPt) ;
 208
 209     fhPtInCone  = new TH2F
 210     ("hPtInCone","p_{T} in isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 211     fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
 212     fhPtInCone->SetXTitle("p_{T} (GeV/c)");
 213     outputContainer->Add(fhPtInCone) ;
 214
 215     fhFRConeSumPt  = new TH2F
 216     ("hFRConePtSum","#Sigma p_{T} in the froward region isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 217     fhFRConeSumPt->SetYTitle("#Sigma p_{T}");
 218     fhFRConeSumPt->SetXTitle("p_{T} (GeV/c)");
 219     outputContainer->Add(fhFRConeSumPt) ;
 220
 221     fhPtInFRCone  = new TH2F
 222     ("hPtInFRCone","p_{T} in froward region isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 223     fhPtInFRCone->SetYTitle("p_{T in cone} (GeV/c)");
 224     fhPtInFRCone->SetXTitle("p_{T} (GeV/c)");
 225     outputContainer->Add(fhPtInFRCone) ;
 226
 227     fhPtIso  = new TH1F("hPt","Number of isolated particles",nptbins,ptmin,ptmax);
 228     fhPtIso->SetYTitle("N");
 229     fhPtIso->SetXTitle("p_{T}(GeV/c)");
 230     outputContainer->Add(fhPtIso) ;
 231
 232     fhPhiIso  = new TH2F
 233     ("hPhi","Number of isolated particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 234     fhPhiIso->SetYTitle("#phi");
 235     fhPhiIso->SetXTitle("p_{T} (GeV/c)");
 236     outputContainer->Add(fhPhiIso) ;
 237
 238     fhEtaIso  = new TH2F
 239     ("hEta","Number of isolated particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 240     fhEtaIso->SetYTitle("#eta");
 241     fhEtaIso->SetXTitle("p_{T} (GeV/c)");
 242     outputContainer->Add(fhEtaIso) ;
 243
 244     fhPtNoIso  = new TH1F("hPtNoIso","Number of not isolated leading particles",nptbins,ptmin,ptmax);
 245     fhPtNoIso->SetYTitle("N");
 246     fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
 247     outputContainer->Add(fhPtNoIso) ;
 248
 249     fhPtDecayIso  = new TH1F("hPtDecayIso","Number of isolated #pi^{0} decay particles",nptbins,ptmin,ptmax);
 250     fhPtNoIso->SetYTitle("N");
 251     fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
 252     outputContainer->Add(fhPtDecayIso) ;
 253
 254     fhPtDecayNoIso  = new TH1F("hPtDecayNoIso","Number of not isolated leading pi0 decay particles",nptbins,ptmin,ptmax);
 255     fhPtNoIso->SetYTitle("N");
 256     fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
 257     outputContainer->Add(fhPtDecayNoIso) ;
 258
 259     if(IsDataMC()){
 260
 261       fhPtIsoPrompt  = new TH1F("hPtMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax);
 262       fhPtIsoPrompt->SetYTitle("N");
 263       fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
 264       outputContainer->Add(fhPtIsoPrompt) ;
 265
 266       fhPhiIsoPrompt  = new TH2F
 267       ("hPhiMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 268       fhPhiIsoPrompt->SetYTitle("#phi");
 269       fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 270       outputContainer->Add(fhPhiIsoPrompt) ;
 271
 272       fhEtaIsoPrompt  = new TH2F
 273       ("hEtaMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 274       fhEtaIsoPrompt->SetYTitle("#eta");
 275       fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 276       outputContainer->Add(fhEtaIsoPrompt) ;
 277
 278       fhPtIsoFragmentation  = new TH1F("hPtMCFragmentation","Number of isolated #gamma",nptbins,ptmin,ptmax);
 279       fhPtIsoFragmentation->SetYTitle("N");
 280       fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
 281       outputContainer->Add(fhPtIsoFragmentation) ;
 282
 283       fhPhiIsoFragmentation  = new TH2F
 284       ("hPhiMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 285       fhPhiIsoFragmentation->SetYTitle("#phi");
 286       fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 287       outputContainer->Add(fhPhiIsoFragmentation) ;
 288
 289       fhEtaIsoFragmentation  = new TH2F
 290       ("hEtaMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 291       fhEtaIsoFragmentation->SetYTitle("#eta");
 292       fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 293       outputContainer->Add(fhEtaIsoFragmentation) ;
 294
 295       fhPtIsoPi0Decay  = new TH1F("hPtMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
 296       fhPtIsoPi0Decay->SetYTitle("N");
 297       fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
 298       outputContainer->Add(fhPtIsoPi0Decay) ;
 299
 300       fhPhiIsoPi0Decay  = new TH2F
 301       ("hPhiMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 302       fhPhiIsoPi0Decay->SetYTitle("#phi");
 303       fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 304       outputContainer->Add(fhPhiIsoPi0Decay) ;
 305
 306       fhEtaIsoPi0Decay  = new TH2F
 307       ("hEtaMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 308       fhEtaIsoPi0Decay->SetYTitle("#eta");
 309       fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 310       outputContainer->Add(fhEtaIsoPi0Decay) ;
 311
 312       fhPtIsoEtaDecay  = new TH1F("hPtMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax);
 313       fhPtIsoEtaDecay->SetYTitle("N");
 314       fhPtIsoEtaDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 315       outputContainer->Add(fhPtIsoEtaDecay) ;
 316
 317       fhPhiIsoEtaDecay  = new TH2F
 318       ("hPhiMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 319       fhPhiIsoEtaDecay->SetYTitle("#phi");
 320       fhPhiIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 321       outputContainer->Add(fhPhiIsoEtaDecay) ;
 322
 323       fhEtaIsoEtaDecay  = new TH2F
 324       ("hEtaMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 325       fhEtaIsoEtaDecay->SetYTitle("#eta");
 326       fhEtaIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 327       outputContainer->Add(fhEtaIsoEtaDecay) ;
 328
 329       fhPtIsoOtherDecay  = new TH1F("hPtMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax);
 330       fhPtIsoOtherDecay->SetYTitle("N");
 331       fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 332       outputContainer->Add(fhPtIsoOtherDecay) ;
 333
 334       fhPhiIsoOtherDecay  = new TH2F
 335       ("hPhiMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 336       fhPhiIsoOtherDecay->SetYTitle("#phi");
 337       fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 338       outputContainer->Add(fhPhiIsoOtherDecay) ;
 339
 340       fhEtaIsoOtherDecay  = new TH2F
 341       ("hEtaMCOtherDecay","Number of isolated #gamma non-#pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 342       fhEtaIsoOtherDecay->SetYTitle("#eta");
 343       fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 344       outputContainer->Add(fhEtaIsoOtherDecay) ;
 345
 346       fhPtIsoConversion  = new TH1F("hPtMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax);
 347       fhPtIsoConversion->SetYTitle("N");
 348       fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
 349       outputContainer->Add(fhPtIsoConversion) ;
 350
 351       fhPhiIsoConversion  = new TH2F
 352       ("hPhiMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 353       fhPhiIsoConversion->SetYTitle("#phi");
 354       fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 355       outputContainer->Add(fhPhiIsoConversion) ;
 356
 357       fhEtaIsoConversion  = new TH2F
 358       ("hEtaMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 359       fhEtaIsoConversion->SetYTitle("#eta");
 360       fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 361       outputContainer->Add(fhEtaIsoConversion) ;
 362
 363       fhPtIsoUnknown  = new TH1F("hPtMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax);
 364       fhPtIsoUnknown->SetYTitle("N");
 365       fhPtIsoUnknown->SetXTitle("p_{T}(GeV/c)");
 366       outputContainer->Add(fhPtIsoUnknown) ;
 367
 368       fhPhiIsoUnknown  = new TH2F
 369       ("hPhiMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 370       fhPhiIsoUnknown->SetYTitle("#phi");
 371       fhPhiIsoUnknown->SetXTitle("p_{T} (GeV/c)");
 372       outputContainer->Add(fhPhiIsoUnknown) ;
 373
 374       fhEtaIsoUnknown  = new TH2F
 375       ("hEtaMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 376       fhEtaIsoUnknown->SetYTitle("#eta");
 377       fhEtaIsoUnknown->SetXTitle("p_{T} (GeV/c)");
 378       outputContainer->Add(fhEtaIsoUnknown) ;
 379
 380       fhPtNoIsoPi0Decay  = new TH1F
 381       ("hPtNoIsoPi0Decay","Number of not isolated leading #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
 382       fhPtNoIsoPi0Decay->SetYTitle("N");
 383       fhPtNoIsoPi0Decay->SetXTitle("p_{T} (GeV/c)");
 384       outputContainer->Add(fhPtNoIsoPi0Decay) ;
 385
 386       fhPtNoIsoEtaDecay  = new TH1F
 387       ("hPtNoIsoEtaDecay","Number of not isolated leading #gamma from eta decay",nptbins,ptmin,ptmax);
 388       fhPtNoIsoEtaDecay->SetYTitle("N");
 389       fhPtNoIsoEtaDecay->SetXTitle("p_{T} (GeV/c)");
 390       outputContainer->Add(fhPtNoIsoEtaDecay) ;
 391
 392       fhPtNoIsoOtherDecay  = new TH1F
 393       ("hPtNoIsoOtherDecay","Number of not isolated leading #gamma from other decay",nptbins,ptmin,ptmax);
 394       fhPtNoIsoOtherDecay->SetYTitle("N");
 395       fhPtNoIsoOtherDecay->SetXTitle("p_{T} (GeV/c)");
 396       outputContainer->Add(fhPtNoIsoOtherDecay) ;
 397
 398       fhPtNoIsoPrompt  = new TH1F
 399       ("hPtNoIsoPrompt","Number of not isolated leading prompt #gamma",nptbins,ptmin,ptmax);
 400       fhPtNoIsoPrompt->SetYTitle("N");
 401       fhPtNoIsoPrompt->SetXTitle("p_{T} (GeV/c)");
 402       outputContainer->Add(fhPtNoIsoPrompt) ;
 403
 404       fhPtIsoMCPhoton  = new TH1F
 405       ("hPtIsoMCPhoton","Number of isolated leading  #gamma",nptbins,ptmin,ptmax);
 406       fhPtIsoMCPhoton->SetYTitle("N");
 407       fhPtIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
 408       outputContainer->Add(fhPtIsoMCPhoton) ;
 409
 410       fhPtNoIsoMCPhoton  = new TH1F
 411       ("hPtNoIsoMCPhoton","Number of not isolated leading #gamma",nptbins,ptmin,ptmax);
 412       fhPtNoIsoMCPhoton->SetYTitle("N");
 413       fhPtNoIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
 414       outputContainer->Add(fhPtNoIsoMCPhoton) ;
 415     }//Histos with MC
 416
 417   }
 418
 419   if(fMakeSeveralIC){
 420     const Int_t buffersize = 255;
 421                 char name[buffersize];
 422                 char title[buffersize];
 423                 for(Int_t icone = 0; icone<fNCones; icone++){
 424                   snprintf(name, buffersize,"hPtSum_Cone_%d",icone);
 425                   snprintf(title, buffersize,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 426                   fhPtSumIsolated[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 427                   fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 428                   fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
 429                   outputContainer->Add(fhPtSumIsolated[icone]) ;
 430
 431                   if(IsDataMC()){
 432                     snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);
 433                     snprintf(title, buffersize,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 434                     fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 435                     fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 436                     fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
 437                     outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
 438
 439                     snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);
 440                     snprintf(title, buffersize,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 441                     fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 442                     fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 443                     fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
 444                     outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
 445
 446                     snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);
 447                     snprintf(title, buffersize,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 448                     fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 449                     fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 450                     fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
 451                     outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
 452
 453         snprintf(name, buffersize,"hPtSumEtaDecay_Cone_%d",icone);
 454                     snprintf(title, buffersize,"Candidate EtaDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 455                     fhPtSumIsolatedEtaDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 456                     fhPtSumIsolatedEtaDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 457                     fhPtSumIsolatedEtaDecay[icone]->SetXTitle("p_{T} (GeV/c)");
 458                     outputContainer->Add(fhPtSumIsolatedEtaDecay[icone]) ;
 459
 460                     snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);
 461                     snprintf(title, buffersize,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 462                     fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 463                     fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 464                     fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
 465                     outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
 466
 467                     snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);
 468                     snprintf(title, buffersize,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 469                     fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 470                     fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 471                     fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
 472                     outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
 473
 474                     snprintf(name, buffersize,"hPtSumUnknown_Cone_%d",icone);
 475                     snprintf(title, buffersize,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 476                     fhPtSumIsolatedUnknown[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 477                     fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 478                     fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
 479                     outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
 480
 481                   }//Histos with MC
 482
 483                   for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){
 484                     snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
 485                     snprintf(title, buffersize,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 486                     fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 487                     fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 488                     outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
 489
 490                     snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
 491                     snprintf(title, buffersize,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 492                     fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 493                     fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 494                     outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
 495
 496                     if(IsDataMC()){
 497                       snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);
 498                       snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 499                       fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 500                       fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 501                       outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
 502
 503                       snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);
 504                       snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 505                       fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 506                       fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 507                       outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
 508
 509                       snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);
 510                       snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 511                       fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 512                       fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 513                       outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
 514
 515                       snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);
 516                       snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 517                       fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 518                       fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 519                       outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
 520
 521                       snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
 522                       snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 523                       fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 524                       fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 525                       outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
 526
 527                       snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
 528                       snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 529                       fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 530                       fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 531                       outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
 532
 533           snprintf(name, buffersize,"hPtThresMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
 534                       snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 535                       fhPtThresIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 536                       fhPtThresIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 537                       outputContainer->Add(fhPtThresIsolatedEtaDecay[icone][ipt]) ;
 538
 539                       snprintf(name, buffersize,"hPtFracMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
 540                       snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 541                       fhPtFracIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 542                       fhPtFracIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 543                       outputContainer->Add(fhPtFracIsolatedEtaDecay[icone][ipt]) ;
 544
 545
 546                       snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
 547                       snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 548                       fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 549                       fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 550                       outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
 551
 552                       snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
 553                       snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 554                       fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 555                       fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 556                       outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
 557
 558                       snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);
 559                       snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 560                       fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 561                       fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 562                       outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
 563
 564                       snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);
 565                       snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 566                       fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 567                       fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 568                       outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
 569
 570                       snprintf(name, buffersize,"hPtThresMCUnknown_Cone_%d_Pt%d",icone,ipt);
 571                       snprintf(title, buffersize,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 572                       fhPtThresIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 573                       fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 574                       outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
 575
 576                       snprintf(name, buffersize,"hPtFracMCUnknown_Cone_%d_Pt%d",icone,ipt);
 577                       snprintf(title, buffersize,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 578                       fhPtFracIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 579                       fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 580                       outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
 581
 582                     }//Histos with MC
 583
 584                   }//icone loop
 585                 }//ipt loop
 586   }
 587
 588
 589   return outputContainer ;
 590
 591 }
 592
 593 //____________________________________________
 594 void AliAnaParticleIsolation::InitParameters()
 595 {
 596
 597   //Initialize the parameters of the analysis.
 598   SetInputAODName("PWG4Particle");
 599   SetAODObjArrayName("IsolationCone");
 600   AddToHistogramsName("AnaIsolation_");
 601
 602   fCalorimeter = "PHOS" ;
 603   fReMakeIC = kFALSE ;
 604   fMakeSeveralIC = kFALSE ;
 605
 606   //----------- Several IC-----------------
 607   fNCones          = 5 ;
 608   fNPtThresFrac    = 5 ;
 609   fConeSizes[0]    = 0.1;  fConeSizes[1]     = 0.2;  fConeSizes[2]    = 0.3; fConeSizes[3]    = 0.4;  fConeSizes[4]    = 0.5;
 610   fPtThresholds[0] = 1.;   fPtThresholds[1] = 2.;    fPtThresholds[2] = 3.;  fPtThresholds[3] = 4.;   fPtThresholds[4] = 5.;
 611   fPtFractions[0]  = 0.05; fPtFractions[1]  = 0.075; fPtFractions[2]  = 0.1; fPtFractions[3]  = 1.25; fPtFractions[4]  = 1.5;
 612
 613   //------------- Histograms settings -------
 614   fHistoNPtSumBins = 100 ;
 615   fHistoPtSumMax   = 50 ;
 616   fHistoPtSumMin   = 0.  ;
 617
 618   fHistoNPtInConeBins = 100 ;
 619   fHistoPtInConeMax   = 50 ;
 620   fHistoPtInConeMin   = 0.  ;
 621
 622 }
 623
 624 //__________________________________________________
 625 void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
 626 {
 627   //Do analysis and fill aods
 628   //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
 629
 630   if(!GetInputAODBranch()){
 631     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 632     abort();
 633   }
 634
 635   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
 636     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());
 637     abort();
 638   }
 639
 640   Int_t n = 0, nfrac = 0;
 641   Bool_t  isolated  = kFALSE ;
 642   Float_t coneptsum = 0 ;
 643   TObjArray * pl    = 0x0; ;
 644
 645   //Select the calorimeter for candidate isolation with neutral particles
 646   if(fCalorimeter == "PHOS")
 647     pl = GetPHOSClusters();
 648   else if (fCalorimeter == "EMCAL")
 649     pl = GetEMCALClusters();
 650
 651   //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
 652   Double_t ptLeading = 0. ;
 653   Int_t    idLeading = -1 ;
 654   TLorentzVector mom ;
 655   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 656   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
 657
 658   for(Int_t iaod = 0; iaod < naod; iaod++){
 659     AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 660
 661     //If too small or too large pt, skip
 662     if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ;
 663
 664     //check if it is low pt trigger particle, then adjust the isolation method
 665     if(aodinput->Pt() < GetIsolationCut()->GetPtThreshold() || aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold())
 666       continue ; //trigger should not come from underlying event
 667
 668     //vertex cut in case of mixing
 669     if (GetMixedEvent()) {
 670       Int_t evt=-1;
 671       Int_t id =-1;
 672       if     (aodinput->GetCaloLabel(0)  >= 0 ){
 673         id=aodinput->GetCaloLabel(0);
 674         if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
 675       }
 676       else if(aodinput->GetTrackLabel(0) >= 0 ){
 677         id=aodinput->GetTrackLabel(0);
 678         if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
 679       }
 680       else continue;
 681
 682       if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
 683         return ;
 684     }
 685
 686     //find the leading particles with highest momentum
 687     if ((aodinput->Pt())>ptLeading) {
 688       ptLeading = aodinput->Pt() ;
 689       idLeading = iaod ;
 690     }
 691     aodinput->SetLeadingParticle(kFALSE);
 692   }//finish searching for leading trigger particle
 693
 694   // Check isolation of leading particle
 695   if(idLeading < 0) return;
 696
 697   AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(idLeading));
 698   aodinput->SetLeadingParticle(kTRUE);
 699
 700   //After cuts, study isolation
 701   n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
 702   GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,GetReader(), kTRUE, aodinput, GetAODObjArrayName(), n,nfrac,coneptsum, isolated);
 703   aodinput->SetIsolated(isolated);
 704
 705   if(GetDebug() > 1) {
 706     if(isolated)printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
 707     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
 708   }
 709
 710 }
 711
 712 //_________________________________________________________
 713 void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
 714 {
 715   //Do analysis and fill histograms
 716   Int_t   n = 0, nfrac = 0;
 717   Bool_t  isolated  = kFALSE ;
 718   Float_t coneptsum = 0 ;
 719
 720   //Loop on stored AOD
 721   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 722   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
 723
 724   //Get vertex for photon momentum calculation
 725   Double_t vertex[]={0,0,0} ; //vertex ;
 726   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
 727           GetReader()->GetVertex(vertex);
 728   }
 729
 730   for(Int_t iaod = 0; iaod < naod ; iaod++){
 731     AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 732
 733     if(!aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
 734
 735     Bool_t  isolation  = aod->IsIsolated();
 736     Bool_t  decay      = aod->IsTagged();
 737     Float_t pt         = aod->Pt();
 738     Float_t phi        = aod->Phi();
 739     Float_t eta        = aod->Eta();
 740     Float_t conesize   = GetIsolationCut()->GetConeSize();
 741
 742     //Recover reference arrays with clusters and tracks
 743     TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
 744     TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
 745     //If too small or too large pt, skip
 746     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
 747
 748     // --- In case of redoing isolation from delta AOD ----
 749     if(fMakeSeveralIC) {
 750       //Analysis of multiple IC at same time
 751       MakeSeveralICAnalysis(aod);
 752     }
 753     else if(fReMakeIC){
 754       //In case a more strict IC is needed in the produced AOD
 755       n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
 756       GetIsolationCut()->MakeIsolationCut(reftracks, refclusters, GetReader(), kFALSE, aod, "", n,nfrac,coneptsum, isolated);
 757       fhConeSumPt->Fill(pt,coneptsum);
 758       if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
 759     }
 760     // --- -------------------------------------------- ----
 761
 762     //Fill pt distribution of particles in cone
 763     //Tracks
 764     coneptsum = 0;
 765     Double_t sumptFR = 0. ;
 766     TObjArray * trackList   = GetCTSTracks() ;
 767     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++){
 768       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
 769       //fill the histograms at forward range
 770       if(!track){
 771         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
 772         continue;
 773       }
 774       Double_t dPhi = phi - track->Phi() + TMath::PiOver2();
 775       Double_t dEta = eta - track->Eta();
 776       Double_t arg  = dPhi*dPhi + dEta*dEta;
 777       if(TMath::Sqrt(arg) < conesize){
 778         fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
 779         sumptFR+=track->Pt();
 780       }
 781       dPhi = phi - track->Phi() - TMath::PiOver2();
 782       arg  = dPhi*dPhi + dEta*dEta;
 783       if(TMath::Sqrt(arg) < conesize){
 784         fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
 785         sumptFR+=track->Pt();
 786       }
 787     }
 788     fhFRConeSumPt->Fill(pt,sumptFR);
 789     if(reftracks){
 790       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++){
 791         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
 792         fhPtInCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
 793         coneptsum+=track->Pt();
 794       }
 795     }
 796
 797     //CaloClusters
 798     if(refclusters){
 799       TLorentzVector mom ;
 800       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++){
 801         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
 802         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
 803
 804         fhPtInCone->Fill(pt, mom.Pt());
 805         coneptsum+=mom.Pt();
 806       }
 807     }
 808
 809     if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);
 810
 811     if(!fReMakeIC) fhConeSumPt->Fill(pt,coneptsum);
 812
 813     if(isolation){
 814
 815       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED, fill histograms\n", iaod);
 816
 817       fhPtIso  ->Fill(pt);
 818       fhPhiIso ->Fill(pt,phi);
 819       fhEtaIso ->Fill(pt,eta);
 820       if (decay) fhPtDecayIso->Fill(pt);
 821
 822       if(IsDataMC()){
 823         Int_t tag =aod->GetTag();
 824
 825         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 826         {
 827           fhPtIsoMCPhoton  ->Fill(pt);
 828         }
 829
 830         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
 831           fhPtIsoPrompt  ->Fill(pt);
 832           fhPhiIsoPrompt ->Fill(pt,phi);
 833           fhEtaIsoPrompt ->Fill(pt,eta);
 834         }
 835         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
 836         {
 837           fhPtIsoFragmentation  ->Fill(pt);
 838           fhPhiIsoFragmentation ->Fill(pt,phi);
 839           fhEtaIsoFragmentation ->Fill(pt,eta);
 840         }
 841         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
 842         {
 843           fhPtIsoPi0Decay  ->Fill(pt);
 844           fhPhiIsoPi0Decay ->Fill(pt,phi);
 845           fhEtaIsoPi0Decay ->Fill(pt,eta);
 846         }
 847         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
 848         {
 849           fhPtIsoEtaDecay  ->Fill(pt);
 850           fhPhiIsoEtaDecay ->Fill(pt,phi);
 851           fhEtaIsoEtaDecay ->Fill(pt,eta);
 852         }
 853         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 854         {
 855           fhPtIsoOtherDecay  ->Fill(pt);
 856           fhPhiIsoOtherDecay ->Fill(pt,phi);
 857           fhEtaIsoOtherDecay ->Fill(pt,eta);
 858         }
 859         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
 860         {
 861           fhPtIsoConversion  ->Fill(pt);
 862           fhPhiIsoConversion ->Fill(pt,phi);
 863           fhEtaIsoConversion ->Fill(pt,eta);
 864         }
 865         else // anything else
 866         {
 867           fhPtIsoUnknown  ->Fill(pt);
 868           fhPhiIsoUnknown ->Fill(pt,phi);
 869           fhEtaIsoUnknown ->Fill(pt,eta);
 870         }
 871       }//Histograms with MC
 872
 873     }//Isolated histograms
 874
 875     if(!isolation)
 876     {
 877       fhPtNoIso  ->Fill(pt);
 878       if (decay) fhPtDecayNoIso->Fill(pt);
 879
 880       if(IsDataMC()){
 881         Int_t tag =aod->GetTag();
 882
 883         if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 884         {
 885           fhPtNoIsoMCPhoton->Fill(pt);
 886         }
 887
 888         if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
 889         {
 890           fhPtNoIsoPi0Decay->Fill(pt);
 891         }
 892         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
 893         {
 894           fhPtNoIsoEtaDecay->Fill(pt);
 895         }
 896         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 897         {
 898           fhPtNoIsoOtherDecay->Fill(pt);
 899         }
 900         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
 901         {
 902           fhPtNoIsoPrompt->Fill(pt);
 903         }
 904
 905       }
 906     }
 907
 908   }// aod loop
 909
 910 }
 911
 912
 913 //_____________________________________________________________________________________
 914 void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph)
 915 {
 916   //Isolation Cut Analysis for both methods and different pt cuts and cones
 917   Float_t ptC = ph->Pt();
 918   Int_t tag   = ph->GetTag();
 919
 920   //Keep original setting used when filling AODs, reset at end of analysis
 921   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
 922   Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();
 923   Float_t rorg       = GetIsolationCut()->GetConeSize();
 924
 925   Float_t coneptsum = 0 ;
 926   Int_t n[10][10];//[fNCones][fNPtThresFrac];
 927   Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
 928   Bool_t  isolated   = kFALSE;
 929
 930   //Loop on cone sizes
 931   for(Int_t icone = 0; icone<fNCones; icone++){
 932     GetIsolationCut()->SetConeSize(fConeSizes[icone]);
 933     coneptsum = 0 ;
 934
 935     //Loop on ptthresholds
 936     for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
 937       n[icone][ipt]=0;
 938       nfrac[icone][ipt]=0;
 939       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
 940       GetIsolationCut()->MakeIsolationCut(ph->GetObjArray(GetAODObjArrayName()+"Tracks"),
 941                                           ph->GetObjArray(GetAODObjArrayName()+"Clusters"),
 942                                           GetReader(), kFALSE, ph, "",n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
 943
 944       //Normal ptThreshold cut
 945       if(n[icone][ipt] == 0) {
 946         fhPtThresIsolated[icone][ipt]->Fill(ptC);
 947         if(IsDataMC()){
 948           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtThresIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
 949           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtThresIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
 950           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
 951           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtThresIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
 952           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtThresIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
 953           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtThresIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
 954           else  fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC) ;
 955         }
 956       }
 957
 958       //Pt threshold on pt cand/ pt in cone fraction
 959       if(nfrac[icone][ipt] == 0) {
 960         fhPtFracIsolated[icone][ipt]->Fill(ptC);
 961         if(IsDataMC()){
 962           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtFracIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
 963           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtFracIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
 964           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
 965           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtFracIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
 966           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtFracIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
 967           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtFracIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
 968           else  fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ;
 969         }
 970       }
 971     }//pt thresh loop
 972
 973     //Sum in cone histograms
 974     fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
 975     if(IsDataMC()){
 976       if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtSumIsolatedPrompt[icone]       ->Fill(ptC,coneptsum) ;
 977       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtSumIsolatedConversion[icone]   ->Fill(ptC,coneptsum) ;
 978       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
 979       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtSumIsolatedPi0Decay[icone]     ->Fill(ptC,coneptsum) ;
 980       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtSumIsolatedEtaDecay[icone]     ->Fill(ptC,coneptsum) ;
 981       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtSumIsolatedOtherDecay[icone]   ->Fill(ptC,coneptsum) ;
 982       else  fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
 983     }
 984
 985   }//cone size loop
 986
 987   //Reset original parameters for AOD analysis
 988   GetIsolationCut()->SetPtThreshold(ptthresorg);
 989   GetIsolationCut()->SetPtFraction(ptfracorg);
 990   GetIsolationCut()->SetConeSize(rorg);
 991
 992 }
 993
 994 //_____________________________________________________________
 995 void AliAnaParticleIsolation::Print(const Option_t * opt) const
 996 {
 997
 998   //Print some relevant parameters set for the analysis
 999   if(! opt)
1000     return;
1001
1002   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1003   AliAnaPartCorrBaseClass::Print(" ");
1004
1005   printf("ReMake Isolation          = %d \n",  fReMakeIC) ;
1006   printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;
1007   printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;
1008
1009   if(fMakeSeveralIC){
1010     printf("N Cone Sizes       =     %d\n", fNCones) ;
1011     printf("Cone Sizes          =    \n") ;
1012     for(Int_t i = 0; i < fNCones; i++)
1013       printf("  %1.2f;",  fConeSizes[i]) ;
1014     printf("    \n") ;
1015
1016     printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
1017     printf(" pT thresholds         =    \n") ;
1018     for(Int_t i = 0; i < fNPtThresFrac; i++)
1019       printf("   %2.2f;",  fPtThresholds[i]) ;
1020
1021     printf("    \n") ;
1022
1023     printf(" pT fractions         =    \n") ;
1024     for(Int_t i = 0; i < fNPtThresFrac; i++)
1025       printf("   %2.2f;",  fPtFractions[i]) ;
1026
1027   }
1028
1029   printf("Histograms: %3.1f < pT sum < %3.1f,  Nbin = %d\n", fHistoPtSumMin,  fHistoPtSumMax,  fHistoNPtSumBins);
1030   printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);
1031
1032   printf("    \n") ;
1033
1034 }
1035