[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliIsolationCut.cxx

/**************************************************************************
 * 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.                  *
 **************************************************************************/

//_________________________________________________________________________
// 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 <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"
#include "AliCaloPID.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     = kOnlyCharged;
  fICMethod       = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
  
}

//________________________________________________________________________________
void  AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS, 
                                        const TObjArray * plNe, 
                                        const AliCaloTrackReader * reader, 
                                        const AliCaloPID * pid, 
                                        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() ;
  if(phiC<0) phiC+=TMath::TwoPi();
  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) || 
         track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)   ) 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 = Radius(etaC, phiC, eta, phi);
      
      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))
  {
    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( pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
    
      //Assume that come from vertex in straight line
      calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
      
      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 = Radius(etaC, phiC, eta, phi);
      
      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 selected isolation criteria
  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(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
    if(fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) 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") ;
  
} 

//___________________________________________________________________________
Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC, 
                                const Float_t eta , const Float_t phi) const
{
  // Calculate the distance to trigger from any particle

  Float_t dEta = etaC-eta;
  Float_t dPhi = phiC-phi;
  
  if(TMath::Abs(dPhi) >= TMath::Pi()) 
    dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
  
  return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
  
}
Commit	Line	Data
a5fb4114	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
a5fb4114	15
	16	//_________________________________________________________________________
	17	// Class containing methods for the isolation cut.
	18	// An AOD candidate (AliAODPWG4ParticleCorrelation type)
	19	// is passed. Look in a cone around the candidate and study
	20	// the hadronic activity inside to decide if the candidate is isolated
	21	//
	22	//
	23	//*-- Author: Gustavo Conesa (LNF-INFN)
	24
	25	//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
	26	//-Yaxian Mao (check the candidate particle is the leading particle or not at the same hemishere)
	27
	28	//////////////////////////////////////////////////////////////////////////////
	29
	30
	31	// --- ROOT system ---
a5fb4114	32	#include <TLorentzVector.h>
	33	#include <TObjArray.h>
	34
	35	// --- AliRoot system ---
	36	#include "AliIsolationCut.h"
	37	#include "AliAODPWG4ParticleCorrelation.h"
	38	#include "AliAODTrack.h"
	39	#include "AliVCluster.h"
	40	#include "AliCaloTrackReader.h"
	41	#include "AliMixedEvent.h"
ac5111f9	42	#include "AliCaloPID.h"
a5fb4114	43
	44	ClassImp(AliIsolationCut)
	45
c5693f62	46	//____________________________________
	47	AliIsolationCut::AliIsolationCut() :
	48	TObject(),
	49	fConeSize(0.),
	50	fPtThreshold(0.),
	51	fSumPtThreshold(0.),
	52	fPtFraction(0.),
	53	fICMethod(0),
	54	fPartInCone(0)
a5fb4114	55	{
	56	//default ctor
	57
	58	//Initialize parameters
	59	InitParameters();
c5693f62	60
a5fb4114	61	}
a5fb4114	62
c5693f62	63	//____________________________________________
a5fb4114	64	TString AliIsolationCut::GetICParametersList()
	65	{
	66	//Put data member values in string to keep in output container
	67
	68	TString parList ; //this will be list of parameters used for this analysis.
	69	const Int_t buffersize = 255;
	70	char onePar[buffersize] ;
	71
	72	snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
	73	parList+=onePar ;
	74	snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
	75	parList+=onePar ;
	76	snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
	77	parList+=onePar ;
	78	snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
	79	parList+=onePar ;
	80	snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
	81	parList+=onePar ;
	82	snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
	83	parList+=onePar ;
1a31a9ab	84
a5fb4114	85	return parList;
	86	}
	87
f18cb329	88	//____________________________________
a5fb4114	89	void AliIsolationCut::InitParameters()
	90	{
	91	//Initialize the parameters of the analysis.
	92
	93	fConeSize = 0.4 ;
	94	fPtThreshold = 1. ;
	95	fSumPtThreshold = 0.5 ;
	96	fPtFraction = 0.1 ;
	97	fPartInCone = kOnlyCharged;
	98	fICMethod = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
	99
	100	}
	101
c5693f62	102	//________________________________________________________________________________
	103	void AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS,
	104	const TObjArray * plNe,
	105	const AliCaloTrackReader * reader,
ac5111f9	106	const AliCaloPID * pid,
c5693f62	107	const Bool_t bFillAOD,
	108	AliAODPWG4ParticleCorrelation *pCandidate,
	109	const TString & aodArrayRefName,
b960c7eb	110	Int_t & n,
	111	Int_t & nfrac,
	112	Float_t & coneptsum,
	113	Bool_t & isolated) const
