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