]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrDep/AliAnaParticleIsolation.cxx
move from AliAODTrack to AliVTrack
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleIsolation.cxx
index 315e5840d0356c6f295ef757c9e47406aa73b134..0a8bc03b021afbdc42c66f4f27d65d5fe85c137b 100755 (executable)
@@ -22,6 +22,8 @@
 //  (see AliRoot versions previous Release 4-09)
 //
 // -- Author: Gustavo Conesa (LNF-INFN) 
+
+//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
 //////////////////////////////////////////////////////////////////////////////
   
   
 #include "AliNeutralMesonSelection.h"
 #include "AliAODPWG4ParticleCorrelation.h"
 #include "AliMCAnalysisUtils.h"
-#include "AliAODTrack.h"
-#include "AliAODCaloCluster.h"
+#include "AliVTrack.h"
+#include "AliVCluster.h"
 
 ClassImp(AliAnaParticleIsolation)
   
 //____________________________________________________________________________
   AliAnaParticleIsolation::AliAnaParticleIsolation() : 
-    AliAnaPartCorrBaseClass(), fCalorimeter(""), fVertex(),
+    AliAnaPartCorrBaseClass(), fCalorimeter(""), 
     fReMakeIC(0), fMakeSeveralIC(0), fMakeInvMass(0),
     fhPtIso(0),fhPhiIso(0),fhEtaIso(0), fhConeSumPt(0), fhPtInCone(0),
+    fhFRConeSumPt(0), fhPtInFRCone(0),
     //Several IC
     fNCones(0),fNPtThresFrac(0), fConeSizes(),  fPtThresholds(),  fPtFractions(), 
     //MC
@@ -73,9 +76,7 @@ ClassImp(AliAnaParticleIsolation)
   
   //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 ;  
@@ -116,186 +117,20 @@ ClassImp(AliAnaParticleIsolation)
 
 }
 
-//____________________________________________________________________________
-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 ;
+  //delete [] fConeSizes ; 
+  //delete [] fPtThresholds ; 
+  //delete [] fPtFractions ; 
+
 }
 
 //_________________________________________________________________________
-Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4Particle * part1) const
+Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4Particle * part1)
 {
   // Search if there is a companion decay photon to the candidate 
   // and discard it in such case
@@ -317,6 +152,55 @@ Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4P
   return kFALSE;
 }
 