a5fb4114	114	{
	115	//Search in cone around a candidate particle if it is isolated
	116	Float_t phiC = pCandidate->Phi() ;
	117	if(phiC<0) phiC+=TMath::TwoPi();
	118	Float_t etaC = pCandidate->Eta() ;
	119	Float_t ptC = pCandidate->Pt() ;
	120	Float_t pt = -100. ;
	121	Float_t eta = -100. ;
	122	Float_t phi = -100. ;
	123	Float_t rad = -100. ;
	124
	125	n = 0 ;
	126	nfrac = 0 ;
	127	coneptsum = 0.;
	128	isolated = kFALSE;
1a31a9ab	129
a5fb4114	130	//Initialize the array with refrences
	131	TObjArray * refclusters = 0x0;
	132	TObjArray * reftracks = 0x0;
	133	Int_t ntrackrefs = 0;
	134	Int_t nclusterrefs = 0;
ac5111f9	135
a5fb4114	136	//Check charged particles in cone.
b960c7eb	137	if(plCTS &&
	138	(fPartInCone==kOnlyCharged \|\| fPartInCone==kNeutralAndCharged))
	139	{
a5fb4114	140	TVector3 p3;
b960c7eb	141	for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
	142	{
	143
a5fb4114	144	AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ;
b960c7eb	145
a5fb4114	146	//Do not count the candidate (pion, conversion photon) or the daughters of the candidate
b960c7eb	147	if(track->GetID() == pCandidate->GetTrackLabel(0) \|\| track->GetID() == pCandidate->GetTrackLabel(1) \|\|
	148	track->GetID() == pCandidate->GetTrackLabel(2) \|\| track->GetID() == pCandidate->GetTrackLabel(3) ) continue ;
	149
a5fb4114	150	p3.SetXYZ(track->Px(),track->Py(),track->Pz());
	151	pt = p3.Pt();
	152	eta = p3.Eta();
	153	phi = p3.Phi() ;
	154	if(phi<0) phi+=TMath::TwoPi();
	155
b960c7eb	156	// Only loop the particle at the same side of candidate
a5fb4114	157	if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
b960c7eb	158
	159	// If at the same side has particle larger than candidate,
	160	// then candidate can not be the leading, skip such events
	161	if(pt > ptC)
	162	{
a5fb4114	163	n = -1;
	164	nfrac = -1;
	165	coneptsum = -1;
	166	isolated = kFALSE;
b960c7eb	167	if(bFillAOD && reftracks)
b960c7eb	168	{
a5fb4114	169	reftracks->Clear();
	170	delete reftracks;
	171	}
b960c7eb	172
a5fb4114	173	return ;
a5fb4114	174	}
f18cb329	175
a5fb4114	176	//Check if there is any particle inside cone with pt larger than fPtThreshold
f18cb329	177
f18cb329	178	rad = Radius(etaC, phiC, eta, phi);
a5fb4114	179
b960c7eb	180	if(rad < fConeSize)
	181	{
	182	if(bFillAOD)
	183	{
a5fb4114	184	ntrackrefs++;
b960c7eb	185	if(ntrackrefs == 1)
b960c7eb	186	{
a5fb4114	187	reftracks = new TObjArray(0);
	188	//reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
	189	TString tempo(aodArrayRefName) ;
	190	tempo += "Tracks" ;
	191	reftracks->SetName(tempo);
	192	reftracks->SetOwner(kFALSE);
	193	}
	194	reftracks->Add(track);
	195	}
b960c7eb	196
a5fb4114	197	//printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
	198	coneptsum+=pt;
	199	if(pt > fPtThreshold ) n++;
	200	if(pt > fPtFraction*ptC ) nfrac++;
b960c7eb	201
a5fb4114	202	} // Inside cone
	203	}// charged particle loop
	204	}//Tracks
	205
	206	//Check neutral particles in cone.
b960c7eb	207	if(plNe &&
	208	(fPartInCone==kOnlyNeutral \|\| fPartInCone==kNeutralAndCharged))
	209	{
a5fb4114	210	TLorentzVector mom ;
b960c7eb	211
	212	for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
	213	{
a5fb4114	214	AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
	215
	216	//Get the index where the cluster comes, to retrieve the corresponding vertex
	217	Int_t evtIndex = 0 ;
b960c7eb	218	if (reader->GetMixedEvent())
a5fb4114	219	evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
b960c7eb	220
a5fb4114	221
a5fb4114	222	//Do not count the candidate (photon or pi0) or the daughters of the candidate
ac5111f9	223	if(calo->GetID() == pCandidate->GetCaloLabel(0) \|\|
b960c7eb	224	calo->GetID() == pCandidate->GetCaloLabel(1) ) continue ;
a5fb4114	225
ac5111f9	226	//Skip matched clusters with tracks
	227	if( pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
	228
	229	//Assume that come from vertex in straight line
	230	calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
a5fb4114	231
	232	pt = mom.Pt();
	233	eta = mom.Eta();
	234	phi = mom.Phi() ;
	235	if(phi<0) phi+=TMath::TwoPi();
a5fb4114	236
b960c7eb	237	// Only loop the particle at the same side of candidate
ac5111f9	238	if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
ac5111f9	239
b960c7eb	240	// If at the same side has particle larger than candidate,
	241	// then candidate can not be the leading, skip such events
	242	if(pt > ptC)
	243	{
a5fb4114	244	n = -1;
	245	nfrac = -1;
	246	coneptsum = -1;
	247	isolated = kFALSE;
b960c7eb	248
	249	if(bFillAOD)
	250	{
	251	if(reftracks)
	252	{
a5fb4114	253	reftracks ->Clear();
	254	delete reftracks;
	255	}
b960c7eb	256
	257	if(refclusters)
	258	{
a5fb4114	259	refclusters->Clear();
	260	delete refclusters;
	261	}
	262	}
	263	return ;
	264	}
	265
	266	//Check if there is any particle inside cone with pt larger than fPtThreshold
f18cb329	267
	268	rad = Radius(etaC, phiC, eta, phi);
	269
b960c7eb	270	if(rad < fConeSize)
	271	{
	272	if(bFillAOD)
	273	{
a5fb4114	274	nclusterrefs++;
b960c7eb	275	if(nclusterrefs==1)
b960c7eb	276	{
a5fb4114	277	refclusters = new TObjArray(0);
	278	//refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
	279	TString tempo(aodArrayRefName) ;
	280	tempo += "Clusters" ;
	281	refclusters->SetName(tempo);
	282	refclusters->SetOwner(kFALSE);
	283	}
	284	refclusters->Add(calo);
	285	}
b960c7eb	286
a5fb4114	287	//printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
a5fb4114	288	coneptsum+=pt;
b960c7eb	289	if(pt > fPtThreshold ) n++;
a5fb4114	290	//if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
b960c7eb	291	if(fPtFraction*ptC<fPtThreshold)
b960c7eb	292	{
1a31a9ab	293	if(pt>fPtThreshold) nfrac++ ;
a5fb4114	294	}
b960c7eb	295	else
b960c7eb	296	{
1a31a9ab	297	if(pt>fPtFraction*ptC) nfrac++;
a5fb4114	298	}
b960c7eb	299
a5fb4114	300	}//in cone
b960c7eb	301
a5fb4114	302	}// neutral particle loop
b960c7eb	303
a5fb4114	304	}//neutrals
1a31a9ab	305
a5fb4114	306	//printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
	307
	308	//Add reference arrays to AOD when filling AODs only
b960c7eb	309	if(bFillAOD)
b960c7eb	310	{
a5fb4114	311	if(refclusters) pCandidate->AddObjArray(refclusters);
	312	if(reftracks) pCandidate->AddObjArray(reftracks);
	313	}
b960c7eb	314
	315	//Check isolation, depending on selected isolation criteria
	316	if( fICMethod == kPtThresIC)
	317	{
a5fb4114	318	if(n==0) isolated = kTRUE ;
a5fb4114	319	}
b960c7eb	320	else if( fICMethod == kSumPtIC)
b960c7eb	321	{
a5fb4114	322	if(coneptsum < fSumPtThreshold)
	323	isolated = kTRUE ;
	324	}
b960c7eb	325	else if( fICMethod == kPtFracIC)
b960c7eb	326	{
a5fb4114	327	if(nfrac==0) isolated = kTRUE ;
a5fb4114	328	}
b960c7eb	329	else if( fICMethod == kSumPtFracIC)
b960c7eb	330	{
a5fb4114	331	//when the fPtFraction*ptC < fSumPtThreshold then consider the later case
	332	if(fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold) isolated = kTRUE ;
	333	if(fPtFractionptC > fSumPtThreshold && coneptsum < fPtFractionptC) isolated = kTRUE ;
	334	}
	335
	336	}
	337
c5693f62	338	//_____________________________________________________
a5fb4114	339	void AliIsolationCut::Print(const Option_t * opt) const
	340	{
	341
	342	//Print some relevant parameters set for the analysis
	343	if(! opt)
	344	return;
	345
	346	printf("** Print %s %s ** \n", GetName(), GetTitle() ) ;
	347
c5693f62	348	printf("IC method = %d\n", fICMethod ) ;
c5693f62	349	printf("Cone Size = %1.2f\n", fConeSize ) ;
a5fb4114	350	printf("pT threshold = %2.1f\n", fPtThreshold) ;
c5693f62	351	printf("pT fraction = %3.1f\n", fPtFraction ) ;
c5693f62	352	printf("particle type in cone = %d\n", fPartInCone ) ;
a5fb4114	353	printf(" \n") ;
	354
	355	}
f18cb329	356
	357	//___________________________________________________________________________
	358	Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC,
	359	const Float_t eta , const Float_t phi) const
	360	{
	361	// Calculate the distance to trigger from any particle
	362
	363	Float_t dEta = etaC-eta;
	364	Float_t dPhi = phiC-phi;
	365
	366	if(TMath::Abs(dPhi) >= TMath::Pi())
	367	dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
	368
	369	return TMath::Sqrt( dEtadEta + dPhidPhi );
	370
	371	}
	372
	373
	374