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