PWG4/PartCorrDep/AliAnaParticleIsolation.cxx

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