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