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