+//________________________________________________________________________
+TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
+{ 
+       //Save parameters used for analysis
+        TString parList ; //this will be list of parameters used for this analysis.
+   const Int_t buffersize = 255;
+        char onePar[buffersize] ;
+        
+        snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
+        parList+=onePar ;      
+        snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
+        parList+=onePar ;
+        snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
+        parList+=onePar ;
+        snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
+        parList+=onePar ;
+        snprintf(onePar, buffersize,"fMakeInvMass=%d (Flag for rejection of candidates with a pi0 inv mass pair) \n",fMakeInvMass) ;
+        parList+=onePar ;
+        
+        if(fMakeSeveralIC){
+          snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
+          parList+=onePar ;
+          snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
+          parList+=onePar ;
+          
+          for(Int_t icone = 0; icone < fNCones ; icone++){
+            snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
+            parList+=onePar ;  
+          }
+          for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
+            snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
+            parList+=onePar ;  
+          }
+          for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
+            snprintf(onePar, buffersize,"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() ;
+        
+        return new TObjString(parList) ;
+       
+}
+
 //________________________________________________________________________
 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
 {  
@@ -325,9 +209,9 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
   TList * outputContainer = new TList() ; 
   outputContainer->SetName("IsolatedParticleHistos") ; 
   
-  Int_t nptbins  = GetHistoNPtBins();
-  Int_t nphibins = GetHistoNPhiBins();
-  Int_t netabins = GetHistoNEtaBins();
+  Int_t nptbins  = GetHistoPtBins();
+  Int_t nphibins = GetHistoPhiBins();
+  Int_t netabins = GetHistoEtaBins();
   Float_t ptmax  = GetHistoPtMax();
   Float_t phimax = GetHistoPhiMax();
   Float_t etamax = GetHistoEtaMax();
@@ -356,6 +240,18 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
     fhPtInCone->SetXTitle("p_{T} (GeV/c)");
     outputContainer->Add(fhPtInCone) ;
     
+    fhFRConeSumPt  = new TH2F
+    ("hFRConePtSum","#Sigma p_{T} in the froward region isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+    fhFRConeSumPt->SetYTitle("#Sigma p_{T}");
+    fhFRConeSumPt->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhFRConeSumPt) ;
+    
+    fhPtInFRCone  = new TH2F
+    ("hPtInFRCone","p_{T} in froward region isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
+    fhPtInFRCone->SetYTitle("p_{T in cone} (GeV/c)");
+    fhPtInFRCone->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhPtInFRCone) ;    
+    
     fhPtIso  = new TH1F("hPt","Isolated Number of particles",nptbins,ptmin,ptmax); 
     fhPtIso->SetYTitle("N");
     fhPtIso->SetXTitle("p_{T}(GeV/c)");
@@ -481,54 +377,55 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
   }
   
   if(fMakeSeveralIC){
-               char name[128];
-               char title[128];
+    const Int_t buffersize = 255;
+               char name[buffersize];
+               char title[buffersize];
                for(Int_t icone = 0; icone<fNCones; icone++){
-                 sprintf(name,"hPtSum_Cone_%d",icone);
-                 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                 snprintf(name, buffersize,"hPtSum_Cone_%d",icone);
+                 snprintf(title, buffersize,"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,"hPtSumPrompt_Cone_%d",icone);
-                   sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);
+                   snprintf(title, buffersize,"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,"hPtSumFragmentation_Cone_%d",icone);
-                   sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);
+                   snprintf(title, buffersize,"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,"hPtSumPi0Decay_Cone_%d",icone);
-                   sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);
+                   snprintf(title, buffersize,"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,"hPtSumOtherDecay_Cone_%d",icone);
-                   sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);
+                   snprintf(title, buffersize,"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,"hPtSumConversion_Cone_%d",icone);
-                   sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);
+                   snprintf(title, buffersize,"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,"hPtSumUnknown_Cone_%d",icone);
-                   sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+                   snprintf(name, buffersize,"hPtSumUnknown_Cone_%d",icone);
+                   snprintf(title, buffersize,"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)");
@@ -537,87 +434,87 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
                  }//Histos with MC
                  
                  for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){ 
-                   sprintf(name,"hPtThres_Cone_%d_Pt%d",icone,ipt);
-                   sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                   snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
+                   snprintf(title, buffersize,"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,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
-                   sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+                   snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
+                   snprintf(title, buffersize,"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,"hPtThresMCPrompt_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);
+                     snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtFracMCPrompt_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);
+                     snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtThresMCFragmentation_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);
+                     snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtFracMCFragmentation_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);
+                     snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtThresMCPi0Decay_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);
+                     snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtFracMCPi0Decay_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);
+                     snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtThresMCOtherDecay_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);
+                     snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtFracMCOtherDecay_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);
+                     snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtThresMCConversion_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);
+                     snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtFracMCConversion_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);
+                     snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtThresMCUnknown_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);
+                     snprintf(name, buffersize,"hPtThresMCUnknown_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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,"hPtFracMCUnknown_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);
+                     snprintf(name, buffersize,"hPtFracMCUnknown_Cone_%d_Pt%d",icone,ipt);
+                     snprintf(title, buffersize,"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]) ;  
@@ -634,53 +531,9 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
+       delete nmsHistos;
   }
   
-  //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 ;
   
 }
@@ -704,16 +557,13 @@ void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
   Int_t n = 0, nfrac = 0;
   Bool_t isolated = kFALSE ; 
   Float_t coneptsum = 0 ;
-  TObjArray * pl = new TObjArray()
+  TObjArray * pl = 0x0; 
   
   //Select the calorimeter for candidate isolation with neutral particles
   if(fCalorimeter == "PHOS")
-    pl = GetAODPHOS();
+    pl = GetPHOSClusters();
   else if (fCalorimeter == "EMCAL")
-    pl = GetAODEMCAL();
-  
-  //Get vertex for photon momentum calculation
-  if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(fVertex);
+    pl = GetEMCALClusters();
   
   //Loop on AOD branch, filled previously in AliAnaPhoton
   TLorentzVector mom ;
@@ -725,6 +575,17 @@ void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
     //If too small or too large pt, skip
     if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ; 
     
