]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
isolation task/cuts modifs (N. Arbor)
authorrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 25 Mar 2011 16:31:36 +0000 (16:31 +0000)
committerrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 25 Mar 2011 16:31:36 +0000 (16:31 +0000)
PWG4/PartCorrBase/AliIsolationCut.cxx
PWG4/PartCorrBase/AliIsolationCut.h
PWG4/PartCorrDep/AliAnaParticleIsolation.cxx
PWG4/PartCorrDep/AliAnaParticleIsolation.h

index cb470a4823ee739b97443b79ebc6a5961037b107..a3f96238bb9166c402f92799a689716e4dae36b2 100755 (executable)
-/**************************************************************************
- * 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 is 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:  $ */
-
-//_________________________________________________________________________
-// Class containing methods for the isolation cut. 
-// An AOD candidate (AliAODPWG4ParticleCorrelation type)
-// is passed. Look in a cone around the candidate and study
-// the hadronic activity inside to decide if the candidate is isolated
-//
-//
-//*-- Author: Gustavo Conesa (LNF-INFN) 
-
-//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
-//-Yaxian Mao (check the candidate particle is the leading particle or not at the same hemishere)
-
-//////////////////////////////////////////////////////////////////////////////
-  
-  
-// --- ROOT system --- 
-//#include <Riostream.h>
-#include <TLorentzVector.h>
-#include <TObjArray.h>
-
-// --- AliRoot system --- 
-#include "AliIsolationCut.h" 
-#include "AliAODPWG4ParticleCorrelation.h"
-#include "AliAODTrack.h"
-#include "AliVCluster.h"
-#include "AliCaloTrackReader.h"
-#include "AliMixedEvent.h"
-
-ClassImp(AliIsolationCut)
-  
-//____________________________________________________________________________
-  AliIsolationCut::AliIsolationCut() : 
-    TObject(),
-    fConeSize(0.),fPtThreshold(0.), fSumPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)
-{
-  //default ctor
-  
-  //Initialize parameters
-  InitParameters();
-
-}
-
-//____________________________________________________________________________
-TString AliIsolationCut::GetICParametersList()
-{
-  //Put data member values in string to keep in output container
-  
-  TString parList ; //this will be list of parameters used for this analysis.
-  const Int_t buffersize = 255;
-  char onePar[buffersize] ;
-  
-  snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
-  parList+=onePar ;    
-  snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
-  parList+=onePar ;
-  snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
-  parList+=onePar ;
-  snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
-  parList+=onePar ;
-  snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
-  parList+=onePar ;
-  snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
-  parList+=onePar ;
-
-  return parList; 
-}
-
-//____________________________________________________________________________
-void AliIsolationCut::InitParameters()
-{
-  //Initialize the parameters of the analysis.
-  
-  fConeSize       = 0.4 ; 
-  fPtThreshold    = 1.  ; 
-  fSumPtThreshold = 0.5 ; 
-  fPtFraction     = 0.1 ; 
-  fPartInCone     = kNeutralAndCharged;
-  fICMethod       = kPtThresIC; // 0 pt threshol method, 1 cone pt sum method
-  
-}
-
-//__________________________________________________________________
-void  AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS,  TObjArray * const plNe, AliCaloTrackReader * const reader, 
-                                       const Bool_t bFillAOD, AliAODPWG4ParticleCorrelation  *pCandidate, 
-                                       const TString & aodArrayRefName,
-                                       Int_t & n, Int_t & nfrac, Float_t &coneptsum,  Bool_t  &isolated) const
-{  
-  //Search in cone around a candidate particle if it is isolated 
-  Float_t phiC  = pCandidate->Phi() ;
-  Float_t etaC  = pCandidate->Eta() ;
-  Float_t ptC   = pCandidate->Pt() ;
-  Float_t pt    = -100. ;
-  Float_t eta   = -100. ;
-  Float_t phi   = -100. ;
-  Float_t rad   = -100. ;
-  
-  n         = 0 ;
-  nfrac     = 0 ;
-  coneptsum = 0.; 
-  isolated  = kFALSE;
-
-  //Initialize the array with refrences
-  TObjArray * refclusters = 0x0;
-  TObjArray * reftracks   = 0x0;
-  Int_t ntrackrefs   = 0;
-  Int_t nclusterrefs = 0;
-  //Check charged particles in cone.
-  if(plCTS && (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged)){
-    TVector3 p3;
-    for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
-      AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ; 
-      //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
-      if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1)) continue ;
-      p3.SetXYZ(track->Px(),track->Py(),track->Pz());
-      pt   = p3.Pt();
-      eta  = p3.Eta();
-      phi  = p3.Phi() ;
-      if(phi<0) phi+=TMath::TwoPi();
-      
-      //only loop the particle at the same side of candidate
-      if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
-      //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events
-      if(pt > ptC){
-        n         = -1;
-        nfrac     = -1;
-        coneptsum = -1;
-        isolated  = kFALSE;
-        if(bFillAOD && reftracks) {
-          reftracks->Clear(); 
-          delete reftracks;
-        }
-        return ;
-      }
-      //Check if there is any particle inside cone with pt larger than  fPtThreshold
-      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
-      
-      if(rad < fConeSize){
-        if(bFillAOD) {
-          ntrackrefs++;
-          if(ntrackrefs == 1){
-            reftracks = new TObjArray(0);
-            //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
-            TString tempo(aodArrayRefName)  ; 
-            tempo += "Tracks" ; 
-            reftracks->SetName(tempo);
-            reftracks->SetOwner(kFALSE);
-          }
-          reftracks->Add(track);
-        }
-        //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
-        coneptsum+=pt;
-        if(pt > fPtThreshold )    n++;
-        if(pt > fPtFraction*ptC ) nfrac++;  
-      } // Inside cone
-    }// charged particle loop
-  }//Tracks
-  
-  //Check neutral particles in cone.  
-  if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){
-         
-    //Get vertex for photon momentum calculation
-    //Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;
-    //if(reader->GetDataType()!= AliCaloTrackReader::kMC) 
-    //{
-      //if(reader->GetSecondInputAODTree()) reader->GetSecondInputAODVertex(vertex2);
-    //}
-    TLorentzVector mom ;
-    for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
-      AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
-      
-      //Get the index where the cluster comes, to retrieve the corresponding vertex
-      Int_t evtIndex = 0 ; 
-      if (reader->GetMixedEvent()) {
-        evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
-      }
-      
-      //Do not count the candidate (photon or pi0) or the daughters of the candidate
-      if(calo->GetID() == pCandidate->GetCaloLabel(0) || calo->GetID() == pCandidate->GetCaloLabel(1)) continue ;      //Skip matched clusters with tracks
-      
-      if(calo->GetNTracksMatched() > 0) continue ; 
-      
-      //Input from second AOD?
-      //Int_t input = 0;
-      //      if     (pCandidate->GetDetector() == "EMCAL" && reader->GetAODEMCALNormalInputEntries() <= ipr) input = 1 ;
-      //      else if(pCandidate->GetDetector() == "PHOS"  && reader->GetAODPHOSNormalInputEntries()  <= ipr) input = 1;
-      
-      //Get Momentum vector, 
-      //if     (input == 0) 
-      calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;//Assume that come from vertex in straight line
-      //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  
-      
-      pt   = mom.Pt();
-      eta  = mom.Eta();
-      phi  = mom.Phi() ;
-      if(phi<0) phi+=TMath::TwoPi();
-      //only loop the particle at the same side of candidate
-      
-      if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
-      //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events
-      if(pt > ptC){
-        n         = -1;
-        nfrac     = -1;
-        coneptsum = -1;
-        isolated  = kFALSE;
-        if(bFillAOD){
-          if(reftracks){  
-            reftracks  ->Clear();
-            delete reftracks;
-          }
-          if(refclusters){
-            refclusters->Clear(); 
-            delete refclusters;
-          }
-        }
-        return ;
-      }
-      
-      //Check if there is any particle inside cone with pt larger than  fPtThreshold
-      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
-      if(rad < fConeSize){
-        if(bFillAOD) {
-          nclusterrefs++;
-          if(nclusterrefs==1){
-            refclusters = new TObjArray(0);
-            //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
-            TString tempo(aodArrayRefName)  ; 
-            tempo += "Clusters" ; 
-            refclusters->SetName(tempo);
-            refclusters->SetOwner(kFALSE);
-          }
-          refclusters->Add(calo);
-        }
-        //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
-        coneptsum+=pt;
-        if(pt > fPtThreshold )     n++;
-        //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
-        if(fPtFraction*ptC<fPtThreshold) {
-            if(pt>fPtThreshold)    nfrac++ ;
-        }
-        else {
-            if(pt>fPtFraction*ptC) nfrac++; 
-        }
-      }//in cone
-    }// neutral particle loop
-  }//neutrals
-
-  //printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
-  
-  //Add reference arrays to AOD when filling AODs only
-  if(bFillAOD) {
-    if(refclusters)    pCandidate->AddObjArray(refclusters);
-    if(reftracks)        pCandidate->AddObjArray(reftracks);
-  }
-  //Check isolation, depending on method.
-  if( fICMethod == kPtThresIC){
-    if(n==0) isolated = kTRUE ;
-  }
-  else if( fICMethod == kSumPtIC){
-    if(coneptsum < fSumPtThreshold)
-      isolated  =  kTRUE ;
-  }
-  else if( fICMethod == kPtFracIC){
-    if(nfrac==0) isolated = kTRUE ;
-  }
-  else if( fICMethod == kSumPtFracIC){
-    //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
-    if(coneptsum < fPtFraction*ptC && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
-  }
-  
-}
-
-//__________________________________________________________________
-void AliIsolationCut::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("IC method          =     %d\n", fICMethod) ; 
-  printf("Cone Size          =     %1.2f\n", fConeSize) ; 
-  printf("pT threshold       =     %2.1f\n", fPtThreshold) ;
-  printf("pT fraction        =     %3.1f\n", fPtFraction) ;
-  printf("particle type in cone =  %d\n",fPartInCone);
-  printf("    \n") ;
-  
-} 
+/**************************************************************************\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 is 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:  $ */\r
+\r
+//_________________________________________________________________________\r
+// Class containing methods for the isolation cut. \r
+// An AOD candidate (AliAODPWG4ParticleCorrelation type)\r
+// is passed. Look in a cone around the candidate and study\r
+// the hadronic activity inside to decide if the candidate is isolated\r
+//\r
+//\r
+//*-- Author: Gustavo Conesa (LNF-INFN) \r
+\r
+//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)\r
+//-Yaxian Mao (check the candidate particle is the leading particle or not at the same hemishere)\r
+\r
+//////////////////////////////////////////////////////////////////////////////\r
+  \r
+  \r
+// --- ROOT system --- \r
+//#include <Riostream.h>\r
+#include <TLorentzVector.h>\r
+#include <TObjArray.h>\r
+\r
+// --- AliRoot system --- \r
+#include "AliIsolationCut.h" \r
+#include "AliAODPWG4ParticleCorrelation.h"\r
+#include "AliAODTrack.h"\r
+#include "AliVCluster.h"\r
+#include "AliCaloTrackReader.h"\r
+#include "AliMixedEvent.h"\r
+\r
+ClassImp(AliIsolationCut)\r
+  \r
+//____________________________________________________________________________\r
+  AliIsolationCut::AliIsolationCut() : \r
+    TObject(),\r
+    fConeSize(0.),fPtThreshold(0.), fSumPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)\r
\r
+{\r
+  //default ctor\r
+  \r
+  //Initialize parameters\r
+  InitParameters();\r
+\r
+}\r
+\r
+//____________________________________________________________________________\r
+TString AliIsolationCut::GetICParametersList()\r
+{\r
+  //Put data member values in string to keep in output container\r
+  \r
+  TString parList ; //this will be list of parameters used for this analysis.\r
+  const Int_t buffersize = 255;\r
+  char onePar[buffersize] ;\r
+  \r
+  snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;\r
+  parList+=onePar ;    \r
+  snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;\r
+  parList+=onePar ;\r
+  snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;\r
+  parList+=onePar ;\r
+  snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;\r
+  parList+=onePar ;\r
+  snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;\r
+  parList+=onePar ;\r
+  snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;\r
+  parList+=onePar ;\r
+\r
+  return parList; \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliIsolationCut::InitParameters()\r
+{\r
+  //Initialize the parameters of the analysis.\r
+  \r
+  fConeSize       = 0.4 ; \r
+  fPtThreshold    = 1.  ; \r
+  fSumPtThreshold = 0.5 ; \r
+  fPtFraction     = 0.1 ; \r
+  fPartInCone     = kOnlyCharged;\r
+  fICMethod       = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method\r
+  \r
+}\r
+\r
+//__________________________________________________________________\r
+void  AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS,  TObjArray * const plNe, AliCaloTrackReader * const reader, \r
+                                       const Bool_t bFillAOD, AliAODPWG4ParticleCorrelation  *pCandidate, \r
+                                       const TString & aodArrayRefName,\r
+                                       Int_t & n, Int_t & nfrac, Float_t &coneptsum,  Bool_t  &isolated, Bool_t &leading) const\r
+{  \r
+  //Search in cone around a candidate particle if it is isolated \r
+  Float_t phiC  = pCandidate->Phi() ;\r
+  Float_t etaC  = pCandidate->Eta() ;\r
+  Float_t ptC   = pCandidate->Pt() ;\r
+  Float_t pt    = -100. ;\r
+  Float_t eta   = -100. ;\r
+  Float_t phi   = -100. ;\r
+  Float_t rad   = -100. ;\r
+  \r
+  n         = 0 ;\r
+  nfrac     = 0 ;\r
+  coneptsum = 0.; \r
+  isolated  = kFALSE;\r
+  leading = kTRUE;\r
+\r
+  //Initialize the array with refrences\r
+  TObjArray * refclusters = 0x0;\r
+  TObjArray * reftracks   = 0x0;\r
+  Int_t ntrackrefs   = 0;\r
+  Int_t nclusterrefs = 0;\r
+  //Check charged particles in cone.\r
+  if(plCTS && (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged)){\r
+    TVector3 p3;\r
+    for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){\r
+      AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ; \r
+      //Do not count the candidate (pion, conversion photon) or the daughters of the candidate\r
+      if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1)) continue ;\r
+      p3.SetXYZ(track->Px(),track->Py(),track->Pz());\r
+      pt   = p3.Pt();\r
+      eta  = p3.Eta();\r
+      phi  = p3.Phi() ;\r
+      if(phi<0) phi+=TMath::TwoPi();\r
+      \r
+      //only loop the particle at the same side of candidate\r
+      if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;\r
+      //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events\r
+      if(pt > ptC){\r
+        n         = -1;\r
+        nfrac     = -1;\r
+        coneptsum = -1;\r
+        isolated  = kFALSE;\r
+       leading = kFALSE;\r
+        if(bFillAOD && reftracks) {\r
+          reftracks->Clear(); \r
+          delete reftracks;\r
+        }\r
+        return ;\r
+      }\r
+      //Check if there is any particle inside cone with pt larger than  fPtThreshold\r
+      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));\r
+      \r
+      if(rad < fConeSize){\r
+        if(bFillAOD) {\r
+          ntrackrefs++;\r
+          if(ntrackrefs == 1){\r
+            reftracks = new TObjArray(0);\r
+            //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));\r
+            TString tempo(aodArrayRefName)  ; \r
+            tempo += "Tracks" ; \r
+            reftracks->SetName(tempo);\r
+            reftracks->SetOwner(kFALSE);\r
+          }\r
+          reftracks->Add(track);\r
+        }\r
+        //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);\r
+        coneptsum+=pt;\r
+        if(pt > fPtThreshold )    n++;\r
+        if(pt > fPtFraction*ptC ) nfrac++;  \r
+      } // Inside cone\r
+    }// charged particle loop\r
+  }//Tracks\r
+  \r
+  //Check neutral particles in cone.  \r
+  if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){\r
+         \r
+    //Get vertex for photon momentum calculation\r
+    //Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;\r
+    //if(reader->GetDataType()!= AliCaloTrackReader::kMC) \r
+    //{\r
+      //if(reader->GetSecondInputAODTree()) reader->GetSecondInputAODVertex(vertex2);\r
+    //}\r
+    TLorentzVector mom ;\r
+    for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){\r
+      AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;\r
+      \r
+      //Get the index where the cluster comes, to retrieve the corresponding vertex\r
+      Int_t evtIndex = 0 ; \r
+      if (reader->GetMixedEvent()) {\r
+        evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; \r
+      }\r
+      \r
+      //Do not count the candidate (photon or pi0) or the daughters of the candidate\r
+      if(calo->GetID() == pCandidate->GetCaloLabel(0) || calo->GetID() == pCandidate->GetCaloLabel(1)) continue ;      //Skip matched clusters with tracks\r
+      \r
+      if(calo->GetNTracksMatched() > 0) continue ; \r
+      \r
+      //Input from second AOD?\r
+      //Int_t input = 0;\r
+      //      if     (pCandidate->GetDetector() == "EMCAL" && reader->GetAODEMCALNormalInputEntries() <= ipr) input = 1 ;\r
+      //      else if(pCandidate->GetDetector() == "PHOS"  && reader->GetAODPHOSNormalInputEntries()  <= ipr) input = 1;\r
+      \r
+      //Get Momentum vector, \r
+      //if     (input == 0) \r
+      calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;//Assume that come from vertex in straight line\r
+      //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  \r
+      \r
+      pt   = mom.Pt();\r
+      eta  = mom.Eta();\r
+      phi  = mom.Phi() ;\r
+      if(phi<0) phi+=TMath::TwoPi();\r
+      //only loop the particle at the same side of candidate\r
+      \r
+      if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;\r
+      //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events\r
+      if(pt > ptC){\r
+        n         = -1;\r
+        nfrac     = -1;\r
+        coneptsum = -1;\r
+        isolated  = kFALSE;\r
+       leading = kFALSE;\r
+        if(bFillAOD){\r
+          if(reftracks){  \r
+            reftracks  ->Clear();\r
+            delete reftracks;\r
+          }\r
+          if(refclusters){\r
+            refclusters->Clear(); \r
+            delete refclusters;\r
+          }\r
+        }\r
+        return ;\r
+      }\r
+      \r
+      //Check if there is any particle inside cone with pt larger than  fPtThreshold\r
+      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));\r
+      if(rad < fConeSize){\r
+        if(bFillAOD) {\r
+          nclusterrefs++;\r
+          if(nclusterrefs==1){\r
+            refclusters = new TObjArray(0);\r
+            //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));\r
+            TString tempo(aodArrayRefName)  ; \r
+            tempo += "Clusters" ; \r
+            refclusters->SetName(tempo);\r
+            refclusters->SetOwner(kFALSE);\r
+          }\r
+          refclusters->Add(calo);\r
+        }\r
+        //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);\r
+        coneptsum+=pt;\r
+        if(pt > fPtThreshold )     n++;\r
+        //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly\r
+        if(fPtFraction*ptC<fPtThreshold) {\r
+            if(pt>fPtThreshold)    nfrac++ ;\r
+        }\r
+        else {\r
+            if(pt>fPtFraction*ptC) nfrac++; \r
+        }\r
+      }//in cone\r
+    }// neutral particle loop\r
+  }//neutrals\r
+\r
+  //printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);\r
+  \r
+  //Add reference arrays to AOD when filling AODs only\r
+  if(bFillAOD) {\r
+    if(refclusters)    pCandidate->AddObjArray(refclusters);\r
+    if(reftracks)        pCandidate->AddObjArray(reftracks);\r
+  }\r
+  //Check isolation, depending on method.\r
+  if( fICMethod == kPtThresIC){\r
+    if(n==0) isolated = kTRUE ;\r
+  }\r
+  else if( fICMethod == kSumPtIC){\r
+    if(coneptsum < fSumPtThreshold)\r
+      isolated  =  kTRUE ;\r
+  }\r
+  else if( fICMethod == kPtFracIC){\r
+    if(nfrac==0) isolated = kTRUE ;\r
+  }\r
+  else if( fICMethod == kSumPtFracIC){\r
+    //when the fPtFraction*ptC < fSumPtThreshold then consider the later case\r
+    if(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;\r
+    if(fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;\r
+  }\r
+  \r
+}\r
+\r
+//__________________________________________________________________\r
+void AliIsolationCut::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("IC method          =     %d\n", fICMethod) ; \r
+  printf("Cone Size          =     %1.2f\n", fConeSize) ; \r
+  printf("pT threshold       =     %2.1f\n", fPtThreshold) ;\r
+  printf("pT fraction        =     %3.1f\n", fPtFraction) ;\r
+  printf("particle type in cone =  %d\n",fPartInCone);\r
+  printf("    \n") ;\r
+  \r
+} \r
index 65783f97f96724985127174d885654ca82c7633b..38ee9e63d54cc46ce0d9a05593d0c5adcbd0a654 100755 (executable)
@@ -1,89 +1,89 @@
-#ifndef ALIISOLATIONCUT_H
-#define ALIISOLATIONCUT_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice     */
-/* $Id:  */
-
-//_________________________________________________________________________
-// Class containing methods for the isolation cut. 
-// An AOD candidate (AliAODPWG4ParticleCorrelation type)
-// is passed. Look in a cone around the candidate and study
-// the hadronic activity inside to decide if the candidate is isolated
-//
-//
-// 
-// -- Author: Gustavo Conesa (INFN-LNF)
-
-//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
-
-// --- ROOT system --- 
-#include <TObject.h>
-class TObjArray ;
-
-// --- ANALYSIS system ---
-class AliAODPWG4ParticleCorrelation ;
-class AliCaloTrackReader ;
-
-class AliIsolationCut : public TObject {
-  
- public: 
-  AliIsolationCut() ; // default ctor
-  virtual ~AliIsolationCut() {;} //virtual dtalr
- private:
-  AliIsolationCut(const AliIsolationCut & g) ; // cpy ctor
-  AliIsolationCut & operator = (const AliIsolationCut & g) ;//cpy assignment
-
- public: 
-  enum type {kPtThresIC, kSumPtIC, kPtFracIC, kSumPtFracIC};
-  enum partInCone {kNeutralAndCharged=0, kOnlyNeutral=1, kOnlyCharged=2};
-       
-  Float_t    GetConeSize()           const {return fConeSize   ;}
-  Float_t    GetPtThreshold()        const {return fPtThreshold;}
-  Float_t    GetSumPtThreshold()     const {return fSumPtThreshold;}
-  Float_t    GetPtFraction()         const {return fPtFraction ;}
-  Int_t      GetICMethod()           const {return fICMethod   ;}
-  Int_t      GetParticleTypeInCone() const {return fPartInCone ;}
-       
-  TString    GetICParametersList() ; 
-
-  void MakeIsolationCut(TObjArray * const plCTS, TObjArray * const plNe, AliCaloTrackReader * const reader, 
-                       const Bool_t bFillAOD, AliAODPWG4ParticleCorrelation  * pCandidate, const TString &aodObjArrayName,
-                       Int_t &n, Int_t & nfrac, Float_t &ptsum, Bool_t & isolated) const ;  
-  
-  void Print(const Option_t * opt)const;
-  
-  void SetConeSize(Float_t r)        {fConeSize = r ; }
-  void SetPtThreshold(Float_t pt)    {fPtThreshold = pt; }
-  void SetSumPtThreshold(Float_t sumpt)    {fSumPtThreshold = sumpt; }
-  void SetPtFraction(Float_t pt)     {fPtFraction = pt; }
-  void SetICMethod(Int_t i )         {fICMethod = i ; }
-  void SetParticleTypeInCone(Int_t i){fPartInCone = i;}
-       
-  void InitParameters();
-  
-  
- private:
-  
-  Float_t      fConeSize ;    //Size of the isolation cone 
-  Float_t      fPtThreshold ; //Mimium pt of the particles in the cone or sum in cone (UE pt mean in the forward region cone)
-  Float_t      fSumPtThreshold ; //Minium of sum pt of the particles in the cone (UE sum in the forward region cone)
-  Float_t      fPtFraction ;  //Fraction of the momentum of particles in cone or sum in cone
-  Int_t        fICMethod ;    //Isolation cut method to be used
-                              // kPtIC: Pt threshold method
-                              // kSumPtIC: Cone pt sum method
-                              // kPtFracIC:   Pt threshold, fraction of candidate pt, method
-                              // kSumPtFracIC:   Cone pt sum , fraction of cone sum, method
-  Int_t        fPartInCone;   // Type of particles inside cone:
-                                                         // kNeutralAndCharged, kOnlyNeutral, kOnlyCharged
-                                
-       
-       
-  ClassDef(AliIsolationCut,3)
-} ;
-
-
-#endif //ALIISOLATIONCUT_H
-
-
-
+#ifndef ALIISOLATIONCUT_H\r
+#define ALIISOLATIONCUT_H\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice     */\r
+/* $Id:  */\r
+\r
+//_________________________________________________________________________\r
+// Class containing methods for the isolation cut. \r
+// An AOD candidate (AliAODPWG4ParticleCorrelation type)\r
+// is passed. Look in a cone around the candidate and study\r
+// the hadronic activity inside to decide if the candidate is isolated\r
+//\r
+//\r
+// \r
+// -- Author: Gustavo Conesa (INFN-LNF)\r
+\r
+//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)\r
+\r
+// --- ROOT system --- \r
+#include <TObject.h>\r
+class TObjArray ;\r
+\r
+// --- ANALYSIS system ---\r
+class AliAODPWG4ParticleCorrelation ;\r
+class AliCaloTrackReader ;\r
+\r
+class AliIsolationCut : public TObject {\r
+  \r
+ public: \r
+  AliIsolationCut() ; // default ctor\r
+  virtual ~AliIsolationCut() {;} //virtual dtalr\r
+ private:\r
+  AliIsolationCut(const AliIsolationCut & g) ; // cpy ctor\r
+  AliIsolationCut & operator = (const AliIsolationCut & g) ;//cpy assignment\r
+\r
+ public: \r
\r
+  enum type {kPtThresIC, kSumPtIC, kPtFracIC, kSumPtFracIC};\r
+  enum partInCone {kNeutralAndCharged=0, kOnlyNeutral=1, kOnlyCharged=2};\r
+       \r
+  Float_t    GetConeSize()           const {return fConeSize   ;}\r
+  Float_t    GetPtThreshold()        const {return fPtThreshold;}\r
+  Float_t    GetSumPtThreshold()     const {return fSumPtThreshold;}\r
+  Float_t    GetPtFraction()         const {return fPtFraction ;}\r
+  Int_t      GetICMethod()           const {return fICMethod   ;}\r
+  Int_t      GetParticleTypeInCone() const {return fPartInCone ;}\r
+       \r
+  TString    GetICParametersList() ; \r
+\r
+  void MakeIsolationCut(TObjArray * const plCTS, TObjArray * const plNe, AliCaloTrackReader * const reader, \r
+                       const Bool_t bFillAOD, AliAODPWG4ParticleCorrelation  * pCandidate, const TString &aodObjArrayName,\r
+                       Int_t &n, Int_t & nfrac, Float_t &ptsum, Bool_t & isolated, Bool_t & leading) const ;  \r
+  \r
+  void Print(const Option_t * opt)const;\r
+  \r
+  void SetConeSize(Float_t r)        {fConeSize = r ; }\r
+  void SetPtThreshold(Float_t pt)    {fPtThreshold = pt; }\r
+  void SetSumPtThreshold(Float_t sumpt)    {fSumPtThreshold = sumpt; }\r
+  void SetPtFraction(Float_t pt)     {fPtFraction = pt; }\r
+  void SetICMethod(Int_t i )         {fICMethod = i ; }\r
+  void SetParticleTypeInCone(Int_t i){fPartInCone = i;}\r
+       \r
+  void InitParameters();\r
+  \r
+  \r
+ private:\r
+  \r
+  Float_t      fConeSize ;    //Size of the isolation cone \r
+  Float_t      fPtThreshold ; //Mimium pt of the particles in the cone or sum in cone (UE pt mean in the forward region cone)\r
+  Float_t      fSumPtThreshold ; //Minium of sum pt of the particles in the cone (UE sum in the forward region cone)\r
+  Float_t      fPtFraction ;  //Fraction of the momentum of particles in cone or sum in cone\r
+  Int_t        fICMethod ;    //Isolation cut method to be used\r
+                              // kPtIC: Pt threshold method\r
+                              // kSumPtIC: Cone pt sum method\r
+                              // kPtFracIC:   Pt threshold, fraction of candidate pt, method\r
+                              // kSumPtFracIC:   Cone pt sum , fraction of cone sum, method\r
+  Int_t        fPartInCone;   // Type of particles inside cone:\r
+                                                         // kNeutralAndCharged, kOnlyNeutral, kOnlyCharged\r
+                                \r
+       \r
+       \r
+  ClassDef(AliIsolationCut,3)\r
+} ;\r
+\r
+\r
+#endif //ALIISOLATIONCUT_H\r
+\r
+\r
+\r
index 0a8bc03b021afbdc42c66f4f27d65d5fe85c137b..90d6db0af952f6b48f90caae419885fe63312525 100755 (executable)
-/**************************************************************************
- * 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) 
-
-//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
-//////////////////////////////////////////////////////////////////////////////
-  
-  
-// --- ROOT system --- 
-#include <TClonesArray.h>
-#include <TList.h>
-#include <TObjString.h>
-#include <TH2F.h>
-//#include <Riostream.h>
-#include <TClass.h>
-
-// --- Analysis system --- 
-#include "AliAnaParticleIsolation.h" 
-#include "AliCaloTrackReader.h"
-#include "AliIsolationCut.h"
-#include "AliNeutralMesonSelection.h"
-#include "AliAODPWG4ParticleCorrelation.h"
-#include "AliMCAnalysisUtils.h"
-#include "AliVTrack.h"
-#include "AliVCluster.h"
-
-ClassImp(AliAnaParticleIsolation)
-  
-//____________________________________________________________________________
-  AliAnaParticleIsolation::AliAnaParticleIsolation() : 
-    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
-    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();
-       
-  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() 
-{
-  //dtor
-  //do not delete histograms
-  
-  //delete [] fConeSizes ; 
-  //delete [] fPtThresholds ; 
-  //delete [] fPtFractions ; 
-
-}
-
-//_________________________________________________________________________
-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
-  // 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("AliAnaParticleIsolation::CheckInvMass() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
-      return kTRUE ;
-    }
-  }//loop
-  
-  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()
-{  
-  // Create histograms to be saved in output file and 
-  // store them in outputContainer
-  TList * outputContainer = new TList() ; 
-  outputContainer->SetName("IsolatedParticleHistos") ; 
-  
-  Int_t nptbins  = GetHistoPtBins();
-  Int_t nphibins = GetHistoPhiBins();
-  Int_t netabins = GetHistoEtaBins();
-  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 isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
-    fhConeSumPt->SetYTitle("#Sigma p_{T}");
-    fhConeSumPt->SetXTitle("p_{T} (GeV/c)");
-    outputContainer->Add(fhConeSumPt) ;
-    
-    fhPtInCone  = new TH2F
-      ("hPtInCone","p_{T} in isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
-    fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
-    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)");
-    outputContainer->Add(fhPtIso) ; 
-    
-    fhPhiIso  = new TH2F
-      ("hPhi","Isolated Number of particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-    fhPhiIso->SetYTitle("#phi");
-    fhPhiIso->SetXTitle("p_{T} (GeV/c)");
-    outputContainer->Add(fhPhiIso) ; 
-    
-    fhEtaIso  = new TH2F
-      ("hEta","Isolated Number of particles",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-    fhEtaIso->SetYTitle("#eta");
-    fhEtaIso->SetXTitle("p_{T} (GeV/c)");
-    outputContainer->Add(fhEtaIso) ;
-    
-    if(IsDataMC()){
-      
-      fhPtIsoPrompt  = new TH1F("hPtMCPrompt","Isolated Number of #gamma prompt",nptbins,ptmin,ptmax); 
-      fhPtIsoPrompt->SetYTitle("N");
-      fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
-      outputContainer->Add(fhPtIsoPrompt) ; 
-      
-      fhPhiIsoPrompt  = new TH2F
-       ("hPhiMCPrompt","Isolated Number of #gamma prompt",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiIsoPrompt->SetYTitle("#phi");
-      fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhPhiIsoPrompt) ; 
-      
-      fhEtaIsoPrompt  = new TH2F
-       ("hEtaMCPrompt","Isolated Number of #gamma prompt ",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaIsoPrompt->SetYTitle("#eta");
-      fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhEtaIsoPrompt) ;
-      
-      fhPtIsoFragmentation  = new TH1F("hPtMCFragmentation","Isolated Number of #gamma",nptbins,ptmin,ptmax); 
-      fhPtIsoFragmentation->SetYTitle("N");
-      fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
-      outputContainer->Add(fhPtIsoFragmentation) ; 
-      
-      fhPhiIsoFragmentation  = new TH2F
-       ("hPhiMCFragmentation","Isolated Number of #gamma fragmentation",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiIsoFragmentation->SetYTitle("#phi");
-      fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhPhiIsoFragmentation) ; 
-      
-      fhEtaIsoFragmentation  = new TH2F
-       ("hEtaMCFragmentation","Isolated Number of #gamma fragmentation",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaIsoFragmentation->SetYTitle("#eta");
-      fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhEtaIsoFragmentation) ;
-      
-      fhPtIsoPi0Decay  = new TH1F("hPtMCPi0Decay","Isolated Number of #gamma from #pi^{0} decay",nptbins,ptmin,ptmax); 
-      fhPtIsoPi0Decay->SetYTitle("N");
-      fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
-      outputContainer->Add(fhPtIsoPi0Decay) ; 
-      
-      fhPhiIsoPi0Decay  = new TH2F
-       ("hPhiMCPi0Decay","Isolated Number of #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiIsoPi0Decay->SetYTitle("#phi");
-      fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhPhiIsoPi0Decay) ; 
-      
-      fhEtaIsoPi0Decay  = new TH2F
-       ("hEtaMCPi0Decay","Isolated Number of #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaIsoPi0Decay->SetYTitle("#eta");
-      fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhEtaIsoPi0Decay) ;
-      
-      fhPtIsoOtherDecay  = new TH1F("hPtMCOtherDecay","Isolated Number of #gamma from non #pi^{0} decay",nptbins,ptmin,ptmax); 
-      fhPtIsoOtherDecay->SetYTitle("N");
-      fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
-      outputContainer->Add(fhPtIsoOtherDecay) ; 
-      
-      fhPhiIsoOtherDecay  = new TH2F
-       ("hPhiMCOtherDecay","Isolated Number of #gamma from non #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiIsoOtherDecay->SetYTitle("#phi");
-      fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhPhiIsoOtherDecay) ; 
-      
-      fhEtaIsoOtherDecay  = new TH2F
-       ("hEtaMCOtherDecay","Isolated Number of #gamma non #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaIsoOtherDecay->SetYTitle("#eta");
-      fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhEtaIsoOtherDecay) ;
-      
-      fhPtIsoConversion  = new TH1F("hPtMCConversion","Isolated Number of #gamma converted",nptbins,ptmin,ptmax); 
-      fhPtIsoConversion->SetYTitle("N");
-      fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
-      outputContainer->Add(fhPtIsoConversion) ; 
-      
-      fhPhiIsoConversion  = new TH2F
-       ("hPhiMCConversion","Isolated Number of #gamma converted",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiIsoConversion->SetYTitle("#phi");
-      fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhPhiIsoConversion) ; 
-      
-      fhEtaIsoConversion  = new TH2F
-       ("hEtaMCConversion","Isolated Number of #gamma converted",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaIsoConversion->SetYTitle("#eta");
-      fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
-      outputContainer->Add(fhEtaIsoConversion) ;
-      
-      fhPtIsoUnknown  = new TH1F("hPtMCUnknown","Isolated Number of non #gamma particles",nptbins,ptmin,ptmax); 
-      fhPtIsoUnknown->SetYTitle("N");
-      fhPtIsoUnknown->SetXTitle("p_{T}(GeV/c)");
-      outputContainer->Add(fhPtIsoUnknown) ; 
-      
-      fhPhiIsoUnknown  = new TH2F
-       ("hPhiMCUnknown","Isolated Number of non #gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiIsoUnknown->SetYTitle("#phi");
-      fhPhiIsoUnknown->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhPhiIsoUnknown) ; 
-      
-      fhEtaIsoUnknown  = new TH2F
-       ("hEtaMCUnknown","Isolated Number of non #gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaIsoUnknown->SetYTitle("#eta");
-      fhEtaIsoUnknown->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhEtaIsoUnknown) ;
-    }//Histos with MC
-    
-  }
-  
-  if(fMakeSeveralIC){
-    const Int_t buffersize = 255;
-               char name[buffersize];
-               char title[buffersize];
-               for(Int_t icone = 0; icone<fNCones; 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()){
-                   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]) ; 
-                   
-                   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]) ; 
-                   
-                   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]) ; 
-                   
-                   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]) ; 
-                   
-                   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]) ; 
-                   
-                   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)");
-                   outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ; 
-                   
-                 }//Histos with MC
-                 
-                 for(Int_t ipt = 0; ipt<fNPtThresFrac;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]) ; 
-                   
-                   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()){
-                     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]) ; 
-                     
-                     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]) ; 
-                     
-                     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]) ; 
-                     
-                     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]) ; 
-                     
-                     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]) ; 
-                     
-                     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]) ; 
-                     
-                     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]) ; 
-                     
-                     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]) ;
-                     
-                     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]) ; 
-                     
-                     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]) ;
-                     
-                     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]) ; 
-                     
-                     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]) ;  
-                     
-                   }//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)) ;
-       delete nmsHistos;
-  }
-  
-  return outputContainer ;
-  
-}
-
-//__________________________________________________________________
-void  AliAnaParticleIsolation::MakeAnalysisFillAOD() 
-{
-  //Do analysis and fill aods
-  //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
-  
-  if(!GetInputAODBranch()){
-    printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
-    abort();
-  }
-  
-  if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
-       printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
-       abort();
-  }
-       
-  Int_t n = 0, nfrac = 0;
-  Bool_t isolated = kFALSE ; 
-  Float_t coneptsum = 0 ;
-  TObjArray * pl = 0x0; ; 
-  
-  //Select the calorimeter for candidate isolation with neutral particles
-  if(fCalorimeter == "PHOS")
-    pl = GetPHOSClusters();
-  else if (fCalorimeter == "EMCAL")
-    pl = GetEMCALClusters();
-  
-  //Loop on AOD branch, filled previously in AliAnaPhoton
-  TLorentzVector mom ;
-  Int_t naod = GetInputAODBranch()->GetEntriesFast();
-  if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
-  for(Int_t iaod = 0; iaod < naod; iaod++){
-    AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
-    
-    //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);
-    if(decay) continue ;
-
-    //After cuts, study isolation
-    n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
-    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);  
-         
-  }//loop
-  
-  if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");  
-  
-}
-
-//__________________________________________________________________
-void  AliAnaParticleIsolation::MakeAnalysisFillHistograms() 
-{
-  //Do analysis and fill histograms
-  Int_t n = 0, nfrac = 0;
-  Bool_t isolated = kFALSE ; 
-  Float_t coneptsum = 0 ;
-  //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();              
-    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");
-    //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);
-      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(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);    
-    }
-    
-    //Fill pt distribution of particles in cone
-    //Tracks
-    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++){
-        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(!fReMakeIC) fhConeSumPt->Fill(ptcluster,coneptsum);
-    
-    if(isolation){    
-      
-      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);
-        }
-      }//Histograms with MC
-      
-    }//Isolated histograms
-    
-  }// aod loop
-  
-}
-
-//____________________________________________________________________________
-void AliAnaParticleIsolation::InitParameters()
-{
-  
-  //Initialize the parameters of the analysis.
-  SetInputAODName("PWG4Particle");
-  SetAODObjArrayName("IsolationCone");  
-  AddToHistogramsName("AnaIsolation_");
-
-  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->GetObjArray(GetAODObjArrayName()+"Tracks"), 
-                                 ph->GetObjArray(GetAODObjArrayName()+"Clusters"),
-                                 GetReader(), 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(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(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(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;
-         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(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(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
-      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
-      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
-      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
-      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(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() ) ;
-  AliAnaPartCorrBaseClass::Print(" ");
-  
-  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") ;
-  
-} 
-
+/**************************************************************************\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
+//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)\r
+//////////////////////////////////////////////////////////////////////////////\r
+  \r
+  \r
+// --- ROOT system --- \r
+#include <TClonesArray.h>\r
+#include <TList.h>\r
+#include <TObjString.h>\r
+#include <TH2F.h>\r
+//#include <Riostream.h>\r
+#include <TClass.h>\r
+\r
+// --- Analysis system --- \r
+#include "AliAnaParticleIsolation.h" \r
+#include "AliCaloTrackReader.h"\r
+#include "AliIsolationCut.h"\r
+#include "AliNeutralMesonSelection.h"\r
+#include "AliAODPWG4ParticleCorrelation.h"\r
+#include "AliMCAnalysisUtils.h"\r
+#include "AliVTrack.h"\r
+#include "AliVCluster.h"\r
+\r
+ClassImp(AliAnaParticleIsolation)\r
+  \r
+//____________________________________________________________________________\r
+  AliAnaParticleIsolation::AliAnaParticleIsolation() : \r
+    AliAnaPartCorrBaseClass(), fCalorimeter(""), \r
+    fReMakeIC(0), fMakeSeveralIC(0), fMakeInvMass(0),\r
+    fhPtIso(0),fhPhiIso(0),fhEtaIso(0), fhPtNoIso(0),fhPtInvMassDecayIso(0),fhPtInvMassDecayNoIso(0), fhConeSumPt(0), fhPtInCone(0),\r
+    fhFRConeSumPt(0), fhPtInFRCone(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
+    fhPtNoIsoPi0Decay(0),fhPtNoIsoPrompt(0),fhPtIsoMCPhoton(0),fhPtNoIsoMCPhoton(0),\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
+  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() \r
+{\r
+  //dtor\r
+  //do not delete histograms\r
+  \r
+  //delete [] fConeSizes ; \r
+  //delete [] fPtThresholds ; \r
+  //delete [] fPtFractions ; \r
+\r
+}\r
+\r
+//_________________________________________________________________________\r
+Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4Particle * part1)\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("AliAnaParticleIsolation::CheckInvMass() - Selected gamma pair: pt %f, phi %f, eta%f\n",(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
+TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()\r
+{ \r
+       //Save parameters used for analysis\r
+        TString parList ; //this will be list of parameters used for this analysis.\r
+   const Int_t buffersize = 255;\r
+        char onePar[buffersize] ;\r
+        \r
+        snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;\r
+        parList+=onePar ;      \r
+        snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;\r
+        parList+=onePar ;\r
+        snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;\r
+        parList+=onePar ;\r
+        snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;\r
+        parList+=onePar ;\r
+        snprintf(onePar, buffersize,"fMakeInvMass=%d (Flag for rejection of candidates with a pi0 inv mass pair) \n",fMakeInvMass) ;\r
+        parList+=onePar ;\r
+        \r
+        if(fMakeSeveralIC){\r
+          snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;\r
+          parList+=onePar ;\r
+          snprintf(onePar, buffersize,"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
+            snprintf(onePar, buffersize,"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
+            snprintf(onePar, buffersize,"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
+            snprintf(onePar, buffersize,"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
+        return new TObjString(parList) ;\r
+       \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  = GetHistoPtBins();\r
+  Int_t nphibins = GetHistoPhiBins();\r
+  Int_t netabins = GetHistoEtaBins();\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 isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
+    fhConeSumPt->SetYTitle("#Sigma p_{T}");\r
+    fhConeSumPt->SetXTitle("p_{T} (GeV/c)");\r
+    outputContainer->Add(fhConeSumPt) ;\r
+    \r
+    fhPtInCone  = new TH2F\r
+      ("hPtInCone","p_{T} in isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);\r
+    fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");\r
+    fhPtInCone->SetXTitle("p_{T} (GeV/c)");\r
+    outputContainer->Add(fhPtInCone) ;\r
+    \r
+    fhFRConeSumPt  = new TH2F\r
+    ("hFRConePtSum","#Sigma p_{T} in the froward region isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);\r
+    fhFRConeSumPt->SetYTitle("#Sigma p_{T}");\r
+    fhFRConeSumPt->SetXTitle("p_{T} (GeV/c)");\r
+    outputContainer->Add(fhFRConeSumPt) ;\r
+    \r
+    fhPtInFRCone  = new TH2F\r
+    ("hPtInFRCone","p_{T} in froward region isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);\r
+    fhPtInFRCone->SetYTitle("p_{T in cone} (GeV/c)");\r
+    fhPtInFRCone->SetXTitle("p_{T} (GeV/c)");\r
+    outputContainer->Add(fhPtInFRCone) ;    \r
+    \r
+    fhPtIso  = new TH1F("hPt","Isolated Number of particles",nptbins,ptmin,ptmax); \r
+    fhPtIso->SetYTitle("N");\r
+    fhPtIso->SetXTitle("p_{T}(GeV/c)");\r
+    outputContainer->Add(fhPtIso) ; \r
+    \r
+    fhPhiIso  = new TH2F\r
+      ("hPhi","Isolated Number of particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
+    fhPhiIso->SetYTitle("#phi");\r
+    fhPhiIso->SetXTitle("p_{T} (GeV/c)");\r
+    outputContainer->Add(fhPhiIso) ; \r
+    \r
+    fhEtaIso  = new TH2F\r
+      ("hEta","Isolated Number of particles",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
+    fhEtaIso->SetYTitle("#eta");\r
+    fhEtaIso->SetXTitle("p_{T} (GeV/c)");\r
+    outputContainer->Add(fhEtaIso) ;\r
+\r
+    fhPtNoIso  = new TH1F("hPtNoIso","Number of not isolated leading particles",nptbins,ptmin,ptmax); \r
+    fhPtNoIso->SetYTitle("N");\r
+    fhPtNoIso->SetXTitle("p_{T}(GeV/c)");\r
+    outputContainer->Add(fhPtNoIso) ;\r
+\r
+    fhPtInvMassDecayIso  = new TH1F("hPtInvMassDecayIso","Number of isolated pi0 decay particles",nptbins,ptmin,ptmax); \r
+    fhPtNoIso->SetYTitle("N");\r
+    fhPtNoIso->SetXTitle("p_{T}(GeV/c)");\r
+    outputContainer->Add(fhPtInvMassDecayIso) ;\r
+\r
+    fhPtInvMassDecayNoIso  = new TH1F("hPtInvMassDecayNoIso","Number of not isolated leading pi0 decay particles",nptbins,ptmin,ptmax); \r
+    fhPtNoIso->SetYTitle("N");\r
+    fhPtNoIso->SetXTitle("p_{T}(GeV/c)");\r
+    outputContainer->Add(fhPtInvMassDecayNoIso) ;\r
+    \r
+    if(IsDataMC()){\r
+      \r
+      fhPtIsoPrompt  = new TH1F("hPtMCPrompt","Isolated Number of #gamma prompt",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
+       ("hPhiMCPrompt","Isolated Number of #gamma prompt",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
+       ("hEtaMCPrompt","Isolated Number of #gamma prompt ",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("hPtMCFragmentation","Isolated Number of #gamma",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
+       ("hPhiMCFragmentation","Isolated Number of #gamma fragmentation",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
+       ("hEtaMCFragmentation","Isolated Number of #gamma fragmentation",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("hPtMCPi0Decay","Isolated Number of #gamma from #pi^{0} decay",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
+       ("hPhiMCPi0Decay","Isolated Number of #gamma from #pi^{0} decay",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
+       ("hEtaMCPi0Decay","Isolated Number of #gamma from #pi^{0} decay",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("hPtMCOtherDecay","Isolated Number of #gamma from non #pi^{0} decay",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
+       ("hPhiMCOtherDecay","Isolated Number of #gamma from non #pi^{0} decay",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
+       ("hEtaMCOtherDecay","Isolated Number of #gamma non #pi^{0} decay",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("hPtMCConversion","Isolated Number of #gamma converted",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
+       ("hPhiMCConversion","Isolated Number of #gamma converted",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
+       ("hEtaMCConversion","Isolated Number of #gamma converted",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("hPtMCUnknown","Isolated Number of non #gamma particles",nptbins,ptmin,ptmax); \r
+      fhPtIsoUnknown->SetYTitle("N");\r
+      fhPtIsoUnknown->SetXTitle("p_{T}(GeV/c)");\r
+      outputContainer->Add(fhPtIsoUnknown) ; \r
+      \r
+      fhPhiIsoUnknown  = new TH2F\r
+       ("hPhiMCUnknown","Isolated Number of non #gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
+      fhPhiIsoUnknown->SetYTitle("#phi");\r
+      fhPhiIsoUnknown->SetXTitle("p_{T} (GeV/c)");\r
+      outputContainer->Add(fhPhiIsoUnknown) ; \r
+      \r
+      fhEtaIsoUnknown  = new TH2F\r
+       ("hEtaMCUnknown","Isolated Number of non #gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
+      fhEtaIsoUnknown->SetYTitle("#eta");\r
+      fhEtaIsoUnknown->SetXTitle("p_{T} (GeV/c)");\r
+      outputContainer->Add(fhEtaIsoUnknown) ;\r
+\r
+      fhPtNoIsoPi0Decay  = new TH1F\r
+       ("hPtNoIsoPi0Decay","Number of not isolated leading #gamma from Pi0 decay",nptbins,ptmin,ptmax); \r
+      fhPtNoIsoPi0Decay->SetYTitle("N");\r
+      fhPtNoIsoPi0Decay->SetXTitle("p_{T} (GeV/c)");\r
+      outputContainer->Add(fhPtNoIsoPi0Decay) ;\r
+\r
+      fhPtNoIsoPrompt  = new TH1F\r
+       ("hPtNoIsoPrompt","Number of not isolated leading prompt #gamma",nptbins,ptmin,ptmax); \r
+      fhPtNoIsoPrompt->SetYTitle("N");\r
+      fhPtNoIsoPrompt->SetXTitle("p_{T} (GeV/c)");\r
+      outputContainer->Add(fhPtNoIsoPrompt) ;\r
+\r
+      fhPtIsoMCPhoton  = new TH1F\r
+       ("hPtIsoMCPhoton","Number of isolated leading  #gamma",nptbins,ptmin,ptmax); \r
+      fhPtIsoMCPhoton->SetYTitle("N");\r
+      fhPtIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");\r
+      outputContainer->Add(fhPtIsoMCPhoton) ;\r
+\r
+      fhPtNoIsoMCPhoton  = new TH1F\r
+       ("hPtNoIsoMCPhoton","Number of not isolated leading  #gamma",nptbins,ptmin,ptmax); \r
+      fhPtNoIsoMCPhoton->SetYTitle("N");\r
+      fhPtNoIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");\r
+      outputContainer->Add(fhPtNoIsoMCPhoton) ;\r
+    }//Histos with MC\r
+    \r
+  }\r
+  \r
+  if(fMakeSeveralIC){\r
+    const Int_t buffersize = 255;\r
+               char name[buffersize];\r
+               char title[buffersize];\r
+               for(Int_t icone = 0; icone<fNCones; icone++){\r
+                 snprintf(name, buffersize,"hPtSum_Cone_%d",icone);\r
+                 snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtSumUnknown_Cone_%d",icone);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);\r
+                   snprintf(title, buffersize,"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
+                   snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);\r
+                   snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtThresMCUnknown_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+                     snprintf(name, buffersize,"hPtFracMCUnknown_Cone_%d_Pt%d",icone,ipt);\r
+                     snprintf(title, buffersize,"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
+       delete nmsHistos;\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
+    printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());\r
+    abort();\r
+  }\r
+  \r
+  if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){\r
+       printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());\r
+       abort();\r
+  }\r
+       \r
+  Int_t n = 0, nfrac = 0;\r
+  Bool_t isolated = kFALSE ;\r
+  Bool_t leading = kTRUE ;\r
+  Bool_t decay = kFALSE;\r
+  Float_t coneptsum = 0 ;\r
+  TObjArray * pl = 0x0; ; \r
+  \r
+  //Select the calorimeter for candidate isolation with neutral particles\r
+  if(fCalorimeter == "PHOS")\r
+    pl = GetPHOSClusters();\r
+  else if (fCalorimeter == "EMCAL")\r
+    pl = GetEMCALClusters();\r
+  \r
+  //Loop on AOD branch, filled previously in AliAnaPhoton\r
+  TLorentzVector mom ;\r
+  Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
+  if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);\r
+  for(Int_t iaod = 0; iaod < naod; iaod++){\r
+    AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
+    \r
+    //If too small or too large pt, skip\r
+    if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ; \r
+    \r
+    //check if it is low pt trigger particle, then adjust the isolation method\r
+    if(aodinput->Pt() < GetIsolationCut()->GetPtThreshold() || aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold()) \r
+      continue ; //trigger should not from underlying event\r
+ //   if(GetIsolationCut()->GetICMethod()!=AliIsolationCut::kPtThresIC && aodinput->Pt()!=0.){ \r
+//      //low pt trigger use pt threshold IC instead of other method\r
+//      if (aodinput->Pt()*(GetIsolationCut()->GetPtFraction())<GetIsolationCut()->GetPtThreshold() || aodinput->Pt()*(GetIsolationCut()->GetPtFraction())<GetIsolationCut()->GetSumPtThreshold()) {\r
+//       // printf("change the IC method to PtThresIC\n") ;\r
+//        GetIsolationCut()->SetICMethod(AliIsolationCut::kPtThresIC); \r
+//      }\r
+//    }\r
+    \r
+    //Check invariant mass, if pi0, tag decay.\r
+    if(fMakeInvMass) decay = CheckInvMass(iaod,aodinput);\r
+\r
+    //After cuts, study isolation\r
+    n=0; nfrac = 0; isolated = kFALSE; leading = kTRUE; coneptsum = 0;\r
+    GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,GetReader(), kTRUE, aodinput, GetAODObjArrayName(), n,nfrac,coneptsum, isolated, leading);\r
+    aodinput->SetIsolated(isolated);\r
+    aodinput->SetLeadingParticle(leading);\r
+    aodinput->SetTagged(decay);\r
+    if(GetDebug() > 1 && isolated) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",iaod);  \r
+         \r
+  }//loop\r
+  \r
+  if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");  \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
+  Bool_t leading = kTRUE ;\r
+  Float_t coneptsum = 0 ;\r
+  //Loop on stored AOD \r
+  Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
+  if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);\r
+  \r
+  //Get vertex for photon momentum calculation\r
+  Double_t vertex[]={0,0,0} ; //vertex ;\r
+  //Double_t vertex2[]={0,0,0} ; //vertex ;\r
+  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
+         GetReader()->GetVertex(vertex);\r
+         //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);\r
+  }    \r
+       \r
+  for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+    AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
+    \r
+    Bool_t  isolation  = aod->IsIsolated(); \r
+    Bool_t lead        = aod->IsLeadingParticle();\r
+    Bool_t dec         = aod->IsTagged();\r
+    Float_t ptcluster  = aod->Pt();\r
+    Float_t phicluster = aod->Phi();\r
+    Float_t etacluster = aod->Eta();\r
+   // Float_t phiForward = aod->Phi()+TMath::PiOver2() ;\r
+    Float_t conesize   = GetIsolationCut()->GetConeSize();\r
+    \r
+    //Recover reference arrays with clusters and tracks\r
+    TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");\r
+    TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");\r
+    //If too small or too large pt, skip\r
+    if(aod->Pt() < GetMinPt() || aod->Pt() > 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; leading = kTRUE; coneptsum = 0;\r
+      GetIsolationCut()->MakeIsolationCut(reftracks, refclusters, GetReader(), kFALSE, aod, "", n,nfrac,coneptsum, isolated, leading);\r
+      fhConeSumPt->Fill(ptcluster,coneptsum);    \r
+      if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);    \r
+    }\r
+    \r
+    //Fill pt distribution of particles in cone\r
+    //Tracks\r
+    coneptsum = 0;\r
+    Double_t sumptFR = 0. ;\r
+    TObjArray * trackList   = GetCTSTracks() ;\r
+    for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++){\r
+      AliVTrack* track = (AliVTrack *) trackList->At(itrack);\r
+      //fill the histograms at forward range\r
+      if(!track){\r
+        printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");\r
+        continue;\r
+      }\r
+      Double_t dPhi = phicluster - track->Phi() + TMath::PiOver2();\r
+      Double_t dEta = etacluster - track->Eta();\r
+      Double_t arg  = dPhi*dPhi + dEta*dEta;\r
+      if(TMath::Sqrt(arg) < conesize){\r
+        fhPtInFRCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));\r
+        sumptFR+=track->Pt();\r
+      }    \r
+      dPhi = phicluster - track->Phi() - TMath::PiOver2();\r
+      arg  = dPhi*dPhi + dEta*dEta;\r
+      if(TMath::Sqrt(arg) < conesize){\r
+        fhPtInFRCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));\r
+        sumptFR+=track->Pt();\r
+      }      \r
+    }\r
+    fhFRConeSumPt->Fill(ptcluster,sumptFR);\r
+    if(reftracks){  \r
+      for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++){\r
+        AliVTrack* track = (AliVTrack *) reftracks->At(itrack);\r
+        fhPtInCone->Fill(ptcluster,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));\r
+        coneptsum+=track->Pt();\r
+      }\r
+    }\r
+\r
+    //CaloClusters\r
+    if(refclusters){    \r
+      TLorentzVector mom ;\r
+      for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++){\r
+        AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);\r
+        Int_t input = 0;\r
+        //                     if     (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo) input = 1 ;\r
+        //                     else if(fCalorimeter == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= icalo) input = 1;\r
+        \r
+        //Get Momentum vector, \r
+        if     (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line\r
+        //                     else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  \r
+        \r
+        fhPtInCone->Fill(ptcluster, mom.Pt());\r
+        coneptsum+=mom.Pt();\r
+      }\r
+    }\r
+         \r
+    if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);\r
+    \r
+    if(!fReMakeIC) fhConeSumPt->Fill(ptcluster,coneptsum);\r
+    \r
+    if(isolation && lead){    \r
+      \r
+      if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED, fill histograms\n", iaod);\r
+      \r
+      fhPtIso  ->Fill(ptcluster);\r
+      fhPhiIso ->Fill(ptcluster,phicluster);\r
+      fhEtaIso ->Fill(ptcluster,etacluster);\r
+      if (dec) fhPtInvMassDecayIso->Fill(ptcluster);\r
+      \r
+      if(IsDataMC()){\r
+        Int_t tag =aod->GetTag();\r
+        \r
+        if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){\r
+          fhPtIsoPrompt  ->Fill(ptcluster);\r
+          fhPhiIsoPrompt ->Fill(ptcluster,phicluster);\r
+          fhEtaIsoPrompt ->Fill(ptcluster,etacluster);\r
+        }\r
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))\r
+        {\r
+          fhPtIsoFragmentation  ->Fill(ptcluster);\r
+          fhPhiIsoFragmentation ->Fill(ptcluster,phicluster);\r
+          fhEtaIsoFragmentation ->Fill(ptcluster,etacluster);\r
+        }\r
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))\r
+        {\r
+          fhPtIsoPi0Decay  ->Fill(ptcluster);\r
+          fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster);\r
+          fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster);\r
+        }\r
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))\r
+        {\r
+          fhPtIsoOtherDecay  ->Fill(ptcluster);\r
+          fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster);\r
+          fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster);\r
+        }\r
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))\r
+        {\r
+          fhPtIsoConversion  ->Fill(ptcluster);\r
+          fhPhiIsoConversion ->Fill(ptcluster,phicluster);\r
+          fhEtaIsoConversion ->Fill(ptcluster,etacluster);\r
+        }\r
+       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))\r
+        {\r
+          fhPtIsoMCPhoton  ->Fill(ptcluster);\r
+        }\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
+    if(!isolation && lead)\r
+      {\r
+       fhPtNoIso  ->Fill(ptcluster);\r
+       if (dec) fhPtInvMassDecayNoIso->Fill(ptcluster);\r
+\r
+       if(IsDataMC()){\r
+        Int_t tag =aod->GetTag();\r
+       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))\r
+         {\r
+            fhPtNoIsoPi0Decay->Fill(ptcluster);\r
+         }\r
+       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))\r
+         {\r
+           fhPtNoIsoPrompt->Fill(ptcluster);\r
+         }\r
+       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))\r
+         {\r
+           fhPtNoIsoMCPhoton->Fill(ptcluster);\r
+         }\r
+       \r
+       }\r
+      }\r
+    \r
+  }// aod loop\r
+  \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaParticleIsolation::InitParameters()\r
+{\r
+  \r
+  //Initialize the parameters of the analysis.\r
+  SetInputAODName("PWG4Particle");\r
+  SetAODObjArrayName("IsolationCone");  \r
+  AddToHistogramsName("AnaIsolation_");\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
+  Bool_t  leading   = kTRUE;\r
+\r
+  //Loop on cone sizes\r
+  for(Int_t icone = 0; icone<fNCones; icone++){\r
+    GetIsolationCut()->SetConeSize(fConeSizes[icone]);\r
+    coneptsum = 0 ;\r
+\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(ph->GetObjArray(GetAODObjArrayName()+"Tracks"), \r
+                                 ph->GetObjArray(GetAODObjArrayName()+"Clusters"),\r
+                                         GetReader(), kFALSE, ph, "",n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated, leading);\r
+\r
+      //Normal ptThreshold cut\r
+      if(n[icone][ipt] == 0) {\r
+       fhPtThresIsolated[icone][ipt]->Fill(ptC);\r
+       if(IsDataMC()){\r
+         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(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(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
+         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(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(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;\r
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;\r
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;\r
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;\r
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(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
+  AliAnaPartCorrBaseClass::Print(" ");\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
index ff3c3369023f3095f9dab642f7594368a274aeeb..911a8d4b45bb5ebbd8d79c5eb24d477f99798162 100755 (executable)
-#ifndef ALIANAPARTICLEISOLATION_H
-#define ALIANAPARTICLEISOLATION_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice     */
-/* $Id: AliAnaParticleIsolation.h 27413 2008-07-18 13:28:12Z gconesab $ */
-
-//_________________________________________________________________________
-
-// Class for the 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 (INFN-LNF)
-
-// --- ROOT system ---
-class TH2F;
-class TList ;
-class TObjString;
-
-// --- ANALYSIS system ---
-#include "AliAnaPartCorrBaseClass.h"
-class AliAODPWG4Particle;
-class AliAODPWG4ParticleCorrelation ;
-
-
-class AliAnaParticleIsolation : public AliAnaPartCorrBaseClass {
-
- public:   
-  AliAnaParticleIsolation() ; // default ctor
-  virtual ~AliAnaParticleIsolation() ; //virtual dtor
-
- private:
-  AliAnaParticleIsolation(const AliAnaParticleIsolation & g) ; // cpy ctor
-  AliAnaParticleIsolation & operator = (const AliAnaParticleIsolation & g) ;//cpy assignment
-
- public:
-
-  Bool_t CheckInvMass(const Int_t icalo,const AliAODPWG4Particle * ph) ;
-  
-  TObjString * GetAnalysisCuts();
-  TList      * GetCreateOutputObjects();
-  
-  void MakeAnalysisFillAOD()  ;
-  
-  void MakeAnalysisFillHistograms() ; 
-  
-  void MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation * ph); 
-  
-  void Print(const Option_t * opt)const;
-  
-  TString GetCalorimeter()   const {return fCalorimeter ; }
-  void SetCalorimeter(TString & det)    {fCalorimeter = det ; }
-  
-  Int_t   GetNCones()              const {return fNCones ; }
-  Int_t   GetNPtThresFrac()        const {return fNPtThresFrac ; }
-  Float_t GetConeSizes(Int_t i)    const {return fConeSizes[i] ; }
-  Float_t GetPtThresholds(Int_t i) const {return fPtThresholds[i] ; }
-  Float_t GetPtFractions(Int_t i)  const {return fPtFractions[i] ; }
-  
-  void InitParameters();
-  
-  void SetNCones(Int_t ncs)                  {fNCones = ncs ; }
-  void SetNPtThresFrac(Int_t npt)            {fNPtThresFrac = npt; }
-  void SetConeSizes(Int_t i, Float_t r)      {fConeSizes[i] = r ; }
-  void SetPtThresholds(Int_t i, Float_t pt)  {fPtThresholds[i] = pt; }
-  void SetPtFractions(Int_t i, Float_t pt)   {fPtFractions[i] = pt; } 
-  
-  Bool_t IsReIsolationOn() const {return fReMakeIC ; }
-  void SwitchOnReIsolation()  { fReMakeIC = kTRUE;}
-  void SwitchOffReIsolation() { fReMakeIC = kFALSE;}
-  
-  Bool_t IsSeveralIsolationOn() const {return fMakeSeveralIC ; }
-  void SwitchOnSeveralIsolation()  { fMakeSeveralIC = kTRUE;}
-  void SwitchOffSeveralIsolation() { fMakeSeveralIC = kFALSE;}
-  
-  Bool_t IsInvariantMassOn() const {return fMakeInvMass ; }
-  void SwitchOnInvariantMass()  { fMakeInvMass = kTRUE;}
-  void SwitchOffInvariantMass() { fMakeInvMass = kFALSE;}
-  
-  //Histogrammes setters and getters
-  virtual void SetHistoPtSumRangeAndNBins(Float_t min, Float_t max, Int_t n) {
-    fHistoNPtSumBins = n ;
-    fHistoPtSumMax = max ;
-    fHistoPtSumMin = min ;
-  }
-  
-  Int_t   GetHistoNPtSumBins() const { return fHistoNPtSumBins ; }
-  Float_t GetHistoPtSumMin()   const { return fHistoPtSumMin ; }
-  Float_t GetHistoPtSumMax()   const { return fHistoPtSumMax ; }
-  
-  virtual void SetHistoPtInConeRangeAndNBins(Float_t min, Float_t max, Int_t n) {
-    fHistoNPtInConeBins = n ;
-    fHistoPtInConeMax = max ;
-    fHistoPtInConeMin = min ;
-  }
-  
-  Int_t   GetHistoNPtInConeBins() const { return fHistoNPtInConeBins ; }
-  Float_t GetHistoPtInConeMin()   const { return fHistoPtInConeMin ; }
-  Float_t GetHistoPtInConeMax()   const { return fHistoPtInConeMax ; }
-  
-  
- private:
-  
-  TString  fCalorimeter ;   // Calorimeter where neutral particles in cone for isolation are;
-  Bool_t   fReMakeIC ;      // Do isolation analysis
-  Bool_t   fMakeSeveralIC ; // Do analysis for different IC
-  Bool_t   fMakeInvMass;    // Select candidate if no pair from decay
-  
-  //Histograms  
-  
-  TH1F * fhPtIso ;     //! Number of isolated particles
-  TH2F * fhPhiIso ;    //! Phi of isolated particles
-  TH2F * fhEtaIso ;    //! eta of isolated particles
-  TH2F * fhConeSumPt ; //! Sum Pt in the cone
-  TH2F * fhPtInCone ;  //! Particle Pt in the cone
-  TH2F * fhFRConeSumPt ; //! Sum Pt in the forward region cone (phi +90)
-  TH2F * fhPtInFRCone ;  //! Particle Pt in the forward region cone (phi +90 ) 
-  
-  //Prompt photon analysis data members for multiple cones and pt thresholds 
-  Int_t       fNCones ;          //! Number of cone sizes to test
-  Int_t       fNPtThresFrac ;    //! Number of ptThres and ptFrac to test
-  
-  Float_t   fConeSizes[5] ;    //! Array with cones to test
-  Float_t   fPtThresholds[5] ; //! Array with pt thresholds to test
-  Float_t   fPtFractions[5] ;  //! Array with pt thresholds to test
-  
-  TH1F* fhPtThresIsolated[5][5] ; //! Isolated particle with pt threshold 
-  TH1F* fhPtFracIsolated[5][5] ;  //! Isolated particle with pt threshold 
-  TH2F* fhPtSumIsolated[5] ;      //! Isolated particle with threshold on cone pt sum
-  
-  //MC
-  TH1F * fhPtIsoPrompt;   //! Number of isolated prompt gamma 
-  TH2F * fhPhiIsoPrompt;  //! Phi of isolated prompt gamma
-  TH2F * fhEtaIsoPrompt;  //! eta of isolated prompt gamma
-  TH1F * fhPtThresIsolatedPrompt[5][5];   //! Isolated prompt gamma with pt threshold 
-  TH1F * fhPtFracIsolatedPrompt[5][5];    //! Isolated prompt gamma with pt frac
-  TH2F * fhPtSumIsolatedPrompt[5];        //! Isolated prompt gamma with threshold on cone pt sume
-  TH1F * fhPtIsoFragmentation;   //! Number of isolated fragmentation gamma 
-  TH2F * fhPhiIsoFragmentation;  //! Phi of isolated fragmentation gamma
-  TH2F * fhEtaIsoFragmentation;  //! eta of isolated fragmentation gamma
-  TH1F * fhPtThresIsolatedFragmentation[5][5];  //! Isolated fragmentation gamma with pt threshold 
-  TH1F * fhPtFracIsolatedFragmentation[5][5];   //! Isolated fragmentation gamma with pt frac
-  TH2F * fhPtSumIsolatedFragmentation[5];       //! Isolated fragmentation gamma with threshold on cone pt sume
-  TH1F * fhPtIsoPi0Decay;   //!Number of isolated Pi0Decay gamma 
-  TH2F * fhPhiIsoPi0Decay;  //! Phi of isolated Pi0Decay gamma
-  TH2F * fhEtaIsoPi0Decay;  //! eta of isolated Pi0Decay gamma
-  TH1F * fhPtThresIsolatedPi0Decay[5][5];  //! Isolated Pi0Decay gamma with pt threshold 
-  TH1F * fhPtFracIsolatedPi0Decay[5][5];   //! Isolated Pi0Decay gamma with pt frac
-  TH2F * fhPtSumIsolatedPi0Decay[5];       //! Isolated Pi0Decay gamma with threshold on cone pt sume
-  TH1F * fhPtIsoOtherDecay;   //! Number of isolated OtherDecay gamma 
-  TH2F * fhPhiIsoOtherDecay;  //! Phi of isolated OtherDecay gamma
-  TH2F * fhEtaIsoOtherDecay;  //! eta of isolated OtherDecay gamma
-  TH1F * fhPtThresIsolatedOtherDecay[5][5];  //! Isolated OtherDecay gamma with pt threshold 
-  TH1F * fhPtFracIsolatedOtherDecay[5][5];   //! Isolated OtherDecay gamma with pt frac
-  TH2F * fhPtSumIsolatedOtherDecay[5];       //! Isolated OtherDecay gamma with threshold on cone pt sume      
-  TH1F * fhPtIsoConversion;   //! Number of isolated Conversion gamma 
-  TH2F * fhPhiIsoConversion;  //! Phi of isolated Conversion gamma
-  TH2F * fhEtaIsoConversion;  //! eta of isolated Conversion gamma
-  TH1F * fhPtThresIsolatedConversion[5][5];  //! Isolated Conversion gamma with pt threshold 
-  TH1F * fhPtFracIsolatedConversion[5][5];   //! Isolated Conversion gamma with pt frac
-  TH2F * fhPtSumIsolatedConversion[5];       //! Isolated Conversion gamma with threshold on cone pt sume
-  TH1F * fhPtIsoUnknown;   //! Number of isolated Unknown 
-  TH2F * fhPhiIsoUnknown;  //! Phi of isolated Unknown
-  TH2F * fhEtaIsoUnknown;  //! eta of isolated Unknown
-  TH1F * fhPtThresIsolatedUnknown[5][5];  //! Isolated Unknown gamma with pt threshold 
-  TH1F * fhPtFracIsolatedUnknown[5][5];   //! Isolated Unknown gamma with pt frac
-  TH2F * fhPtSumIsolatedUnknown[5];       //! Isolated Unknown gamma with threshold on cone pt sume
-  
-  //Histograms settings
-  Int_t   fHistoNPtSumBins; // Number of bins in PtSum histograms
-  Float_t fHistoPtSumMax;   // PtSum maximum in histogram
-  Float_t fHistoPtSumMin;        // PtSum minimum in histogram
-  Int_t   fHistoNPtInConeBins; // Number of bins in PtInCone histogram
-  Float_t fHistoPtInConeMax;   // PtInCone maximum in histogram
-  Float_t fHistoPtInConeMin;   // PtInCone maximum in histogram 
-  
-  ClassDef(AliAnaParticleIsolation,3)
-    } ;
-
-
-#endif //ALIANAPARTICLEISOLATION_H
-
-
-
+#ifndef ALIANAPARTICLEISOLATION_H\r
+#define ALIANAPARTICLEISOLATION_H\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice     */\r
+/* $Id: AliAnaParticleIsolation.h 27413 2008-07-18 13:28:12Z gconesab $ */\r
+\r
+//_________________________________________________________________________\r
+\r
+// Class for the 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 (INFN-LNF)\r
+\r
+// --- ROOT system ---\r
+class TH2F;\r
+class TList ;\r
+class TObjString;\r
+\r
+// --- ANALYSIS system ---\r
+#include "AliAnaPartCorrBaseClass.h"\r
+class AliAODPWG4Particle;\r
+class AliAODPWG4ParticleCorrelation ;\r
+\r
+\r
+class AliAnaParticleIsolation : public AliAnaPartCorrBaseClass {\r
+\r
+ public:   \r
+  AliAnaParticleIsolation() ; // default ctor\r
+  virtual ~AliAnaParticleIsolation() ; //virtual dtor\r
+\r
+ private:\r
+  AliAnaParticleIsolation(const AliAnaParticleIsolation & g) ; // cpy ctor\r
+  AliAnaParticleIsolation & operator = (const AliAnaParticleIsolation & g) ;//cpy assignment\r
+\r
+ public:\r
+\r
+  Bool_t CheckInvMass(const Int_t icalo,const AliAODPWG4Particle * ph) ;\r
+  \r
+  TObjString * GetAnalysisCuts();\r
+  TList      * GetCreateOutputObjects();\r
+  \r
+  void MakeAnalysisFillAOD()  ;\r
+  \r
+  void MakeAnalysisFillHistograms() ; \r
+  \r
+  void MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation * ph); \r
+  \r
+  void Print(const Option_t * opt)const;\r
+  \r
+  TString GetCalorimeter()   const {return fCalorimeter ; }\r
+  void SetCalorimeter(TString & det)    {fCalorimeter = det ; }\r
+  \r
+  Int_t   GetNCones()              const {return fNCones ; }\r
+  Int_t   GetNPtThresFrac()        const {return fNPtThresFrac ; }\r
+  Float_t GetConeSizes(Int_t i)    const {return fConeSizes[i] ; }\r
+  Float_t GetPtThresholds(Int_t i) const {return fPtThresholds[i] ; }\r
+  Float_t GetPtFractions(Int_t i)  const {return fPtFractions[i] ; }\r
+  \r
+  void InitParameters();\r
+  \r
+  void SetNCones(Int_t ncs)                  {fNCones = ncs ; }\r
+  void SetNPtThresFrac(Int_t npt)            {fNPtThresFrac = npt; }\r
+  void SetConeSizes(Int_t i, Float_t r)      {fConeSizes[i] = r ; }\r
+  void SetPtThresholds(Int_t i, Float_t pt)  {fPtThresholds[i] = pt; }\r
+  void SetPtFractions(Int_t i, Float_t pt)   {fPtFractions[i] = pt; } \r
+  \r
+  Bool_t IsReIsolationOn() const {return fReMakeIC ; }\r
+  void SwitchOnReIsolation()  { fReMakeIC = kTRUE;}\r
+  void SwitchOffReIsolation() { fReMakeIC = kFALSE;}\r
+  \r
+  Bool_t IsSeveralIsolationOn() const {return fMakeSeveralIC ; }\r
+  void SwitchOnSeveralIsolation()  { fMakeSeveralIC = kTRUE;}\r
+  void SwitchOffSeveralIsolation() { fMakeSeveralIC = kFALSE;}\r
+  \r
+  Bool_t IsInvariantMassOn() const {return fMakeInvMass ; }\r
+  void SwitchOnInvariantMass()  { fMakeInvMass = kTRUE;}\r
+  void SwitchOffInvariantMass() { fMakeInvMass = kFALSE;}\r
+  \r
+  //Histogrammes setters and getters\r
+  virtual void SetHistoPtSumRangeAndNBins(Float_t min, Float_t max, Int_t n) {\r
+    fHistoNPtSumBins = n ;\r
+    fHistoPtSumMax = max ;\r
+    fHistoPtSumMin = min ;\r
+  }\r
+  \r
+  Int_t   GetHistoNPtSumBins() const { return fHistoNPtSumBins ; }\r
+  Float_t GetHistoPtSumMin()   const { return fHistoPtSumMin ; }\r
+  Float_t GetHistoPtSumMax()   const { return fHistoPtSumMax ; }\r
+  \r
+  virtual void SetHistoPtInConeRangeAndNBins(Float_t min, Float_t max, Int_t n) {\r
+    fHistoNPtInConeBins = n ;\r
+    fHistoPtInConeMax = max ;\r
+    fHistoPtInConeMin = min ;\r
+  }\r
+  \r
+  Int_t   GetHistoNPtInConeBins() const { return fHistoNPtInConeBins ; }\r
+  Float_t GetHistoPtInConeMin()   const { return fHistoPtInConeMin ; }\r
+  Float_t GetHistoPtInConeMax()   const { return fHistoPtInConeMax ; }\r
+  \r
+  \r
+ private:\r
+  \r
+  TString  fCalorimeter ;   // Calorimeter where neutral particles in cone for isolation are;\r
+  Bool_t   fReMakeIC ;      // Do isolation analysis\r
+  Bool_t   fMakeSeveralIC ; // Do analysis for different IC\r
+  Bool_t   fMakeInvMass;    // Select candidate if no pair from decay\r
+  \r
+  //Histograms  \r
+  \r
+  TH1F * fhPtIso ;     //! Number of isolated particles\r
+  TH2F * fhPhiIso ;    //! Phi of isolated particles\r
+  TH2F * fhEtaIso ;    //! eta of isolated particles\r
+  TH1F * fhPtNoIso ;   //! Number of not isolated leading particles\r
+  TH1F * fhPtInvMassDecayIso ;   //! Number of isolated Pi0 decay particles (invariant mass tag)\r
+  TH1F * fhPtInvMassDecayNoIso ;   //! Number of not isolated Pi0 decay leading particles (invariant mass tag)\r
+  TH2F * fhConeSumPt ; //! Sum Pt in the cone\r
+  TH2F * fhPtInCone ;  //! Particle Pt in the cone\r
+  TH2F * fhFRConeSumPt ; //! Sum Pt in the forward region cone (phi +90)\r
+  TH2F * fhPtInFRCone ;  //! Particle Pt in the forward region cone (phi +90 ) \r
+  \r
+  //Prompt photon analysis data members for multiple cones and pt thresholds \r
+  Int_t       fNCones ;          //! Number of cone sizes to test\r
+  Int_t       fNPtThresFrac ;    //! Number of ptThres and ptFrac to test\r
+  \r
+  Float_t   fConeSizes[5] ;    //! Array with cones to test\r
+  Float_t   fPtThresholds[5] ; //! Array with pt thresholds to test\r
+  Float_t   fPtFractions[5] ;  //! Array with pt thresholds to test\r
+  \r
+  TH1F* fhPtThresIsolated[5][5] ; //! Isolated particle with pt threshold \r
+  TH1F* fhPtFracIsolated[5][5] ;  //! Isolated particle with pt threshold \r
+  TH2F* fhPtSumIsolated[5] ;      //! Isolated particle with threshold on cone pt sum\r
+  \r
+  //MC\r
+  TH1F * fhPtIsoPrompt;   //! Number of isolated prompt gamma \r
+  TH2F * fhPhiIsoPrompt;  //! Phi of isolated prompt gamma\r
+  TH2F * fhEtaIsoPrompt;  //! eta of isolated prompt gamma\r
+  TH1F * fhPtThresIsolatedPrompt[5][5];   //! Isolated prompt gamma with pt threshold \r
+  TH1F * fhPtFracIsolatedPrompt[5][5];    //! Isolated prompt gamma with pt frac\r
+  TH2F * fhPtSumIsolatedPrompt[5];        //! Isolated prompt gamma with threshold on cone pt sume\r
+  TH1F * fhPtIsoFragmentation;   //! Number of isolated fragmentation gamma \r
+  TH2F * fhPhiIsoFragmentation;  //! Phi of isolated fragmentation gamma\r
+  TH2F * fhEtaIsoFragmentation;  //! eta of isolated fragmentation gamma\r
+  TH1F * fhPtThresIsolatedFragmentation[5][5];  //! Isolated fragmentation gamma with pt threshold \r
+  TH1F * fhPtFracIsolatedFragmentation[5][5];   //! Isolated fragmentation gamma with pt frac\r
+  TH2F * fhPtSumIsolatedFragmentation[5];       //! Isolated fragmentation gamma with threshold on cone pt sume\r
+  TH1F * fhPtIsoPi0Decay;   //!Number of isolated Pi0Decay gamma \r
+  TH2F * fhPhiIsoPi0Decay;  //! Phi of isolated Pi0Decay gamma\r
+  TH2F * fhEtaIsoPi0Decay;  //! eta of isolated Pi0Decay gamma\r
+  TH1F * fhPtThresIsolatedPi0Decay[5][5];  //! Isolated Pi0Decay gamma with pt threshold \r
+  TH1F * fhPtFracIsolatedPi0Decay[5][5];   //! Isolated Pi0Decay gamma with pt frac\r
+  TH2F * fhPtSumIsolatedPi0Decay[5];       //! Isolated Pi0Decay gamma with threshold on cone pt sume\r
+  TH1F * fhPtIsoOtherDecay;   //! Number of isolated OtherDecay gamma \r
+  TH2F * fhPhiIsoOtherDecay;  //! Phi of isolated OtherDecay gamma\r
+  TH2F * fhEtaIsoOtherDecay;  //! eta of isolated OtherDecay gamma\r
+  TH1F * fhPtThresIsolatedOtherDecay[5][5];  //! Isolated OtherDecay gamma with pt threshold \r
+  TH1F * fhPtFracIsolatedOtherDecay[5][5];   //! Isolated OtherDecay gamma with pt frac\r
+  TH2F * fhPtSumIsolatedOtherDecay[5];       //! Isolated OtherDecay gamma with threshold on cone pt sume      \r
+  TH1F * fhPtIsoConversion;   //! Number of isolated Conversion gamma \r
+  TH2F * fhPhiIsoConversion;  //! Phi of isolated Conversion gamma\r
+  TH2F * fhEtaIsoConversion;  //! eta of isolated Conversion gamma\r
+  TH1F * fhPtThresIsolatedConversion[5][5];  //! Isolated Conversion gamma with pt threshold \r
+  TH1F * fhPtFracIsolatedConversion[5][5];   //! Isolated Conversion gamma with pt frac\r
+  TH2F * fhPtSumIsolatedConversion[5];       //! Isolated Conversion gamma with threshold on cone pt sume\r
+  TH1F * fhPtIsoUnknown;   //! Number of isolated Unknown \r
+  TH2F * fhPhiIsoUnknown;  //! Phi of isolated Unknown\r
+  TH2F * fhEtaIsoUnknown;  //! eta of isolated Unknown\r
+  TH1F * fhPtThresIsolatedUnknown[5][5];  //! Isolated Unknown gamma with pt threshold \r
+  TH1F * fhPtFracIsolatedUnknown[5][5];   //! Isolated Unknown gamma with pt frac\r
+  TH2F * fhPtSumIsolatedUnknown[5];       //! Isolated Unknown gamma with threshold on cone pt sume\r
+\r
+  TH1F * fhPtNoIsoPi0Decay;  //! Number of not isolated leading Pi0Decay gamma \r
+  TH1F * fhPtNoIsoPrompt;   //! Number of not isolated leading prompt gamma \r
+  TH1F * fhPtIsoMCPhoton;  //! Number of isolated leading gamma \r
+  TH1F * fhPtNoIsoMCPhoton;   //! Number of not isolated leading gamma \r
+\r
+  //Histograms settings\r
+  Int_t   fHistoNPtSumBins; // Number of bins in PtSum histograms\r
+  Float_t fHistoPtSumMax;   // PtSum maximum in histogram\r
+  Float_t fHistoPtSumMin;        // PtSum minimum in histogram\r
+  Int_t   fHistoNPtInConeBins; // Number of bins in PtInCone histogram\r
+  Float_t fHistoPtInConeMax;   // PtInCone maximum in histogram\r
+  Float_t fHistoPtInConeMin;   // PtInCone maximum in histogram \r
+  \r
+  ClassDef(AliAnaParticleIsolation,3)\r
+    } ;\r
+\r
+\r
+#endif //ALIANAPARTICLEISOLATION_H\r
+\r
+\r
+\r