]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrDep/AliAnaParticleIsolation.cxx
Several changes:
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleIsolation.cxx
index 3e937102b3234635e0ff49264bb33a361988787e..0572f9b5c6e91159e7ff42524f88312fe3657260 100755 (executable)
 /**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes hereby granted      *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-/* $Id: AliAnaParticleIsolation.cxx 28688 2008-09-11 15:04:07Z gconesab $ */\r
-\r
-//_________________________________________________________________________\r
-// Class for analysis of particle isolation\r
-// Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)\r
-//\r
-//  Class created from old AliPHOSGammaJet \r
-//  (see AliRoot versions previous Release 4-09)\r
-//\r
-// -- Author: Gustavo Conesa (LNF-INFN) \r
-//////////////////////////////////////////////////////////////////////////////\r
-  \r
-  \r
-// --- ROOT system --- \r
-#include <TList.h>\r
-\r
-// --- Analysis system --- \r
-#include "AliAnaParticleIsolation.h" \r
-#include "AliLog.h"\r
-#include "AliCaloTrackReader.h"\r
-#include "AliIsolationCut.h"\r
-#include "AliNeutralMesonSelection.h"\r
-#include "AliAODPWG4ParticleCorrelation.h"\r
-#include "AliMCAnalysisUtils.h"\r
-\r
-ClassImp(AliAnaParticleIsolation)\r
-  \r
-//____________________________________________________________________________\r
-  AliAnaParticleIsolation::AliAnaParticleIsolation() : \r
-    AliAnaPartCorrBaseClass(), fCalorimeter(""), fVertex(),\r
-       fReMakeIC(0), fMakeSeveralIC(0), fMakeInvMass(0),\r
-    fhPtIso(0),fhPhiIso(0),fhEtaIso(0), fhConeSumPt(0), fhPtInCone(0),\r
-    //Several IC\r
-    fNCones(0),fNPtThresFrac(0), fConeSizes(),  fPtThresholds(),  fPtFractions(), \r
-    //MC\r
-       fhPtIsoPrompt(0),fhPhiIsoPrompt(0),fhEtaIsoPrompt(0), \r
-    fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),\r
-       fhPtIsoFragmentation(0),fhPhiIsoFragmentation(0),fhEtaIsoFragmentation(0), \r
-    fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),\r
-       fhPtIsoPi0Decay(0),fhPhiIsoPi0Decay(0),fhEtaIsoPi0Decay(0),\r
-    fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),\r
-       fhPtIsoOtherDecay(0),fhPhiIsoOtherDecay(0),fhEtaIsoOtherDecay(0), \r
-    fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),\r
-       fhPtIsoConversion(0),fhPhiIsoConversion(0),fhEtaIsoConversion(0), \r
-    fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),\r
-       fhPtIsoUnknown(0),fhPhiIsoUnknown(0),fhEtaIsoUnknown(0), \r
-    fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown(),\r
-       //Histograms settings\r
-       fHistoNPtSumBins(0),    fHistoPtSumMax(0.),    fHistoPtSumMin(0.),\r
-       fHistoNPtInConeBins(0), fHistoPtInConeMax(0.), fHistoPtInConeMin(0.)\r
-{\r
-  //default ctor\r
-  \r
-  //Initialize parameters\r
-  InitParameters();\r
-  \r
-  fVertex[0] = 0.;  fVertex[1] = 0.;  fVertex[2] = 0.;  \r
-         \r
-  for(Int_t i = 0; i < 5 ; i++){ \r
-    fConeSizes[i] = 0 ; \r
-    fhPtSumIsolated[i] = 0 ;  \r
-    \r
-    fhPtSumIsolatedPrompt[i] = 0 ;  \r
-    fhPtSumIsolatedFragmentation[i] = 0 ;  \r
-    fhPtSumIsolatedPi0Decay[i] = 0 ;  \r
-    fhPtSumIsolatedOtherDecay[i] = 0 ;  \r
-    fhPtSumIsolatedConversion[i] = 0 ;  \r
-    fhPtSumIsolatedUnknown[i] = 0 ;  \r
-    \r
-    for(Int_t j = 0; j < 5 ; j++){ \r
-      fhPtThresIsolated[i][j] = 0 ;  \r
-      fhPtFracIsolated[i][j] = 0 ; \r
-      \r
-      fhPtThresIsolatedPrompt[i][j] = 0 ;  \r
-      fhPtThresIsolatedFragmentation[i][j] = 0 ; \r
-      fhPtThresIsolatedPi0Decay[i][j] = 0 ;  \r
-      fhPtThresIsolatedOtherDecay[i][j] = 0 ;  \r
-      fhPtThresIsolatedConversion[i][j] = 0 ;  \r
-      fhPtThresIsolatedUnknown[i][j] = 0 ;  \r
-  \r
-      fhPtFracIsolatedPrompt[i][j] = 0 ;  \r
-      fhPtFracIsolatedFragmentation[i][j] = 0 ;  \r
-      fhPtFracIsolatedPi0Decay[i][j] = 0 ;  \r
-      fhPtFracIsolatedOtherDecay[i][j] = 0 ;  \r
-      fhPtFracIsolatedConversion[i][j] = 0 ;\r
-      fhPtFracIsolatedUnknown[i][j] = 0 ;  \r
\r
-    }  \r
-  } \r
-  \r
-  for(Int_t i = 0; i < 5 ; i++){ \r
-    fPtFractions[i]=  0 ; \r
-    fPtThresholds[i]= 0 ; \r
-  } \r
-\r
-\r
-}\r
-\r
-//____________________________________________________________________________\r
-AliAnaParticleIsolation::AliAnaParticleIsolation(const AliAnaParticleIsolation & g) : \r
-  AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),fVertex(),\r
-  fReMakeIC(g.fReMakeIC), fMakeSeveralIC(g.fMakeSeveralIC),  fMakeInvMass(g.fMakeInvMass),\r
-  fhPtIso(g.fhPtIso),fhPhiIso(g.fhPhiIso),fhEtaIso(g.fhEtaIso), \r
-  fhConeSumPt(g.fhConeSumPt), fhPtInCone(g.fhPtInCone),\r
-  //Several IC\r
-  fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(),  fPtFractions(), \r
-  fhPtThresIsolated(),  fhPtFracIsolated(), fhPtSumIsolated(),\r
-  //MC\r
-  fhPtIsoPrompt(g.fhPtIsoPrompt),fhPhiIsoPrompt(g.fhPhiIsoPrompt),fhEtaIsoPrompt(g.fhEtaIsoPrompt), \r
-  fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),\r
-  fhPtIsoFragmentation(g.fhPtIsoFragmentation),fhPhiIsoFragmentation(g.fhPhiIsoFragmentation),fhEtaIsoFragmentation(g.fhEtaIsoFragmentation),\r
-  fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),\r
-  fhPtIsoPi0Decay(g.fhPtIsoPi0Decay),fhPhiIsoPi0Decay(g.fhPhiIsoPi0Decay),fhEtaIsoPi0Decay(g.fhEtaIsoPi0Decay), \r
-  fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),\r
-  fhPtIsoOtherDecay(g.fhPtIsoOtherDecay),fhPhiIsoOtherDecay(g.fhPhiIsoOtherDecay),fhEtaIsoOtherDecay(g.fhEtaIsoOtherDecay), \r
-  fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),\r
-  fhPtIsoConversion(g. fhPtIsoConversion),fhPhiIsoConversion(g.fhPhiIsoConversion),fhEtaIsoConversion(g.fhEtaIsoConversion), \r
-  fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),\r
-  fhPtIsoUnknown(g.fhPtIsoUnknown),fhPhiIsoUnknown(g.fhPhiIsoUnknown),fhEtaIsoUnknown(g.fhEtaIsoUnknown), \r
-  fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown(),\r
-  //Histograms\r
-  fHistoNPtSumBins(g.fHistoNPtSumBins), fHistoPtSumMax(g.fHistoPtSumMax), fHistoPtSumMin(g.fHistoPtSumMax),\r
-  fHistoNPtInConeBins(g.fHistoNPtInConeBins), fHistoPtInConeMax(g.fHistoPtInConeMax), fHistoPtInConeMin(g.fHistoPtInConeMin)\r
-{  \r
-  // cpy ctor\r
-  \r
-       fVertex[0] = g.fVertex[0];  \r
-       fVertex[1] = g.fVertex[1];  \r
-       fVertex[2] = g.fVertex[2];  \r
-\r
-       //Several IC\r
-       for(Int_t i = 0; i < fNCones ; i++){\r
-               fConeSizes[i] =  g.fConeSizes[i];\r
-               fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; \r
-               \r
-               fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i]; \r
-               fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i]; \r
-               fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i]; \r
-               fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i]; \r
-               fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i]; \r
-               fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i]; \r
-               \r
-               for(Int_t j = 0; j < fNPtThresFrac ; j++){\r
-                       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j]; \r
-                       fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];\r
-                       \r
-                       fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j]; \r
-                       fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j]; \r
-                       fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j]; \r
-                       fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j]; \r
-                       fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j]; \r
-                       fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j]; \r
-                       \r
-                       fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j]; \r
-                       fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j]; \r
-                       fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j]; \r
-                       fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j]; \r
-                       fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j]; \r
-                       fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j]; \r
-                       \r
-               } \r
-       }\r
-       \r
-       for(Int_t i = 0; i < fNPtThresFrac ; i++){\r
-               fPtFractions[i]=   g.fPtFractions[i];\r
-               fPtThresholds[i]=   g.fPtThresholds[i];\r
-       }\r
-       \r
-}\r
-\r
-//_________________________________________________________________________\r
-AliAnaParticleIsolation & AliAnaParticleIsolation::operator = (const AliAnaParticleIsolation & g)\r
-{\r
-       // assignment operator\r
-  \r
-       if(&g == this) return *this;\r
-\r
-       fReMakeIC = g.fReMakeIC ;\r
-       fMakeSeveralIC = g.fMakeSeveralIC ;\r
-       fMakeInvMass = g.fMakeInvMass ;\r
-       fCalorimeter = g.fCalorimeter ;\r
-       fVertex[0] = g.fVertex[0];  \r
-       fVertex[1] = g.fVertex[1];  \r
-       fVertex[2] = g.fVertex[2]; \r
-        \r
-       fhConeSumPt = g.fhConeSumPt ;\r
-       fhPtInCone = g.fhPtInCone;\r
-       \r
-       fhPtIso = g.fhPtIso ; \r
-       fhPhiIso = g.fhPhiIso ;\r
-       fhEtaIso = g.fhEtaIso ;\r
-       \r
-       fhPtIsoPrompt = g.fhPtIsoPrompt;\r
-       fhPhiIsoPrompt = g.fhPhiIsoPrompt;\r
-       fhEtaIsoPrompt = g.fhEtaIsoPrompt; \r
-       fhPtIsoFragmentation = g.fhPtIsoFragmentation;\r
-       fhPhiIsoFragmentation = g.fhPhiIsoFragmentation;\r
-       fhEtaIsoFragmentation = g.fhEtaIsoFragmentation; \r
-       fhPtIsoPi0Decay = g.fhPtIsoPi0Decay;\r
-       fhPhiIsoPi0Decay = g.fhPhiIsoPi0Decay;\r
-       fhEtaIsoPi0Decay = g.fhEtaIsoPi0Decay; \r
-       fhPtIsoOtherDecay = g.fhPtIsoOtherDecay;\r
-       fhPhiIsoOtherDecay = g.fhPhiIsoOtherDecay;\r
-       fhEtaIsoOtherDecay = g.fhEtaIsoOtherDecay; \r
-       fhPtIsoConversion = g. fhPtIsoConversion;\r
-       fhPhiIsoConversion = g.fhPhiIsoConversion;\r
-       fhEtaIsoConversion = g.fhEtaIsoConversion; \r
-       fhPtIsoUnknown = g.fhPtIsoUnknown;\r
-       fhPhiIsoUnknown = g.fhPhiIsoUnknown;\r
-       fhEtaIsoUnknown = g.fhEtaIsoUnknown; \r
-       \r
-       //Several IC\r
-       fNCones = g.fNCones ;\r
-       fNPtThresFrac = g.fNPtThresFrac ; \r
-       \r
-       for(Int_t i = 0; i < fNCones ; i++){\r
-               fConeSizes[i] =  g.fConeSizes[i];\r
-               fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;\r
-               \r
-               fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i]; \r
-               fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i]; \r
-               fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i]; \r
-               fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i]; \r
-               fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i]; \r
-               fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i]; \r
-               \r
-               for(Int_t j = 0; j < fNPtThresFrac ; j++){\r
-                       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;\r
-                       fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;\r
-                       \r
-                       fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j]; \r
-                       fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j]; \r
-                       fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j]; \r
-                       fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j]; \r
-                       fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j]; \r
-                       fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j]; \r
-                       \r
-                       fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j]; \r
-                       fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j]; \r
-                       fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j]; \r
-                       fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j]; \r
-                       fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j]; \r
-                       fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j]; \r
-                       \r
-               }\r
-       }\r
-       \r
-       for(Int_t i = 0; i < fNPtThresFrac ; i++){\r
-               fPtThresholds[i]=   g.fPtThresholds[i];\r
-               fPtFractions[i]=   g.fPtFractions[i];\r
-       }\r
-       \r
-       \r
-       fHistoNPtSumBins    = g.fHistoNPtSumBins;\r
-       fHistoPtSumMax      = g.fHistoPtSumMax; \r
-       fHistoPtSumMin      = g.fHistoPtSumMax;\r
-       fHistoNPtInConeBins = g.fHistoNPtInConeBins;\r
-       fHistoPtInConeMax   = g.fHistoPtInConeMax;\r
-       fHistoPtInConeMin   = g.fHistoPtInConeMin;      \r
-       \r
-       return *this;\r
-       \r
-}\r
-\r
-//____________________________________________________________________________\r
-AliAnaParticleIsolation::~AliAnaParticleIsolation() \r
-{\r
-  //dtor\r
-  //do not delete histograms\r
-\r
- delete [] fConeSizes ; \r
- delete [] fPtThresholds ; \r
- delete [] fPtFractions ; \r
- delete [] fVertex ;\r
-}\r
-\r
-//_________________________________________________________________________\r
-Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4ParticleCorrelation * part1) const\r
-{\r
-       // Search if there is a companion decay photon to the candidate \r
-       // and discard it in such case\r
-       // Use it only if isolation candidates are photons\r
-       // Make sure that no selection on photon pt is done in the input aod photon list.\r
-       TLorentzVector mom1 = *(part1->Momentum());\r
-       TLorentzVector mom2 ;\r
-       for(Int_t jaod = 0; jaod < GetInputAODBranch()->GetEntriesFast(); jaod++){\r
-               if(iaod == jaod) continue ;\r
-               AliAODPWG4ParticleCorrelation * part2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jaod));\r
-               mom2 = *(part2->Momentum());\r
-               //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
-               if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){\r
-                       if(GetDebug() > 1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
-                       return kTRUE ;\r
-               }\r
-       }//loop\r
-       \r
-  return kFALSE;\r
-}\r
-\r
-//________________________________________________________________________\r
-TList *  AliAnaParticleIsolation::GetCreateOutputObjects()\r
-{  \r
-       // Create histograms to be saved in output file and \r
-       // store them in outputContainer\r
-       TList * outputContainer = new TList() ; \r
-       outputContainer->SetName("IsolatedParticleHistos") ; \r
-       \r
-       Int_t nptbins  = GetHistoNPtBins();\r
-       Int_t nphibins = GetHistoNPhiBins();\r
-       Int_t netabins = GetHistoNEtaBins();\r
-       Float_t ptmax  = GetHistoPtMax();\r
-       Float_t phimax = GetHistoPhiMax();\r
-       Float_t etamax = GetHistoEtaMax();\r
-       Float_t ptmin  = GetHistoPtMin();\r
-       Float_t phimin = GetHistoPhiMin();\r
-       Float_t etamin = GetHistoEtaMin();      \r
-       \r
-       Int_t nptsumbins    = fHistoNPtSumBins;\r
-       Float_t ptsummax    = fHistoPtSumMax;\r
-       Float_t ptsummin    = fHistoPtSumMin;   \r
-       Int_t nptinconebins = fHistoNPtInConeBins;\r
-       Float_t ptinconemax = fHistoPtInConeMax;\r
-       Float_t ptinconemin = fHistoPtInConeMin;\r
-       \r
-       if(!fMakeSeveralIC){\r
-               \r
-               fhConeSumPt  = new TH2F\r
-               ("hConePtSum","#Sigma p_{T}  in cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-               fhConeSumPt->SetYTitle("#Sigma p_{T}");\r
-               fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");\r
-               outputContainer->Add(fhConeSumPt) ;\r
-               \r
-               fhPtInCone  = new TH2F\r
-               ("hPtInCone","p_{T} in cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);\r
-               fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");\r
-               fhPtInCone->SetXTitle("p_{T #gamma} (GeV/c)");\r
-               outputContainer->Add(fhPtInCone) ;\r
-               \r
-               fhPtIso  = new TH1F("hPtIso","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-               fhPtIso->SetYTitle("N");\r
-               fhPtIso->SetXTitle("p_{T #gamma}(GeV/c)");\r
-               outputContainer->Add(fhPtIso) ; \r
-               \r
-               fhPhiIso  = new TH2F\r
-               ("hPhiIso","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-               fhPhiIso->SetYTitle("#phi");\r
-               fhPhiIso->SetXTitle("p_{T #gamma} (GeV/c)");\r
-               outputContainer->Add(fhPhiIso) ; \r
-               \r
-               fhEtaIso  = new TH2F\r
-               ("hEtaIso","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-               fhEtaIso->SetYTitle("#eta");\r
-               fhEtaIso->SetXTitle("p_{T #gamma} (GeV/c)");\r
-               outputContainer->Add(fhEtaIso) ;\r
-               \r
-               if(IsDataMC()){\r
-                       \r
-                       fhPtIsoPrompt  = new TH1F("hPtIsoPrompt","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-                       fhPtIsoPrompt->SetYTitle("N");\r
-                       fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");\r
-                       outputContainer->Add(fhPtIsoPrompt) ; \r
-                       \r
-                       fhPhiIsoPrompt  = new TH2F\r
-                       ("hPhiIsoPrompt","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-                       fhPhiIsoPrompt->SetYTitle("#phi");\r
-                       fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhPhiIsoPrompt) ; \r
-                       \r
-                       fhEtaIsoPrompt  = new TH2F\r
-                       ("hEtaIsoPrompt","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-                       fhEtaIsoPrompt->SetYTitle("#eta");\r
-                       fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhEtaIsoPrompt) ;\r
-                       \r
-                       fhPtIsoFragmentation  = new TH1F("hPtIsoFragmentation","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-                       fhPtIsoFragmentation->SetYTitle("N");\r
-                       fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");\r
-                       outputContainer->Add(fhPtIsoFragmentation) ; \r
-                       \r
-                       fhPhiIsoFragmentation  = new TH2F\r
-                       ("hPhiIsoFragmentation","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-                       fhPhiIsoFragmentation->SetYTitle("#phi");\r
-                       fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhPhiIsoFragmentation) ; \r
-                       \r
-                       fhEtaIsoFragmentation  = new TH2F\r
-                       ("hEtaIsoFragmentation","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-                       fhEtaIsoFragmentation->SetYTitle("#eta");\r
-                       fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhEtaIsoFragmentation) ;\r
-                       \r
-                       fhPtIsoPi0Decay  = new TH1F("hPtIsoPi0Decay","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-                       fhPtIsoPi0Decay->SetYTitle("N");\r
-                       fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");\r
-                       outputContainer->Add(fhPtIsoPi0Decay) ; \r
-                       \r
-                       fhPhiIsoPi0Decay  = new TH2F\r
-                       ("hPhiIsoPi0Decay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-                       fhPhiIsoPi0Decay->SetYTitle("#phi");\r
-                       fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhPhiIsoPi0Decay) ; \r
-                       \r
-                       fhEtaIsoPi0Decay  = new TH2F\r
-                       ("hEtaIsoPi0Decay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-                       fhEtaIsoPi0Decay->SetYTitle("#eta");\r
-                       fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhEtaIsoPi0Decay) ;\r
-                       \r
-                       fhPtIsoOtherDecay  = new TH1F("hPtIsoOtherDecay","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-                       fhPtIsoOtherDecay->SetYTitle("N");\r
-                       fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");\r
-                       outputContainer->Add(fhPtIsoOtherDecay) ; \r
-                       \r
-                       fhPhiIsoOtherDecay  = new TH2F\r
-                       ("hPhiIsoOtherDecay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-                       fhPhiIsoOtherDecay->SetYTitle("#phi");\r
-                       fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhPhiIsoOtherDecay) ; \r
-                       \r
-                       fhEtaIsoOtherDecay  = new TH2F\r
-                       ("hEtaIsoOtherDecay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-                       fhEtaIsoOtherDecay->SetYTitle("#eta");\r
-                       fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhEtaIsoOtherDecay) ;\r
-                       \r
-                       fhPtIsoConversion  = new TH1F("hPtIsoConversion","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-                       fhPtIsoConversion->SetYTitle("N");\r
-                       fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");\r
-                       outputContainer->Add(fhPtIsoConversion) ; \r
-                       \r
-                       fhPhiIsoConversion  = new TH2F\r
-                       ("hPhiIsoConversion","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-                       fhPhiIsoConversion->SetYTitle("#phi");\r
-                       fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhPhiIsoConversion) ; \r
-                       \r
-                       fhEtaIsoConversion  = new TH2F\r
-                       ("hEtaIsoConversion","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-                       fhEtaIsoConversion->SetYTitle("#eta");\r
-                       fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhEtaIsoConversion) ;\r
-                       \r
-                       fhPtIsoUnknown  = new TH1F("hPtIsoUnknown","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); \r
-                       fhPtIsoUnknown->SetYTitle("N");\r
-                       fhPtIsoUnknown->SetXTitle("p_{T #gamma}(GeV/c)");\r
-                       outputContainer->Add(fhPtIsoUnknown) ; \r
-                       \r
-                       fhPhiIsoUnknown  = new TH2F\r
-                       ("hPhiIsoUnknown","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
-                       fhPhiIsoUnknown->SetYTitle("#phi");\r
-                       fhPhiIsoUnknown->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhPhiIsoUnknown) ; \r
-                       \r
-                       fhEtaIsoUnknown  = new TH2F\r
-                       ("hEtaIsoUnknown","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
-                       fhEtaIsoUnknown->SetYTitle("#eta");\r
-                       fhEtaIsoUnknown->SetXTitle("p_{T #gamma} (GeV/c)");\r
-                       outputContainer->Add(fhEtaIsoUnknown) ;\r
-               }//Histos with MC\r
-               \r
-       }\r
-       \r
-       if(fMakeSeveralIC){\r
-               char name[128];\r
-               char title[128];\r
-               for(Int_t icone = 0; icone<fNCones; icone++){\r
-                       sprintf(name,"hPtSumIsolated_Cone_%d",icone);\r
-                       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                       fhPtSumIsolated[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                       outputContainer->Add(fhPtSumIsolated[icone]) ; \r
-                       \r
-                       if(IsDataMC()){\r
-                               sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);\r
-                               sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                               fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                               fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                               fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ; \r
-                               \r
-                               sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);\r
-                               sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                               fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                               fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                               fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ; \r
-                               \r
-                               sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);\r
-                               sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                               fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                               fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                               fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ; \r
-                               \r
-                               sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);\r
-                               sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                               fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                               fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                               fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ; \r
-                               \r
-                               sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);\r
-                               sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                               fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                               fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                               fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtSumIsolatedConversion[icone]) ; \r
-                               \r
-                               sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);\r
-                               sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);\r
-                               fhPtSumIsolatedUnknown[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
-                               fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");\r
-                               fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ; \r
-                               \r
-                       }//Histos with MC\r
-                       \r
-                       for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){ \r
-                               sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);\r
-                               sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                               fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                               fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtThresIsolated[icone][ipt]) ; \r
-                               \r
-                               sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);\r
-                               sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                               fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                               fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                               outputContainer->Add(fhPtFracIsolated[icone][ipt]) ; \r
-                               \r
-                               if(IsDataMC()){\r
-                                       sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;\r
-                                       \r
-                                       sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;\r
-                                       \r
-                                       sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtThresIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ; \r
-                                       \r
-                                       sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);\r
-                                       sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);\r
-                                       fhPtFracIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);\r
-                                       fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");\r
-                                       outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;  \r
-                                       \r
-                               }//Histos with MC\r
-                               \r
-                       }//icone loop\r
-               }//ipt loop\r
-       }\r
-       \r
-       //Keep neutral meson selection histograms if requiered\r
-       //Setting done in AliNeutralMesonSelection\r
-       if(fMakeInvMass && GetNeutralMesonSelection()){\r
-               TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;\r
-               if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())\r
-                       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;\r
-       }\r
-       \r
-       \r
-       //Save parameters used for analysis\r
-       TString parList ; //this will be list of parameters used for this analysis.\r
-       char onePar[255] ;\r
-       \r
-       sprintf(onePar,"--- AliAnaParticleIsolation ---\n") ;\r
-       parList+=onePar ;       \r
-       sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;\r
-       parList+=onePar ;\r
-       sprintf(onePar,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;\r
-       parList+=onePar ;\r
-       sprintf(onePar,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;\r
-       parList+=onePar ;\r
-       sprintf(onePar,"fMakeInvMass=%d (Flag for rejection of candidates with a pi0 inv mass pair) \n",fMakeInvMass) ;\r
-       parList+=onePar ;\r
-    \r
-       if(fMakeSeveralIC){\r
-               sprintf(onePar,"fNCones =%d (Number of cone sizes) \n",fNCones) ;\r
-               parList+=onePar ;\r
-               sprintf(onePar,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;\r
-               parList+=onePar ;\r
-               \r
-               for(Int_t icone = 0; icone < fNCones ; icone++){\r
-                       sprintf(onePar,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;\r
-                       parList+=onePar ;       \r
-               }\r
-               for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){\r
-                       sprintf(onePar,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;\r
-                       parList+=onePar ;       \r
-               }\r
-               for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){\r
-                       sprintf(onePar,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;\r
-                       parList+=onePar ;       \r
-               }               \r
-       }\r
-       \r
-       //Get parameters set in base class.\r
-       parList += GetBaseParametersList() ;\r
-       \r
-       //Get parameters set in IC class.\r
-       if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;\r
-       \r
-       TObjString *oString= new TObjString(parList) ;\r
-       outputContainer->Add(oString);\r
-       \r
-       \r
-       return outputContainer ;\r
-       \r
-}\r
-\r
-//__________________________________________________________________\r
-void  AliAnaParticleIsolation::MakeAnalysisFillAOD() \r
-{\r
-       //Do analysis and fill aods\r
-       //Search for the isolated photon in fCalorimeter with pt > GetMinPt()\r
-       \r
-       if(!GetInputAODBranch())\r
-               AliFatal(Form("ParticleIsolation::FillAOD: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data()));\r
-       \r
-       Int_t n = 0, nfrac = 0;\r
-       Bool_t isolated = kFALSE ; \r
-       Float_t coneptsum = 0 ;\r
-       TClonesArray * pl = new TClonesArray; \r
-       \r
-       //Select the calorimeter for candidate isolation with neutral particles\r
-       if(fCalorimeter == "PHOS")\r
-               pl = GetAODPHOS();\r
-       else if (fCalorimeter == "EMCAL")\r
-               pl = GetAODEMCAL();\r
-       \r
-       //Get vertex for photon momentum calculation\r
-       if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(fVertex);\r
-    \r
-       //Loop on AOD branch, filled previously in AliAnaPhoton\r
-       TLorentzVector mom ;\r
-       \r
-       for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast(); iaod++){\r
-               AliAODPWG4ParticleCorrelation * aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
-                               \r
-               //If too small or too large pt, skip\r
-               if(aod->Pt() < GetMinPt() || aod->Pt() > GetMaxPt() ) continue ; \r
-               \r
-               //Check invariant mass, if pi0, skip.\r
-               Bool_t decay = kFALSE ;\r
-               if(fMakeInvMass) decay = CheckInvMass(iaod,aod);\r
-               if(decay) continue ;\r
-               \r
-               n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;\r
-               GetIsolationCut()->MakeIsolationCut((TSeqCollection*)GetAODCTS(), (TSeqCollection*)pl, \r
-                                                                                               fVertex, kTRUE, aod, n,nfrac,coneptsum, isolated);\r
-               aod->SetIsolated(isolated);\r
-               \r
-       }//loop\r
-       \r
-       if(GetDebug() > 1) printf("End fill AODs ");  \r
-       \r
-}\r
-\r
-//__________________________________________________________________\r
-void  AliAnaParticleIsolation::MakeAnalysisFillHistograms() \r
-{\r
-       //Do analysis and fill histograms\r
-       Int_t n = 0, nfrac = 0;\r
-       Bool_t isolated = kFALSE ; \r
-       Float_t coneptsum = 0 ;\r
-       //cout<<"MakeAnalysisFillHistograms"<<endl;\r
-       \r
-       //Loop on stored AOD \r
-       Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
-       if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);\r
-       for(Int_t iaod = 0; iaod < naod ; iaod++){\r
-               AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
-               Bool_t isolation = aod->IsIsolated();              \r
-               Float_t ptcluster  = aod->Pt();\r
-               Float_t phicluster = aod->Phi();\r
-               Float_t etacluster = aod->Eta();\r
-               \r
-               //is pt too small skip it\r
-               if(ptcluster < GetMinPt() || ptcluster > GetMaxPt() ) continue ; \r
-                               \r
-               if(fMakeSeveralIC) {\r
-                       //Analysis of multiple IC at same time\r
-                       MakeSeveralICAnalysis(aod);\r
-                       continue;\r
-               }\r
-               else if(fReMakeIC){\r
-                       //In case a more strict IC is needed in the produced AOD\r
-                       n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;\r
-                       GetIsolationCut()->MakeIsolationCut((TSeqCollection*)aod->GetRefIsolationConeTracks(), \r
-                                                                                               (TSeqCollection*)aod->GetRefIsolationConeClusters(), \r
-                                                                                               fVertex, kFALSE, aod,\r
-                                                                                               n,nfrac,coneptsum, isolated);\r
-                       fhConeSumPt->Fill(ptcluster,coneptsum);       \r
-               }\r
-               \r
-               //Fill pt distribution of particles in cone\r
-               //Tracks\r
-               coneptsum=0;\r
-               for(Int_t itrack=0; itrack < (aod->GetRefIsolationConeTracks())->GetEntriesFast(); itrack++){\r
-                       AliAODTrack* track = (AliAODTrack *)((aod->GetRefIsolationConeTracks())->At(itrack));\r
-                       fhPtInCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));\r
-                       coneptsum+=track->Pt();\r
-               }\r
-               //CaloClusters\r
-               TLorentzVector mom ;\r
-               for(Int_t icalo=0; icalo < (aod->GetRefIsolationConeClusters())->GetEntriesFast(); icalo++){\r
-                       AliAODCaloCluster* calo = (AliAODCaloCluster *)((aod->GetRefIsolationConeClusters())->At(icalo));\r
-                       calo->GetMomentum(mom,fVertex);//Assume that come from vertex in straight line\r
-                       fhPtInCone->Fill(ptcluster, mom.Pt());\r
-                       coneptsum+=mom.Pt();\r
-               }\r
-               \r
-               if(!fReMakeIC) fhConeSumPt->Fill(ptcluster,coneptsum);\r
-               \r
-               if(isolation){    \r
-                       fhPtIso  ->Fill(ptcluster);\r
-                       fhPhiIso ->Fill(ptcluster,phicluster);\r
-                       fhEtaIso ->Fill(ptcluster,etacluster);\r
-                       \r
-                       if(IsDataMC()){\r
-                               Int_t tag =aod->GetTag();\r
-                               \r
-                               if(tag == AliMCAnalysisUtils::kMCPrompt){\r
-                                       fhPtIsoPrompt  ->Fill(ptcluster);\r
-                                       fhPhiIsoPrompt ->Fill(ptcluster,phicluster);\r
-                                       fhEtaIsoPrompt ->Fill(ptcluster,etacluster);\r
-                               }\r
-                               else if(tag == AliMCAnalysisUtils::kMCFragmentation)\r
-                               {\r
-                                       fhPtIsoFragmentation  ->Fill(ptcluster);\r
-                                       fhPhiIsoFragmentation ->Fill(ptcluster,phicluster);\r
-                                       fhEtaIsoFragmentation ->Fill(ptcluster,etacluster);\r
-                               }\r
-                               else if(tag == AliMCAnalysisUtils::kMCPi0Decay)\r
-                               {\r
-                                       fhPtIsoPi0Decay  ->Fill(ptcluster);\r
-                                       fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster);\r
-                                       fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster);\r
-                               }\r
-                               else if(tag == AliMCAnalysisUtils::kMCEtaDecay || tag == AliMCAnalysisUtils::kMCOtherDecay)\r
-                               {\r
-                                       fhPtIsoOtherDecay  ->Fill(ptcluster);\r
-                                       fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster);\r
-                                       fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster);\r
-                               }\r
-                               else if(tag == AliMCAnalysisUtils::kMCConversion)\r
-                               {\r
-                                       fhPtIsoConversion  ->Fill(ptcluster);\r
-                                       fhPhiIsoConversion ->Fill(ptcluster,phicluster);\r
-                                       fhEtaIsoConversion ->Fill(ptcluster,etacluster);\r
-                               }\r
-                               else\r
-                               {\r
-                                       fhPtIsoUnknown  ->Fill(ptcluster);\r
-                                       fhPhiIsoUnknown ->Fill(ptcluster,phicluster);\r
-                                       fhEtaIsoUnknown ->Fill(ptcluster,etacluster);\r
-                               }\r
-                       }//Histograms with MC\r
-                       \r
-               }//Isolated histograms\r
-               \r
-       }// aod loop\r
-       \r
-}\r
-\r
-//____________________________________________________________________________\r
-void AliAnaParticleIsolation::InitParameters()\r
-{\r
-  \r
-  //Initialize the parameters of the analysis.\r
-  SetInputAODName("photons");\r
-  \r
-  fCalorimeter = "PHOS" ;\r
-  fReMakeIC = kFALSE ;\r
-  fMakeSeveralIC = kFALSE ;\r
-  fMakeInvMass = kFALSE ;\r
-  \r
- //----------- Several IC-----------------\r
-  fNCones           = 5 ; \r
-  fNPtThresFrac     = 5 ; \r
-  fConeSizes[0]    = 0.1;  fConeSizes[1]     = 0.2;  fConeSizes[2]    = 0.3; fConeSizes[3]    = 0.4;  fConeSizes[4]    = 0.5;\r
-  fPtThresholds[0] = 1.;   fPtThresholds[1] = 2.;    fPtThresholds[2] = 3.;  fPtThresholds[3] = 4.;   fPtThresholds[4] = 5.; \r
-  fPtFractions[0]  = 0.05; fPtFractions[1]  = 0.075; fPtFractions[2]  = 0.1; fPtFractions[3]  = 1.25; fPtFractions[4]  = 1.5; \r
-\r
-//------------- Histograms settings -------\r
-  fHistoNPtSumBins = 100 ;\r
-  fHistoPtSumMax   = 50 ;\r
-  fHistoPtSumMin   = 0.  ;\r
-\r
-  fHistoNPtInConeBins = 100 ;\r
-  fHistoPtInConeMax   = 50 ;\r
-  fHistoPtInConeMin   = 0.  ;\r
-\r
-}\r
-\r
-//__________________________________________________________________\r
-void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph) \r
-{\r
-       //Isolation Cut Analysis for both methods and different pt cuts and cones\r
-       Float_t ptC = ph->Pt(); \r
-       Int_t tag   = ph->GetTag();\r
-       \r
-       //Keep original setting used when filling AODs, reset at end of analysis  \r
-       Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();\r
-       Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();\r
-       Float_t rorg       = GetIsolationCut()->GetConeSize();\r
-       \r
-       Float_t coneptsum = 0 ;  \r
-       Int_t n[10][10];//[fNCones][fNPtThresFrac];\r
-       Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];\r
-       Bool_t  isolated   = kFALSE;\r
-       \r
-       //Loop on cone sizes\r
-       for(Int_t icone = 0; icone<fNCones; icone++){\r
-               GetIsolationCut()->SetConeSize(fConeSizes[icone]);\r
-               coneptsum = 0 ;\r
-               //Loop on ptthresholds\r
-               for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){\r
-                       n[icone][ipt]=0;\r
-                       nfrac[icone][ipt]=0;\r
-                       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);\r
-                       GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(), \r
-                                                                                               (TSeqCollection*)ph->GetRefIsolationConeClusters(), \r
-                                                                                               fVertex, kFALSE, ph,\r
-                                                                                               n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);\r
-                       //Normal ptThreshold cut\r
-                       if(n[icone][ipt] == 0) {\r
-                               fhPtThresIsolated[icone][ipt]->Fill(ptC);\r
-                               if(IsDataMC()){\r
-                                       if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;\r
-                                       else  fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC) ;\r
-                               }\r
-                       }\r
-                       \r
-                       //Pt threshold on pt cand/ pt in cone fraction\r
-                       if(nfrac[icone][ipt] == 0) {\r
-                               fhPtFracIsolated[icone][ipt]->Fill(ptC);\r
-                               if(IsDataMC()){\r
-                                       if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
-                                       else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;\r
-                                       else  fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ;\r
-                               }\r
-                       }\r
-               }//pt thresh loop\r
-               \r
-               //Sum in cone histograms\r
-               fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;\r
-               if(IsDataMC()){\r
-                       if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;\r
-                       else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;\r
-                       else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;\r
-                       else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;\r
-                       else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;\r
-                       else  fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;\r
-               }\r
-               \r
-       }//cone size loop\r
-       \r
-       //Reset original parameters for AOD analysis\r
-       GetIsolationCut()->SetPtThreshold(ptthresorg);\r
-       GetIsolationCut()->SetPtFraction(ptfracorg);\r
-       GetIsolationCut()->SetConeSize(rorg);\r
-       \r
-}\r
-\r
-//__________________________________________________________________\r
-void AliAnaParticleIsolation::Print(const Option_t * opt) const\r
-{\r
-       \r
-       //Print some relevant parameters set for the analysis\r
-       if(! opt)\r
-               return;\r
-       \r
-       printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
-       \r
-       printf("Min Gamma pT       =     %2.2f\n",  GetMinPt()) ;\r
-       printf("Max Gamma pT       =     %3.2f\n",  GetMaxPt()) ;\r
-       \r
-       printf("ReMake Isolation          = %d \n",  fReMakeIC) ;\r
-       printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;\r
-       printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;\r
-       \r
-       if(fMakeSeveralIC){\r
-               printf("N Cone Sizes       =     %d\n", fNCones) ; \r
-               printf("Cone Sizes          =    \n") ;\r
-               for(Int_t i = 0; i < fNCones; i++)\r
-                       printf("  %1.2f;",  fConeSizes[i]) ;\r
-               printf("    \n") ;\r
-               \r
-               printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;\r
-               printf(" pT thresholds         =    \n") ;\r
-               for(Int_t i = 0; i < fNPtThresFrac; i++)\r
-                       printf("   %2.2f;",  fPtThresholds[i]) ;\r
-               \r
-               printf("    \n") ;\r
-               \r
-               printf(" pT fractions         =    \n") ;\r
-               for(Int_t i = 0; i < fNPtThresFrac; i++)\r
-                       printf("   %2.2f;",  fPtFractions[i]) ;\r
-               \r
-       }  \r
-       \r
-       printf("Histograms: %3.1f < pT sum < %3.1f,  Nbin = %d\n", fHistoPtSumMin,  fHistoPtSumMax,  fHistoNPtSumBins);\r
-       printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);\r
-       \r
-       printf("    \n") ;\r
-       \r
-} \r
-\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes hereby granted      *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id: AliAnaParticleIsolation.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
+
+//_________________________________________________________________________
+// Class for analysis of particle isolation
+// Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)
+//
+//  Class created from old AliPHOSGammaJet 
+//  (see AliRoot versions previous Release 4-09)
+//
+// -- Author: Gustavo Conesa (LNF-INFN) 
+//////////////////////////////////////////////////////////////////////////////
+  
+  
+// --- ROOT system --- 
+#include <TClonesArray.h>
+#include <TList.h>
+#include <TObjString.h>
+#include <TH2F.h>
+
+// --- Analysis system --- 
+#include "AliAnaParticleIsolation.h" 
+#include "AliCaloTrackReader.h"
+#include "AliIsolationCut.h"
+#include "AliNeutralMesonSelection.h"
+#include "AliAODPWG4ParticleCorrelation.h"
+#include "AliMCAnalysisUtils.h"
+#include "AliAODTrack.h"
+#include "AliAODCaloCluster.h"
+
+ClassImp(AliAnaParticleIsolation)
+  
+//____________________________________________________________________________
+  AliAnaParticleIsolation::AliAnaParticleIsolation() : 
+    AliAnaPartCorrBaseClass(), fCalorimeter(""), fVertex(),
+    fReMakeIC(0), fMakeSeveralIC(0), fMakeInvMass(0),
+    fhPtIso(0),fhPhiIso(0),fhEtaIso(0), fhConeSumPt(0), fhPtInCone(0),
+    //Several IC
+    fNCones(0),fNPtThresFrac(0), fConeSizes(),  fPtThresholds(),  fPtFractions(), 
+    //MC
+    fhPtIsoPrompt(0),fhPhiIsoPrompt(0),fhEtaIsoPrompt(0), 
+    fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),
+    fhPtIsoFragmentation(0),fhPhiIsoFragmentation(0),fhEtaIsoFragmentation(0), 
+    fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),
+    fhPtIsoPi0Decay(0),fhPhiIsoPi0Decay(0),fhEtaIsoPi0Decay(0),
+    fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),
+    fhPtIsoOtherDecay(0),fhPhiIsoOtherDecay(0),fhEtaIsoOtherDecay(0), 
+    fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),
+    fhPtIsoConversion(0),fhPhiIsoConversion(0),fhEtaIsoConversion(0), 
+    fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),
+    fhPtIsoUnknown(0),fhPhiIsoUnknown(0),fhEtaIsoUnknown(0), 
+    fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown(),
+    //Histograms settings
+    fHistoNPtSumBins(0),    fHistoPtSumMax(0.),    fHistoPtSumMin(0.),
+    fHistoNPtInConeBins(0), fHistoPtInConeMax(0.), fHistoPtInConeMin(0.)
+{
+  //default ctor
+  
+  //Initialize parameters
+  InitParameters();
+  
+  fVertex[0] = 0.;  fVertex[1] = 0.;  fVertex[2] = 0.;  
+  
+  for(Int_t i = 0; i < 5 ; i++){ 
+    fConeSizes[i] = 0 ; 
+    fhPtSumIsolated[i] = 0 ;  
+    
+    fhPtSumIsolatedPrompt[i] = 0 ;  
+    fhPtSumIsolatedFragmentation[i] = 0 ;  
+    fhPtSumIsolatedPi0Decay[i] = 0 ;  
+    fhPtSumIsolatedOtherDecay[i] = 0 ;  
+    fhPtSumIsolatedConversion[i] = 0 ;  
+    fhPtSumIsolatedUnknown[i] = 0 ;  
+    
+    for(Int_t j = 0; j < 5 ; j++){ 
+      fhPtThresIsolated[i][j] = 0 ;  
+      fhPtFracIsolated[i][j] = 0 ; 
+      
+      fhPtThresIsolatedPrompt[i][j] = 0 ;  
+      fhPtThresIsolatedFragmentation[i][j] = 0 ; 
+      fhPtThresIsolatedPi0Decay[i][j] = 0 ;  
+      fhPtThresIsolatedOtherDecay[i][j] = 0 ;  
+      fhPtThresIsolatedConversion[i][j] = 0 ;  
+      fhPtThresIsolatedUnknown[i][j] = 0 ;  
+  
+      fhPtFracIsolatedPrompt[i][j] = 0 ;  
+      fhPtFracIsolatedFragmentation[i][j] = 0 ;  
+      fhPtFracIsolatedPi0Decay[i][j] = 0 ;  
+      fhPtFracIsolatedOtherDecay[i][j] = 0 ;  
+      fhPtFracIsolatedConversion[i][j] = 0 ;
+      fhPtFracIsolatedUnknown[i][j] = 0 ;  
+    }  
+  } 
+  
+  for(Int_t i = 0; i < 5 ; i++){ 
+    fPtFractions[i]=  0 ; 
+    fPtThresholds[i]= 0 ; 
+  } 
+
+
+}
+
+//____________________________________________________________________________
+AliAnaParticleIsolation::AliAnaParticleIsolation(const AliAnaParticleIsolation & g) : 
+  AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),fVertex(),
+  fReMakeIC(g.fReMakeIC), fMakeSeveralIC(g.fMakeSeveralIC),  fMakeInvMass(g.fMakeInvMass),
+  fhPtIso(g.fhPtIso),fhPhiIso(g.fhPhiIso),fhEtaIso(g.fhEtaIso), 
+  fhConeSumPt(g.fhConeSumPt), fhPtInCone(g.fhPtInCone),
+  //Several IC
+  fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(),  fPtFractions(), 
+  fhPtThresIsolated(),  fhPtFracIsolated(), fhPtSumIsolated(),
+  //MC
+  fhPtIsoPrompt(g.fhPtIsoPrompt),fhPhiIsoPrompt(g.fhPhiIsoPrompt),fhEtaIsoPrompt(g.fhEtaIsoPrompt), 
+  fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),
+  fhPtIsoFragmentation(g.fhPtIsoFragmentation),fhPhiIsoFragmentation(g.fhPhiIsoFragmentation),fhEtaIsoFragmentation(g.fhEtaIsoFragmentation),
+  fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),
+  fhPtIsoPi0Decay(g.fhPtIsoPi0Decay),fhPhiIsoPi0Decay(g.fhPhiIsoPi0Decay),fhEtaIsoPi0Decay(g.fhEtaIsoPi0Decay), 
+  fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),
+  fhPtIsoOtherDecay(g.fhPtIsoOtherDecay),fhPhiIsoOtherDecay(g.fhPhiIsoOtherDecay),fhEtaIsoOtherDecay(g.fhEtaIsoOtherDecay), 
+  fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),
+  fhPtIsoConversion(g. fhPtIsoConversion),fhPhiIsoConversion(g.fhPhiIsoConversion),fhEtaIsoConversion(g.fhEtaIsoConversion), 
+  fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),
+  fhPtIsoUnknown(g.fhPtIsoUnknown),fhPhiIsoUnknown(g.fhPhiIsoUnknown),fhEtaIsoUnknown(g.fhEtaIsoUnknown), 
+  fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown(),
+  //Histograms
+  fHistoNPtSumBins(g.fHistoNPtSumBins), fHistoPtSumMax(g.fHistoPtSumMax), fHistoPtSumMin(g.fHistoPtSumMax),
+  fHistoNPtInConeBins(g.fHistoNPtInConeBins), fHistoPtInConeMax(g.fHistoPtInConeMax), fHistoPtInConeMin(g.fHistoPtInConeMin)
+{  
+  // cpy ctor
+  
+  fVertex[0] = g.fVertex[0];  
+  fVertex[1] = g.fVertex[1];  
+  fVertex[2] = g.fVertex[2];  
+  
+  //Several IC
+  for(Int_t i = 0; i < fNCones ; i++){
+    fConeSizes[i] =  g.fConeSizes[i];
+    fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; 
+    
+    fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i]; 
+    fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i]; 
+    fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i]; 
+    fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i]; 
+    fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i]; 
+    fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i]; 
+    
+    for(Int_t j = 0; j < fNPtThresFrac ; j++){
+      fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j]; 
+      fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];
+      
+      fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j]; 
+      fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j]; 
+      fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j]; 
+      fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j]; 
+      fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j]; 
+      fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j]; 
+      
+      fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j]; 
+      fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j]; 
+      fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j]; 
+      fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j]; 
+      fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j]; 
+      fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j]; 
+                       
+    } 
+  }
+  
+  for(Int_t i = 0; i < fNPtThresFrac ; i++){
+    fPtFractions[i]=   g.fPtFractions[i];
+    fPtThresholds[i]=   g.fPtThresholds[i];
+  }
+  
+}
+
+//_________________________________________________________________________
+AliAnaParticleIsolation & AliAnaParticleIsolation::operator = (const AliAnaParticleIsolation & g)
+{
+  // assignment operator
+  
+  if(&g == this) return *this;
+  
+  fReMakeIC = g.fReMakeIC ;
+  fMakeSeveralIC = g.fMakeSeveralIC ;
+  fMakeInvMass = g.fMakeInvMass ;
+  fCalorimeter = g.fCalorimeter ;
+  fVertex[0] = g.fVertex[0];  
+  fVertex[1] = g.fVertex[1];  
+  fVertex[2] = g.fVertex[2]; 
+  
+  fhConeSumPt = g.fhConeSumPt ;
+  fhPtInCone = g.fhPtInCone;
+  
+  fhPtIso = g.fhPtIso ; 
+  fhPhiIso = g.fhPhiIso ;
+  fhEtaIso = g.fhEtaIso ;
+  
+  fhPtIsoPrompt = g.fhPtIsoPrompt;
+  fhPhiIsoPrompt = g.fhPhiIsoPrompt;
+  fhEtaIsoPrompt = g.fhEtaIsoPrompt; 
+  fhPtIsoFragmentation = g.fhPtIsoFragmentation;
+  fhPhiIsoFragmentation = g.fhPhiIsoFragmentation;
+  fhEtaIsoFragmentation = g.fhEtaIsoFragmentation; 
+  fhPtIsoPi0Decay = g.fhPtIsoPi0Decay;
+  fhPhiIsoPi0Decay = g.fhPhiIsoPi0Decay;
+  fhEtaIsoPi0Decay = g.fhEtaIsoPi0Decay; 
+  fhPtIsoOtherDecay = g.fhPtIsoOtherDecay;
+  fhPhiIsoOtherDecay = g.fhPhiIsoOtherDecay;
+  fhEtaIsoOtherDecay = g.fhEtaIsoOtherDecay; 
+  fhPtIsoConversion = g. fhPtIsoConversion;
+  fhPhiIsoConversion = g.fhPhiIsoConversion;
+  fhEtaIsoConversion = g.fhEtaIsoConversion; 
+  fhPtIsoUnknown = g.fhPtIsoUnknown;
+  fhPhiIsoUnknown = g.fhPhiIsoUnknown;
+  fhEtaIsoUnknown = g.fhEtaIsoUnknown; 
+  
+  //Several IC
+  fNCones = g.fNCones ;
+  fNPtThresFrac = g.fNPtThresFrac ; 
+  
+  for(Int_t i = 0; i < fNCones ; i++){
+    fConeSizes[i] =  g.fConeSizes[i];
+    fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
+    
+    fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i]; 
+    fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i]; 
+    fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i]; 
+    fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i]; 
+    fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i]; 
+    fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i]; 
+    
+    for(Int_t j = 0; j < fNPtThresFrac ; j++){
+      fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;
+      fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;
+      
+      fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j]; 
+      fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j]; 
+      fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j]; 
+      fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j]; 
+      fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j]; 
+      fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j]; 
+      
+      fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j]; 
+      fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j]; 
+      fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j]; 
+      fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j]; 
+      fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j]; 
+      fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j]; 
+      
+    }
+  }
+  
+  for(Int_t i = 0; i < fNPtThresFrac ; i++){
+    fPtThresholds[i]=   g.fPtThresholds[i];
+    fPtFractions[i]=   g.fPtFractions[i];
+  }
+  
+  
+  fHistoNPtSumBins    = g.fHistoNPtSumBins;
+  fHistoPtSumMax      = g.fHistoPtSumMax; 
+  fHistoPtSumMin      = g.fHistoPtSumMax;
+  fHistoNPtInConeBins = g.fHistoNPtInConeBins;
+  fHistoPtInConeMax   = g.fHistoPtInConeMax;
+  fHistoPtInConeMin   = g.fHistoPtInConeMin;   
+  
+  return *this;
+  
+}
+
+//____________________________________________________________________________
+AliAnaParticleIsolation::~AliAnaParticleIsolation() 
+{
+  //dtor
+  //do not delete histograms
+  
+  delete [] fConeSizes ; 
+  delete [] fPtThresholds ; 
+  delete [] fPtFractions ; 
+  delete [] fVertex ;
+}
+
+//_________________________________________________________________________
+Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4ParticleCorrelation * part1) const
+{
+  // Search if there is a companion decay photon to the candidate 
+  // and discard it in such case
+  // Use it only if isolation candidates are photons
+  // Make sure that no selection on photon pt is done in the input aod photon list.
+  TLorentzVector mom1 = *(part1->Momentum());
+  TLorentzVector mom2 ;
+  for(Int_t jaod = 0; jaod < GetInputAODBranch()->GetEntriesFast(); jaod++){
+    if(iaod == jaod) continue ;
+    AliAODPWG4ParticleCorrelation * part2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jaod));
+    mom2 = *(part2->Momentum());
+    //Select good pair (good phi, pt cuts, aperture and invariant mass)
+    if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
+      if(GetDebug() > 1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
+      return kTRUE ;
+    }
+  }//loop
+  
+  return kFALSE;
+}
+
+//________________________________________________________________________
+TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
+{  
+  // Create histograms to be saved in output file and 
+  // store them in outputContainer
+  TList * outputContainer = new TList() ; 
+  outputContainer->SetName("IsolatedParticleHistos") ; 
+  
+  Int_t nptbins  = GetHistoNPtBins();
+  Int_t nphibins = GetHistoNPhiBins();
+  Int_t netabins = GetHistoNEtaBins();
+  Float_t ptmax  = GetHistoPtMax();
+  Float_t phimax = GetHistoPhiMax();
+  Float_t etamax = GetHistoEtaMax();
+  Float_t ptmin  = GetHistoPtMin();
+  Float_t phimin = GetHistoPhiMin();
+  Float_t etamin = GetHistoEtaMin();   
+  
+  Int_t nptsumbins    = fHistoNPtSumBins;
+  Float_t ptsummax    = fHistoPtSumMax;
+  Float_t ptsummin    = fHistoPtSumMin;        
+  Int_t nptinconebins = fHistoNPtInConeBins;
+  Float_t ptinconemax = fHistoPtInConeMax;
+  Float_t ptinconemin = fHistoPtInConeMin;
+  
+  if(!fMakeSeveralIC){
+    
+    fhConeSumPt  = new TH2F
+      ("hConePtSum","#Sigma p_{T}  in cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+    fhConeSumPt->SetYTitle("#Sigma p_{T}");
+    fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhConeSumPt) ;
+    
+    fhPtInCone  = new TH2F
+      ("hPtInCone","p_{T} in cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
+    fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
+    fhPtInCone->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPtInCone) ;
+    
+    fhPtIso  = new TH1F("hPtIso","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+    fhPtIso->SetYTitle("N");
+    fhPtIso->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtIso) ; 
+    
+    fhPhiIso  = new TH2F
+      ("hPhiIso","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    fhPhiIso->SetYTitle("#phi");
+    fhPhiIso->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiIso) ; 
+    
+    fhEtaIso  = new TH2F
+      ("hEtaIso","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+    fhEtaIso->SetYTitle("#eta");
+    fhEtaIso->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaIso) ;
+    
+    if(IsDataMC()){
+      
+      fhPtIsoPrompt  = new TH1F("hPtIsoPrompt","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+      fhPtIsoPrompt->SetYTitle("N");
+      fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtIsoPrompt) ; 
+      
+      fhPhiIsoPrompt  = new TH2F
+       ("hPhiIsoPrompt","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+      fhPhiIsoPrompt->SetYTitle("#phi");
+      fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhPhiIsoPrompt) ; 
+      
+      fhEtaIsoPrompt  = new TH2F
+       ("hEtaIsoPrompt","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+      fhEtaIsoPrompt->SetYTitle("#eta");
+      fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhEtaIsoPrompt) ;
+      
+      fhPtIsoFragmentation  = new TH1F("hPtIsoFragmentation","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+      fhPtIsoFragmentation->SetYTitle("N");
+      fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtIsoFragmentation) ; 
+      
+      fhPhiIsoFragmentation  = new TH2F
+       ("hPhiIsoFragmentation","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+      fhPhiIsoFragmentation->SetYTitle("#phi");
+      fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhPhiIsoFragmentation) ; 
+      
+      fhEtaIsoFragmentation  = new TH2F
+       ("hEtaIsoFragmentation","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+      fhEtaIsoFragmentation->SetYTitle("#eta");
+      fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhEtaIsoFragmentation) ;
+      
+      fhPtIsoPi0Decay  = new TH1F("hPtIsoPi0Decay","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+      fhPtIsoPi0Decay->SetYTitle("N");
+      fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtIsoPi0Decay) ; 
+      
+      fhPhiIsoPi0Decay  = new TH2F
+       ("hPhiIsoPi0Decay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+      fhPhiIsoPi0Decay->SetYTitle("#phi");
+      fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhPhiIsoPi0Decay) ; 
+      
+      fhEtaIsoPi0Decay  = new TH2F
+       ("hEtaIsoPi0Decay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+      fhEtaIsoPi0Decay->SetYTitle("#eta");
+      fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhEtaIsoPi0Decay) ;
+      
+      fhPtIsoOtherDecay  = new TH1F("hPtIsoOtherDecay","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+      fhPtIsoOtherDecay->SetYTitle("N");
+      fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtIsoOtherDecay) ; 
+      
+      fhPhiIsoOtherDecay  = new TH2F
+       ("hPhiIsoOtherDecay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+      fhPhiIsoOtherDecay->SetYTitle("#phi");
+      fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhPhiIsoOtherDecay) ; 
+      
+      fhEtaIsoOtherDecay  = new TH2F
+       ("hEtaIsoOtherDecay","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+      fhEtaIsoOtherDecay->SetYTitle("#eta");
+      fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhEtaIsoOtherDecay) ;
+      
+      fhPtIsoConversion  = new TH1F("hPtIsoConversion","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+      fhPtIsoConversion->SetYTitle("N");
+      fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtIsoConversion) ; 
+      
+      fhPhiIsoConversion  = new TH2F
+       ("hPhiIsoConversion","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+      fhPhiIsoConversion->SetYTitle("#phi");
+      fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhPhiIsoConversion) ; 
+      
+      fhEtaIsoConversion  = new TH2F
+       ("hEtaIsoConversion","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+      fhEtaIsoConversion->SetYTitle("#eta");
+      fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhEtaIsoConversion) ;
+      
+      fhPtIsoUnknown  = new TH1F("hPtIsoUnknown","Isolated Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+      fhPtIsoUnknown->SetYTitle("N");
+      fhPtIsoUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtIsoUnknown) ; 
+      
+      fhPhiIsoUnknown  = new TH2F
+       ("hPhiIsoUnknown","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+      fhPhiIsoUnknown->SetYTitle("#phi");
+      fhPhiIsoUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhPhiIsoUnknown) ; 
+      
+      fhEtaIsoUnknown  = new TH2F
+       ("hEtaIsoUnknown","Isolated #phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+      fhEtaIsoUnknown->SetYTitle("#eta");
+      fhEtaIsoUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
+      outputContainer->Add(fhEtaIsoUnknown) ;
+    }//Histos with MC
+    
+  }
+  
+  if(fMakeSeveralIC){
+               char name[128];
+               char title[128];
+               for(Int_t icone = 0; icone<fNCones; icone++){
+                 sprintf(name,"hPtSumIsolated_Cone_%d",icone);
+                 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                 fhPtSumIsolated[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
+                 outputContainer->Add(fhPtSumIsolated[icone]) ; 
+                 
+                 if(IsDataMC()){
+                   sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);
+                   sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                   fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                   fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ; 
+                   
+                   sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);
+                   sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                   fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                   fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ; 
+                   
+                   sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);
+                   sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                   fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                   fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ; 
+                   
+                   sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);
+                   sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                   fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                   fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ; 
+                   
+                   sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);
+                   sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                   fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                   fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtSumIsolatedConversion[icone]) ; 
+                   
+                   sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);
+                   sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   fhPtSumIsolatedUnknown[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+                   fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+                   fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ; 
+                   
+                 }//Histos with MC
+                 
+                 for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){ 
+                   sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);
+                   sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                   fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                   fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
+                   
+                   sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);
+                   sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                   fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                   fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                   outputContainer->Add(fhPtFracIsolated[icone][ipt]) ; 
+                   
+                   if(IsDataMC()){
+                     sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
+                     
+                     sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
+                     
+                     sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtThresIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ; 
+                     
+                     sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);
+                     sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                     fhPtFracIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
+                     fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+                     outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;  
+                     
+                   }//Histos with MC
+                   
+                 }//icone loop
+               }//ipt loop
+  }
+  
+  //Keep neutral meson selection histograms if requiered
+  //Setting done in AliNeutralMesonSelection
+  if(fMakeInvMass && GetNeutralMesonSelection()){
+    TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
+    if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
+      for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
+  }
+  
+  
+  //Save parameters used for analysis
+  TString parList ; //this will be list of parameters used for this analysis.
+  char onePar[255] ;
+  
+  sprintf(onePar,"--- AliAnaParticleIsolation ---\n") ;
+  parList+=onePar ;    
+  sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
+  parList+=onePar ;
+  sprintf(onePar,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
+  parList+=onePar ;
+  sprintf(onePar,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
+  parList+=onePar ;
+  sprintf(onePar,"fMakeInvMass=%d (Flag for rejection of candidates with a pi0 inv mass pair) \n",fMakeInvMass) ;
+  parList+=onePar ;
+  
+  if(fMakeSeveralIC){
+    sprintf(onePar,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
+    parList+=onePar ;
+    sprintf(onePar,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
+    parList+=onePar ;
+    
+    for(Int_t icone = 0; icone < fNCones ; icone++){
+      sprintf(onePar,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
+      parList+=onePar ;        
+    }
+    for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
+      sprintf(onePar,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
+      parList+=onePar ;        
+    }
+    for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
+      sprintf(onePar,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
+      parList+=onePar ;        
+    }          
+  }
+  
+  //Get parameters set in base class.
+  parList += GetBaseParametersList() ;
+  
+  //Get parameters set in IC class.
+  if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
+  
+  TObjString *oString= new TObjString(parList) ;
+  outputContainer->Add(oString);
+  
+  
+  return outputContainer ;
+  
+}
+
+//__________________________________________________________________
+void  AliAnaParticleIsolation::MakeAnalysisFillAOD() 
+{
+  //Do analysis and fill aods
+  //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
+  
+  if(!GetInputAODBranch()){
+    printf("ParticleIsolation::FillAOD: No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
+    abort();
+  }
+  Int_t n = 0, nfrac = 0;
+  Bool_t isolated = kFALSE ; 
+  Float_t coneptsum = 0 ;
+  TRefArray * pl = new TRefArray(); 
+  
+  //Select the calorimeter for candidate isolation with neutral particles
+  if(fCalorimeter == "PHOS")
+    pl = GetAODPHOS();
+  else if (fCalorimeter == "EMCAL")
+    pl = GetAODEMCAL();
+  
+  //Get vertex for photon momentum calculation
+  if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(fVertex);
+  
+  //Loop on AOD branch, filled previously in AliAnaPhoton
+  TLorentzVector mom ;
+  
+  for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast(); iaod++){
+    AliAODPWG4ParticleCorrelation * aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
+    
+    //If too small or too large pt, skip
+    if(aod->Pt() < GetMinPt() || aod->Pt() > GetMaxPt() ) continue ; 
+    
+    //Check invariant mass, if pi0, skip.
+    Bool_t decay = kFALSE ;
+    if(fMakeInvMass) decay = CheckInvMass(iaod,aod);
+    if(decay) continue ;
+    
+    n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
+    GetIsolationCut()->MakeIsolationCut(GetAODCTS(),pl,
+                                       fVertex, kTRUE, aod, n,nfrac,coneptsum, isolated);
+    aod->SetIsolated(isolated);
+    
+  }//loop
+  
+  if(GetDebug() > 1) printf("End fill AODs ");  
+  
+}
+
+//__________________________________________________________________
+void  AliAnaParticleIsolation::MakeAnalysisFillHistograms() 
+{
+  //Do analysis and fill histograms
+  Int_t n = 0, nfrac = 0;
+  Bool_t isolated = kFALSE ; 
+  Float_t coneptsum = 0 ;
+  //cout<<"MakeAnalysisFillHistograms"<<endl;
+  
+  //Loop on stored AOD 
+  Int_t naod = GetInputAODBranch()->GetEntriesFast();
+  if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
+  for(Int_t iaod = 0; iaod < naod ; iaod++){
+    AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
+    Bool_t isolation = aod->IsIsolated();              
+    Float_t ptcluster  = aod->Pt();
+    Float_t phicluster = aod->Phi();
+    Float_t etacluster = aod->Eta();
+    
+    //is pt too small skip it
+    if(ptcluster < GetMinPt() || ptcluster > GetMaxPt() ) continue ; 
+    
+    if(fMakeSeveralIC) {
+      //Analysis of multiple IC at same time
+      MakeSeveralICAnalysis(aod);
+      continue;
+    }
+    else if(fReMakeIC){
+      //In case a more strict IC is needed in the produced AOD
+      n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
+      GetIsolationCut()->MakeIsolationCut(aod->GetRefIsolationConeTracks(), 
+                                         aod->GetRefIsolationConeClusters(), 
+                                         fVertex, kFALSE, aod,
+                                         n,nfrac,coneptsum, isolated);
+      fhConeSumPt->Fill(ptcluster,coneptsum);       
+    }
+    
+    //Fill pt distribution of particles in cone
+    //Tracks
+    coneptsum=0;
+    for(Int_t itrack=0; itrack < (aod->GetRefIsolationConeTracks())->GetEntriesFast(); itrack++){
+      AliAODTrack* track = (AliAODTrack *)((aod->GetRefIsolationConeTracks())->At(itrack));
+      fhPtInCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
+      coneptsum+=track->Pt();
+    }
+    //CaloClusters
+    TLorentzVector mom ;
+    for(Int_t icalo=0; icalo < (aod->GetRefIsolationConeClusters())->GetEntriesFast(); icalo++){
+      AliAODCaloCluster* calo = (AliAODCaloCluster *)((aod->GetRefIsolationConeClusters())->At(icalo));
+      calo->GetMomentum(mom,fVertex);//Assume that come from vertex in straight line
+      fhPtInCone->Fill(ptcluster, mom.Pt());
+      coneptsum+=mom.Pt();
+    }
+    
+    if(!fReMakeIC) fhConeSumPt->Fill(ptcluster,coneptsum);
+    
+    if(isolation){    
+      fhPtIso  ->Fill(ptcluster);
+      fhPhiIso ->Fill(ptcluster,phicluster);
+      fhEtaIso ->Fill(ptcluster,etacluster);
+      
+      if(IsDataMC()){
+       Int_t tag =aod->GetTag();
+       
+       if(tag == AliMCAnalysisUtils::kMCPrompt){
+         fhPtIsoPrompt  ->Fill(ptcluster);
+         fhPhiIsoPrompt ->Fill(ptcluster,phicluster);
+         fhEtaIsoPrompt ->Fill(ptcluster,etacluster);
+       }
+       else if(tag == AliMCAnalysisUtils::kMCFragmentation)
+         {
+           fhPtIsoFragmentation  ->Fill(ptcluster);
+           fhPhiIsoFragmentation ->Fill(ptcluster,phicluster);
+           fhEtaIsoFragmentation ->Fill(ptcluster,etacluster);
+         }
+       else if(tag == AliMCAnalysisUtils::kMCPi0Decay)
+         {
+           fhPtIsoPi0Decay  ->Fill(ptcluster);
+           fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster);
+           fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster);
+         }
+       else if(tag == AliMCAnalysisUtils::kMCEtaDecay || tag == AliMCAnalysisUtils::kMCOtherDecay)
+         {
+           fhPtIsoOtherDecay  ->Fill(ptcluster);
+           fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster);
+           fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster);
+         }
+       else if(tag == AliMCAnalysisUtils::kMCConversion)
+         {
+           fhPtIsoConversion  ->Fill(ptcluster);
+           fhPhiIsoConversion ->Fill(ptcluster,phicluster);
+           fhEtaIsoConversion ->Fill(ptcluster,etacluster);
+         }
+       else
+         {
+           fhPtIsoUnknown  ->Fill(ptcluster);
+           fhPhiIsoUnknown ->Fill(ptcluster,phicluster);
+           fhEtaIsoUnknown ->Fill(ptcluster,etacluster);
+         }
+      }//Histograms with MC
+      
+    }//Isolated histograms
+    
+  }// aod loop
+  
+}
+
+//____________________________________________________________________________
+void AliAnaParticleIsolation::InitParameters()
+{
+  
+  //Initialize the parameters of the analysis.
+  SetInputAODName("photons");
+  
+  fCalorimeter = "PHOS" ;
+  fReMakeIC = kFALSE ;
+  fMakeSeveralIC = kFALSE ;
+  fMakeInvMass = kFALSE ;
+  
+ //----------- Several IC-----------------
+  fNCones           = 5 ; 
+  fNPtThresFrac     = 5 ; 
+  fConeSizes[0]    = 0.1;  fConeSizes[1]     = 0.2;  fConeSizes[2]    = 0.3; fConeSizes[3]    = 0.4;  fConeSizes[4]    = 0.5;
+  fPtThresholds[0] = 1.;   fPtThresholds[1] = 2.;    fPtThresholds[2] = 3.;  fPtThresholds[3] = 4.;   fPtThresholds[4] = 5.; 
+  fPtFractions[0]  = 0.05; fPtFractions[1]  = 0.075; fPtFractions[2]  = 0.1; fPtFractions[3]  = 1.25; fPtFractions[4]  = 1.5; 
+
+//------------- Histograms settings -------
+  fHistoNPtSumBins = 100 ;
+  fHistoPtSumMax   = 50 ;
+  fHistoPtSumMin   = 0.  ;
+
+  fHistoNPtInConeBins = 100 ;
+  fHistoPtInConeMax   = 50 ;
+  fHistoPtInConeMin   = 0.  ;
+
+}
+
+//__________________________________________________________________
+void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph) 
+{
+  //Isolation Cut Analysis for both methods and different pt cuts and cones
+  Float_t ptC = ph->Pt();      
+  Int_t tag   = ph->GetTag();
+  
+  //Keep original setting used when filling AODs, reset at end of analysis  
+  Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
+  Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();
+  Float_t rorg       = GetIsolationCut()->GetConeSize();
+  
+  Float_t coneptsum = 0 ;  
+  Int_t n[10][10];//[fNCones][fNPtThresFrac];
+  Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
+  Bool_t  isolated   = kFALSE;
+  
+  //Loop on cone sizes
+  for(Int_t icone = 0; icone<fNCones; icone++){
+    GetIsolationCut()->SetConeSize(fConeSizes[icone]);
+    coneptsum = 0 ;
+    //Loop on ptthresholds
+    for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
+      n[icone][ipt]=0;
+      nfrac[icone][ipt]=0;
+      GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
+      GetIsolationCut()->MakeIsolationCut(ph->GetRefIsolationConeTracks(), 
+                                         ph->GetRefIsolationConeClusters(), 
+                                         fVertex, kFALSE, ph,
+                                         n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
+      //Normal ptThreshold cut
+      if(n[icone][ipt] == 0) {
+       fhPtThresIsolated[icone][ipt]->Fill(ptC);
+       if(IsDataMC()){
+         if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;
+         else  fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC) ;
+       }
+      }
+      
+      //Pt threshold on pt cand/ pt in cone fraction
+      if(nfrac[icone][ipt] == 0) {
+       fhPtFracIsolated[icone][ipt]->Fill(ptC);
+       if(IsDataMC()){
+         if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;
+         else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;
+         else  fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ;
+       }
+      }
+    }//pt thresh loop
+    
+    //Sum in cone histograms
+    fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
+    if(IsDataMC()){
+      if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
+      else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
+      else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
+      else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
+      else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
+      else  fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
+    }
+    
+  }//cone size loop
+  
+  //Reset original parameters for AOD analysis
+  GetIsolationCut()->SetPtThreshold(ptthresorg);
+  GetIsolationCut()->SetPtFraction(ptfracorg);
+  GetIsolationCut()->SetConeSize(rorg);
+  
+}
+
+//__________________________________________________________________
+void AliAnaParticleIsolation::Print(const Option_t * opt) const
+{
+  
+  //Print some relevant parameters set for the analysis
+  if(! opt)
+    return;
+  
+  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
+  
+  printf("Min Gamma pT       =     %2.2f\n",  GetMinPt()) ;
+  printf("Max Gamma pT       =     %3.2f\n",  GetMaxPt()) ;
+  
+  printf("ReMake Isolation          = %d \n",  fReMakeIC) ;
+  printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;
+  printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;
+  
+  if(fMakeSeveralIC){
+    printf("N Cone Sizes       =     %d\n", fNCones) ; 
+    printf("Cone Sizes          =    \n") ;
+    for(Int_t i = 0; i < fNCones; i++)
+      printf("  %1.2f;",  fConeSizes[i]) ;
+    printf("    \n") ;
+    
+    printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
+    printf(" pT thresholds         =    \n") ;
+    for(Int_t i = 0; i < fNPtThresFrac; i++)
+      printf("   %2.2f;",  fPtThresholds[i]) ;
+    
+    printf("    \n") ;
+    
+    printf(" pT fractions         =    \n") ;
+    for(Int_t i = 0; i < fNPtThresFrac; i++)
+      printf("   %2.2f;",  fPtFractions[i]) ;
+    
+  }  
+  
+  printf("Histograms: %3.1f < pT sum < %3.1f,  Nbin = %d\n", fHistoPtSumMin,  fHistoPtSumMax,  fHistoNPtSumBins);
+  printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);
+  
+  printf("    \n") ;
+  
+} 
+