+    //check if it is low pt trigger particle, then adjust the isolation method
+    if(aodinput->Pt() < GetIsolationCut()->GetPtThreshold() || aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold()) 
+      continue ; //trigger should not from underlying event
+ //   if(GetIsolationCut()->GetICMethod()!=AliIsolationCut::kPtThresIC && aodinput->Pt()!=0.){ 
+//      //low pt trigger use pt threshold IC instead of other method
+//      if (aodinput->Pt()*(GetIsolationCut()->GetPtFraction())<GetIsolationCut()->GetPtThreshold() || aodinput->Pt()*(GetIsolationCut()->GetPtFraction())<GetIsolationCut()->GetSumPtThreshold()) {
+//       // printf("change the IC method to PtThresIC\n") ;
+//        GetIsolationCut()->SetICMethod(AliIsolationCut::kPtThresIC); 
+//      }
+//    }
+    
     //Check invariant mass, if pi0, skip.
     Bool_t decay = kFALSE ;
     if(fMakeInvMass) decay = CheckInvMass(iaod,aodinput);
@@ -732,7 +593,7 @@ void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
 
     //After cuts, study isolation
     n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
-    GetIsolationCut()->MakeIsolationCut(GetAODCTS(),pl,fVertex, kTRUE, aodinput, GetAODObjArrayName(), n,nfrac,coneptsum, isolated);
+    GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,GetReader(), kTRUE, aodinput, GetAODObjArrayName(), n,nfrac,coneptsum, isolated);
     aodinput->SetIsolated(isolated);
     if(GetDebug() > 1 && isolated) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",iaod);  
          
@@ -749,22 +610,34 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
   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("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
+  
+  //Get vertex for photon momentum calculation
+  Double_t vertex[]={0,0,0} ; //vertex ;
+  //Double_t vertex2[]={0,0,0} ; //vertex ;
+  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
+         GetReader()->GetVertex(vertex);
+         //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
+  }    
+       
   for(Int_t iaod = 0; iaod < naod ; iaod++){
     AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
-
-    Bool_t isolation   = aod->IsIsolated();              
+    
+    Bool_t  isolation  = aod->IsIsolated();              
     Float_t ptcluster  = aod->Pt();
     Float_t phicluster = aod->Phi();
     Float_t etacluster = aod->Eta();
+   // Float_t phiForward = aod->Phi()+TMath::PiOver2() ;
+    Float_t conesize   = GetIsolationCut()->GetConeSize();
+    
     //Recover reference arrays with clusters and tracks
-       TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
-    TObjArray * reftracks = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
-  
+    TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
+    TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
+    //If too small or too large pt, skip
+    if(aod->Pt() < GetMinPt() || aod->Pt() > GetMaxPt() ) continue ; 
+    
     if(fMakeSeveralIC) {
       //Analysis of multiple IC at same time
       MakeSeveralICAnalysis(aod);
@@ -773,83 +646,114 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
     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(reftracks, refclusters, fVertex, kFALSE, aod, "", n,nfrac,coneptsum, isolated);
+      GetIsolationCut()->MakeIsolationCut(reftracks, refclusters, GetReader(), kFALSE, aod, "", n,nfrac,coneptsum, isolated);
       fhConeSumPt->Fill(ptcluster,coneptsum);    
-         if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);    
-       }
+      if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);    
+    }
     
     //Fill pt distribution of particles in cone
     //Tracks
-    coneptsum=0;
-       if(reftracks){  
-        for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++){
-            AliAODTrack* track = (AliAODTrack *) reftracks->At(itrack);
-            fhPtInCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
-            coneptsum+=track->Pt();
-        }
-       }
-         
+    coneptsum = 0;
+    Double_t sumptFR = 0. ;
+    TObjArray * trackList   = GetCTSTracks() ;
+    for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++){
+      AliVTrack* track = (AliVTrack *) trackList->At(itrack);
+      //fill the histograms at forward range
+      if(!track){
+        printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
+        continue;
+      }
+      Double_t dPhi = phicluster - track->Phi() + TMath::PiOver2();
+      Double_t dEta = etacluster - track->Eta();
+      Double_t arg  = dPhi*dPhi + dEta*dEta;
+      if(TMath::Sqrt(arg) < conesize){
+        fhPtInFRCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
+        sumptFR+=track->Pt();
+      }    
+      dPhi = phicluster - track->Phi() - TMath::PiOver2();
+      arg  = dPhi*dPhi + dEta*dEta;
+      if(TMath::Sqrt(arg) < conesize){
+        fhPtInFRCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
+        sumptFR+=track->Pt();
+      }      
+    }
+    fhFRConeSumPt->Fill(ptcluster,sumptFR);
+    if(reftracks){  
+      for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++){
+        AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
+        fhPtInCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
+        coneptsum+=track->Pt();
+      }
+    }
+
     //CaloClusters
