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