-       if(refclusters){    
-               TLorentzVector mom ;
-        for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++){
-            AliAODCaloCluster* calo = (AliAODCaloCluster *) refclusters->At(icalo);
-            calo->GetMomentum(mom,fVertex);//Assume that come from vertex in straight line
-            fhPtInCone->Fill(ptcluster, mom.Pt());
-            coneptsum+=mom.Pt();
-       }
+    if(refclusters){    
+      TLorentzVector mom ;
+      for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++){
+        AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
+        Int_t input = 0;
+        //                     if     (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo) input = 1 ;
+        //                     else if(fCalorimeter == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= icalo) input = 1;
+        
+        //Get Momentum vector, 
+        if     (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
+        //                     else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  
+        
+        fhPtInCone->Fill(ptcluster, mom.Pt());
+        coneptsum+=mom.Pt();
+      }
     }
          
-       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);
+    if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);
     
-       if(!fReMakeIC) fhConeSumPt->Fill(ptcluster,coneptsum);
+    if(!fReMakeIC) fhConeSumPt->Fill(ptcluster,coneptsum);
     
     if(isolation){    
-               
-         if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED, fill histograms\n", iaod);
-               
+      
+      if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED, fill histograms\n", iaod);
+      
       fhPtIso  ->Fill(ptcluster);
       fhPhiIso ->Fill(ptcluster,phicluster);
       fhEtaIso ->Fill(ptcluster,etacluster);
       
       if(IsDataMC()){
-       Int_t tag =aod->GetTag();
-       
-       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
-         fhPtIsoPrompt  ->Fill(ptcluster);
-         fhPhiIsoPrompt ->Fill(ptcluster,phicluster);
-         fhEtaIsoPrompt ->Fill(ptcluster,etacluster);
-       }
-       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
-         {
-           fhPtIsoFragmentation  ->Fill(ptcluster);
-           fhPhiIsoFragmentation ->Fill(ptcluster,phicluster);
-           fhEtaIsoFragmentation ->Fill(ptcluster,etacluster);
-         }
-       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
-         {
-           fhPtIsoPi0Decay  ->Fill(ptcluster);
-           fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster);
-           fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster);
-         }
-       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
-         {
-           fhPtIsoOtherDecay  ->Fill(ptcluster);
-           fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster);
-           fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster);
-         }
-       else if(GetMCAnalysisUtils()->CheckTagBit(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);
-         }
+        Int_t tag =aod->GetTag();
+        
+        if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
+          fhPtIsoPrompt  ->Fill(ptcluster);
+          fhPhiIsoPrompt ->Fill(ptcluster,phicluster);
+          fhEtaIsoPrompt ->Fill(ptcluster,etacluster);
+        }
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
+        {
+          fhPtIsoFragmentation  ->Fill(ptcluster);
+          fhPhiIsoFragmentation ->Fill(ptcluster,phicluster);
+          fhEtaIsoFragmentation ->Fill(ptcluster,etacluster);
+        }
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
+        {
+          fhPtIsoPi0Decay  ->Fill(ptcluster);
+          fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster);
+          fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster);
+        }
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
+        {
+          fhPtIsoOtherDecay  ->Fill(ptcluster);
+          fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster);
+          fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster);
+        }
+        else if(GetMCAnalysisUtils()->CheckTagBit(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
@@ -919,7 +823,7 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
       GetIsolationCut()->MakeIsolationCut(ph->GetObjArray(GetAODObjArrayName()+"Tracks"), 
                                  ph->GetObjArray(GetAODObjArrayName()+"Clusters"),
-                                 fVertex, kFALSE, ph, "",n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
+                                 GetReader(), kFALSE, ph, "",n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
 
       //Normal ptThreshold cut
       if(n[icone][ipt] == 0) {