[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 "AliEMCALGeometry.h"
#include "AliEMCALGeoParams.h"
#include "AliCalorimeterUtils.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),
fDebug(-1),
fFracIsThresh(1)
{
  //default ctor
  
  //Initialize parameters
  InitParameters();
  
}


//__________________________________________________________________________________
Float_t AliIsolationCut::GetCellDensity(const AliAODPWG4ParticleCorrelation * pCandidate, 
                                        const AliCaloTrackReader * reader) const 
{
  // Get good cell density (number of active cells over all cells in cone)
  
  Double_t coneCells    = 0.; //number of cells in cone with radius fConeSize
  Double_t coneCellsBad = 0.; //number of bad cells in cone with radius fConeSize
  Double_t cellDensity  = 1.;

  Float_t phiC  = pCandidate->Phi() ;
  if(phiC<0) phiC+=TMath::TwoPi();
  Float_t etaC  = pCandidate->Eta() ;
  
  if(pCandidate->GetDetector()=="EMCAL")
  {
    AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
    AliCalorimeterUtils *cu = reader->GetCaloUtils();
    
    Int_t absId = -999;
    if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
    {
      //Get absolute (col,row) of candidate
      Int_t iEta=-1, iPhi=-1, iRCU = -1;      
      Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
      
      Int_t colC = iEta;
      if (nSupMod % 2) colC =  AliEMCALGeoParams::fgkEMCALCols + iEta ;
      Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
      
      Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
      //loop on cells in a square of side fConeSize to check cells in cone    
      for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
      {
        for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
        {
          if (Radius(colC, rowC, icol, irow) < sqrSize)
          {
            coneCells += 1.;
            
            Int_t cellSM  = -999;
            Int_t cellEta = -999;
            Int_t cellPhi = -999;
            if(icol > AliEMCALGeoParams::fgkEMCALCols-1) 
            {
              cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
              cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
              cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
            }
            if(icol < AliEMCALGeoParams::fgkEMCALCols) 
            {
              cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
              cellEta = icol;
              cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
            }
            
            //Count as bad "cells" out of EMCAL acceptance
            if(icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2 || 
               irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*16./3) //5*nRows+1/3*nRows
            {
              coneCellsBad += 1.;
            }
            //Count as bad "cells" marked as bad in the DataBase
            else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1) 
            {
              coneCellsBad += 1. ;
            }
          }
        }
      }//end of cells loop
    }
    
    else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
    
    if (coneCells > 0.) 
    {
      cellDensity = (coneCells-coneCellsBad)/coneCells;
      //printf("Energy density = %f\n", cellDensity);
    }
  }

  return cellDensity;
  
}

//____________________________________________
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 ;
 snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
  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
  fFracIsThresh   = 1; 
}

//________________________________________________________________________________
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 ptC   = pCandidate->Pt() ;
  Float_t phiC  = pCandidate->Phi() ;
  if(phiC<0) phiC+=TMath::TwoPi();
  Float_t etaC  = pCandidate->Eta() ;
  Float_t pt    = -100. ;
  Float_t eta   = -100. ;
  Float_t phi   = -100. ;
  Float_t rad   = -100. ;
  
  n         = 0 ;
  nfrac     = 0 ;
  coneptsum = 0.; 
  isolated  = kFALSE;

  if(fDebug>0) 
  {
    printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
           pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
    if(plCTS) printf(", nTracks %d"  ,plCTS->GetEntriesFast());
    if(plNe)  printf(", nClusters %d",plNe ->GetEntriesFast());
    
    printf("\n");
  }
  
  //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 ++ )
    {
      AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ; 
      
      if(track)
      {
        //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() ;
      }
      else
      {// Mixed event stored in AliAODPWG4Particles
        AliAODPWG4Particle * trackmix = dynamic_cast<AliAODPWG4Particle*>(plCTS->At(ipr)) ; 
        if(!trackmix)
        {
          printf("AliIsolationCut::MakeIsolationCut() - Wrong track data type, continue\n");
          continue;
        }
        
        pt  = trackmix->Pt();
        eta = trackmix->Eta();
        phi = trackmix->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;
      
        pCandidate->SetLeadingParticle(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(fDebug>0) 
        printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
               
      if(rad < fConeSize)
      {
        if(fDebug>0)  printf(" -  inside candidate cone");

        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);
        }
        
        
        coneptsum+=pt;
        if(pt > fPtThreshold )    n++;
        if(pt > fPtFraction*ptC ) nfrac++;  
        
      } // Inside cone

      if(fDebug>0)  printf("\n");

    }// 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 = dynamic_cast<AliVCluster *>(plNe->At(ipr)) ;
      
      if(calo)
      {
        //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 in case of neutral+charged analysis
        if( fPartInCone == kNeutralAndCharged && 
           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() ;
      }
      else 
      {// Mixed event stored in AliAODPWG4Particles
        AliAODPWG4Particle * calomix = dynamic_cast<AliAODPWG4Particle*>(plNe->At(ipr)) ; 
        if(!calomix)
        {
          printf("AliIsolationCut::MakeIsolationCut() - Wrong calo data type, continue\n");
          continue;
        }
        
        pt  = calomix->Pt();
        eta = calomix->Eta();
        phi = calomix->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;
        
        pCandidate->SetLeadingParticle(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(fDebug>0) 
        printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
      
      if(rad < fConeSize)
      {
        if(fDebug>0)  printf(" - inside candidate cone");

        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);
        }
        
        coneptsum+=pt;
        if(pt > fPtThreshold )   n++;
        //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
        if(fFracIsThresh){
	  if( fPtFraction*ptC<fPtThreshold)
	    {
	      if(pt>fPtThreshold)    nfrac++ ;
	    }
	  else 
	    {
	     if(pt>fPtFraction*ptC) nfrac++; 
	    }
	}
	else {
	  if(pt>fPtFraction*ptC) nfrac++;   
	}
        
      }//in cone
      
      if(fDebug>0)  printf("\n");

    }// neutral particle loop
  
  }//neutrals
  
  
  //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(fFracIsThresh ){
      if(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
      if( fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
    }
    else 
      {
 	if(coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
      }
  }
 else if( fICMethod == kSumDensityIC)
  {    
    // Get good cell density (number of active cells over all cells in cone)
    // and correct energy in cone
    Float_t cellDensity = GetCellDensity(pCandidate,reader);
    if(coneptsum < fSumPtThreshold*cellDensity)
      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("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
  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"
f8d07abf	38	#include "AliEMCALGeometry.h"
	39	#include "AliEMCALGeoParams.h"
	40	#include "AliCalorimeterUtils.h"
a5fb4114	41	#include "AliAODTrack.h"
	42	#include "AliVCluster.h"
	43	#include "AliCaloTrackReader.h"
	44	#include "AliMixedEvent.h"
ac5111f9	45	#include "AliCaloPID.h"
a5fb4114	46
	47	ClassImp(AliIsolationCut)
	48
c5693f62	49	//____________________________________
	50	AliIsolationCut::AliIsolationCut() :
	51	TObject(),
	52	fConeSize(0.),
	53	fPtThreshold(0.),
	54	fSumPtThreshold(0.),
	55	fPtFraction(0.),
	56	fICMethod(0),
03bae431	57	fPartInCone(0),
22c7f802	58	fDebug(-1),
22c7f802	59	fFracIsThresh(1)
a5fb4114	60	{
	61	//default ctor
	62
	63	//Initialize parameters
	64	InitParameters();
c5693f62	65
a5fb4114	66	}
a5fb4114	67
f8d07abf	68
	69	//__________________________________________________________________________________
	70	Float_t AliIsolationCut::GetCellDensity(const AliAODPWG4ParticleCorrelation * pCandidate,
	71	const AliCaloTrackReader * reader) const
	72	{
	73	// Get good cell density (number of active cells over all cells in cone)
	74
	75	Double_t coneCells = 0.; //number of cells in cone with radius fConeSize
	76	Double_t coneCellsBad = 0.; //number of bad cells in cone with radius fConeSize
	77	Double_t cellDensity = 1.;
	78
	79	Float_t phiC = pCandidate->Phi() ;
	80	if(phiC<0) phiC+=TMath::TwoPi();
	81	Float_t etaC = pCandidate->Eta() ;
	82
	83	if(pCandidate->GetDetector()=="EMCAL")
	84	{
	85	AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
	86	AliCalorimeterUtils *cu = reader->GetCaloUtils();
	87
	88	Int_t absId = -999;
	89	if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
	90	{
	91	//Get absolute (col,row) of candidate
	92	Int_t iEta=-1, iPhi=-1, iRCU = -1;
	93	Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
	94
	95	Int_t colC = iEta;
	96	if (nSupMod % 2) colC = AliEMCALGeoParams::fgkEMCALCols + iEta ;
	97	Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
	98
	99	Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
	100	//loop on cells in a square of side fConeSize to check cells in cone
	101	for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
	102	{
	103	for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
	104	{
	105	if (Radius(colC, rowC, icol, irow) < sqrSize)
	106	{
	107	coneCells += 1.;
	108
	109	Int_t cellSM = -999;
	110	Int_t cellEta = -999;
	111	Int_t cellPhi = -999;
	112	if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
	113	{
	114	cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
	115	cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
	116	cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
	117	}
	118	if(icol < AliEMCALGeoParams::fgkEMCALCols)
	119	{
	120	cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
	121	cellEta = icol;
	122	cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
	123	}
	124
	125	//Count as bad "cells" out of EMCAL acceptance
	126	if(icol < 0 \|\| icol > AliEMCALGeoParams::fgkEMCALCols*2 \|\|
	127	irow < 0 \|\| irow > AliEMCALGeoParams::fgkEMCALRows16./3) //5nRows+1/3*nRows
	128	{
	129	coneCellsBad += 1.;
	130	}
	131	//Count as bad "cells" marked as bad in the DataBase
132	else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1)
133	{
134	coneCellsBad += 1. ;
135	}
136	}
137	}
138	}//end of cells loop
139	}
140
141	else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
142
143	if (coneCells > 0.)
144	{
145	cellDensity = (coneCells-coneCellsBad)/coneCells;
8c00ddb0	146	//printf("Energy density = %f\n", cellDensity);
f8d07abf	147	}
	148	}
	149
	150	return cellDensity;
	151
	152	}
	153
c5693f62	154	//____________________________________________
a5fb4114	155	TString AliIsolationCut::GetICParametersList()
	156	{
	157	//Put data member values in string to keep in output container
	158
	159	TString parList ; //this will be list of parameters used for this analysis.
	160	const Int_t buffersize = 255;
	161	char onePar[buffersize] ;
	162
	163	snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
	164	parList+=onePar ;
	165	snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
	166	parList+=onePar ;
	167	snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
	168	parList+=onePar ;
	169	snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
	170	parList+=onePar ;
	171	snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
	172	parList+=onePar ;
	173	snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
	174	parList+=onePar ;
22c7f802	175	snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
	176	parList+=onePar ;
	177
a5fb4114	178	return parList;
	179	}
	180
f18cb329	181	//____________________________________
a5fb4114	182	void AliIsolationCut::InitParameters()
	183	{
	184	//Initialize the parameters of the analysis.
	185
	186	fConeSize = 0.4 ;
	187	fPtThreshold = 1. ;
	188	fSumPtThreshold = 0.5 ;
	189	fPtFraction = 0.1 ;
	190	fPartInCone = kOnlyCharged;
	191	fICMethod = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
22c7f802	192	fFracIsThresh = 1;
a5fb4114	193	}
a5fb4114	194
c5693f62	195	//________________________________________________________________________________
	196	void AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS,
	197	const TObjArray * plNe,
	198	const AliCaloTrackReader * reader,
ac5111f9	199	const AliCaloPID * pid,
c5693f62	200	const Bool_t bFillAOD,
	201	AliAODPWG4ParticleCorrelation *pCandidate,
	202	const TString & aodArrayRefName,
b960c7eb	203	Int_t & n,
	204	Int_t & nfrac,
	205	Float_t & coneptsum,
	206	Bool_t & isolated) const
a5fb4114	207	{
a5fb4114	208	//Search in cone around a candidate particle if it is isolated
03bae431	209	Float_t ptC = pCandidate->Pt() ;
a5fb4114	210	Float_t phiC = pCandidate->Phi() ;
	211	if(phiC<0) phiC+=TMath::TwoPi();
	212	Float_t etaC = pCandidate->Eta() ;
a5fb4114	213	Float_t pt = -100. ;
	214	Float_t eta = -100. ;
	215	Float_t phi = -100. ;
	216	Float_t rad = -100. ;
	217
	218	n = 0 ;
	219	nfrac = 0 ;
	220	coneptsum = 0.;
	221	isolated = kFALSE;
03bae431	222
	223	if(fDebug>0)
	224	{
	225	printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
	226	pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
	227	if(plCTS) printf(", nTracks %d" ,plCTS->GetEntriesFast());
	228	if(plNe) printf(", nClusters %d",plNe ->GetEntriesFast());
	229
	230	printf("\n");
	231	}
1a31a9ab	232
a5fb4114	233	//Initialize the array with refrences
a258315d	234	TObjArray * refclusters = 0x0;
	235	TObjArray * reftracks = 0x0;
	236	Int_t ntrackrefs = 0;
	237	Int_t nclusterrefs = 0;
ac5111f9	238
a5fb4114	239	//Check charged particles in cone.
b960c7eb	240	if(plCTS &&
	241	(fPartInCone==kOnlyCharged \|\| fPartInCone==kNeutralAndCharged))
	242	{
a5fb4114	243	TVector3 p3;
b960c7eb	244	for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
b960c7eb	245	{
a258315d	246	AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ;
b960c7eb	247
a258315d	248	if(track)
	249	{
	250	//Do not count the candidate (pion, conversion photon) or the daughters of the candidate
	251	if(track->GetID() == pCandidate->GetTrackLabel(0) \|\| track->GetID() == pCandidate->GetTrackLabel(1) \|\|
	252	track->GetID() == pCandidate->GetTrackLabel(2) \|\| track->GetID() == pCandidate->GetTrackLabel(3) ) continue ;
	253
	254	p3.SetXYZ(track->Px(),track->Py(),track->Pz());
	255	pt = p3.Pt();
	256	eta = p3.Eta();
	257	phi = p3.Phi() ;
	258	}
	259	else
	260	{// Mixed event stored in AliAODPWG4Particles
	261	AliAODPWG4Particle * trackmix = dynamic_cast<AliAODPWG4Particle*>(plCTS->At(ipr)) ;
	262	if(!trackmix)
	263	{
	264	printf("AliIsolationCut::MakeIsolationCut() - Wrong track data type, continue\n");
	265	continue;
	266	}
	267
	268	pt = trackmix->Pt();
	269	eta = trackmix->Eta();
	270	phi = trackmix->Phi() ;
	271	}
b960c7eb	272
a258315d	273	if( phi < 0 ) phi+=TMath::TwoPi();
a5fb4114	274
b960c7eb	275	// Only loop the particle at the same side of candidate
a5fb4114	276	if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
b960c7eb	277
	278	// If at the same side has particle larger than candidate,
	279	// then candidate can not be the leading, skip such events
	280	if(pt > ptC)
	281	{
a5fb4114	282	n = -1;
	283	nfrac = -1;
	284	coneptsum = -1;
	285	isolated = kFALSE;
3f150b4b	286
	287	pCandidate->SetLeadingParticle(kFALSE);
	288
b960c7eb	289	if(bFillAOD && reftracks)
b960c7eb	290	{
a5fb4114	291	reftracks->Clear();
	292	delete reftracks;
	293	}
b960c7eb	294
a5fb4114	295	return ;
a5fb4114	296	}
f18cb329	297
a5fb4114	298	//Check if there is any particle inside cone with pt larger than fPtThreshold
f18cb329	299
f18cb329	300	rad = Radius(etaC, phiC, eta, phi);
a5fb4114	301
03bae431	302	if(fDebug>0)
	303	printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
	304
b960c7eb	305	if(rad < fConeSize)
b960c7eb	306	{
03bae431	307	if(fDebug>0) printf(" - inside candidate cone");
03bae431	308
b960c7eb	309	if(bFillAOD)
b960c7eb	310	{
a5fb4114	311	ntrackrefs++;
b960c7eb	312	if(ntrackrefs == 1)
b960c7eb	313	{
a5fb4114	314	reftracks = new TObjArray(0);
	315	//reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
	316	TString tempo(aodArrayRefName) ;
	317	tempo += "Tracks" ;
	318	reftracks->SetName(tempo);
	319	reftracks->SetOwner(kFALSE);
	320	}
	321	reftracks->Add(track);
	322	}
b960c7eb	323
03bae431	324
a5fb4114	325	coneptsum+=pt;
	326	if(pt > fPtThreshold ) n++;
	327	if(pt > fPtFraction*ptC ) nfrac++;
b960c7eb	328
a5fb4114	329	} // Inside cone
03bae431	330
	331	if(fDebug>0) printf("\n");
	332
a5fb4114	333	}// charged particle loop
03bae431	334
03bae431	335
a5fb4114	336	}//Tracks
a5fb4114	337
03bae431	338
a5fb4114	339	//Check neutral particles in cone.
b960c7eb	340	if(plNe &&
	341	(fPartInCone==kOnlyNeutral \|\| fPartInCone==kNeutralAndCharged))
	342	{
a5fb4114	343	TLorentzVector mom ;
b960c7eb	344
	345	for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
	346	{
a258315d	347	AliVCluster * calo = dynamic_cast<AliVCluster *>(plNe->At(ipr)) ;
b960c7eb	348
a258315d	349	if(calo)
	350	{
	351	//Get the index where the cluster comes, to retrieve the corresponding vertex
	352	Int_t evtIndex = 0 ;
	353	if (reader->GetMixedEvent())
	354	evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
	355
	356
	357	//Do not count the candidate (photon or pi0) or the daughters of the candidate
	358	if(calo->GetID() == pCandidate->GetCaloLabel(0) \|\|
	359	calo->GetID() == pCandidate->GetCaloLabel(1) ) continue ;
	360
	361	//Skip matched clusters with tracks in case of neutral+charged analysis
	362	if( fPartInCone == kNeutralAndCharged &&
	363	pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
	364
	365	//Assume that come from vertex in straight line
	366	calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
	367
	368	pt = mom.Pt() ;
	369	eta = mom.Eta() ;
	370	phi = mom.Phi() ;
	371	}
	372	else
	373	{// Mixed event stored in AliAODPWG4Particles
	374	AliAODPWG4Particle * calomix = dynamic_cast<AliAODPWG4Particle*>(plNe->At(ipr)) ;
	375	if(!calomix)
	376	{
	377	printf("AliIsolationCut::MakeIsolationCut() - Wrong calo data type, continue\n");
	378	continue;
	379	}
	380
	381	pt = calomix->Pt();
	382	eta = calomix->Eta();
	383	phi = calomix->Phi() ;
	384	}
a5fb4114	385
a258315d	386	if( phi < 0 ) phi+=TMath::TwoPi();
a5fb4114	387
a5fb4114	388
b960c7eb	389	// Only loop the particle at the same side of candidate
ac5111f9	390	if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
ac5111f9	391
b960c7eb	392	// If at the same side has particle larger than candidate,
	393	// then candidate can not be the leading, skip such events
	394	if(pt > ptC)
	395	{
a5fb4114	396	n = -1;
	397	nfrac = -1;
	398	coneptsum = -1;
	399	isolated = kFALSE;
b960c7eb	400
3f150b4b	401	pCandidate->SetLeadingParticle(kFALSE);
3f150b4b	402
b960c7eb	403	if(bFillAOD)
	404	{
	405	if(reftracks)
	406	{
a5fb4114	407	reftracks ->Clear();
	408	delete reftracks;
	409	}
b960c7eb	410
	411	if(refclusters)
	412	{
a5fb4114	413	refclusters->Clear();
	414	delete refclusters;
	415	}
	416	}
	417	return ;
	418	}
	419
	420	//Check if there is any particle inside cone with pt larger than fPtThreshold
f18cb329	421
	422	rad = Radius(etaC, phiC, eta, phi);
	423
03bae431	424	if(fDebug>0)
	425	printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
	426
b960c7eb	427	if(rad < fConeSize)
b960c7eb	428	{
03bae431	429	if(fDebug>0) printf(" - inside candidate cone");
03bae431	430
b960c7eb	431	if(bFillAOD)
b960c7eb	432	{
a5fb4114	433	nclusterrefs++;
b960c7eb	434	if(nclusterrefs==1)
b960c7eb	435	{
a5fb4114	436	refclusters = new TObjArray(0);
	437	//refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
	438	TString tempo(aodArrayRefName) ;
	439	tempo += "Clusters" ;
	440	refclusters->SetName(tempo);
	441	refclusters->SetOwner(kFALSE);
	442	}
	443	refclusters->Add(calo);
	444	}
b960c7eb	445
a5fb4114	446	coneptsum+=pt;
b960c7eb	447	if(pt > fPtThreshold ) n++;
a5fb4114	448	//if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
22c7f802	449	if(fFracIsThresh){
	450	if( fPtFraction*ptC<fPtThreshold)
	451	{
	452	if(pt>fPtThreshold) nfrac++ ;
	453	}
	454	else
	455	{
	456	if(pt>fPtFraction*ptC) nfrac++;
	457	}
	458	}
	459	else {
	460	if(pt>fPtFraction*ptC) nfrac++;
	461	}
b960c7eb	462
a5fb4114	463	}//in cone
b960c7eb	464
03bae431	465	if(fDebug>0) printf("\n");
03bae431	466
a5fb4114	467	}// neutral particle loop
b960c7eb	468
a5fb4114	469	}//neutrals
1a31a9ab	470
a5fb4114	471
a5fb4114	472	//Add reference arrays to AOD when filling AODs only
b960c7eb	473	if(bFillAOD)
b960c7eb	474	{
a5fb4114	475	if(refclusters) pCandidate->AddObjArray(refclusters);
	476	if(reftracks) pCandidate->AddObjArray(reftracks);
	477	}
b960c7eb	478
	479	//Check isolation, depending on selected isolation criteria
	480	if( fICMethod == kPtThresIC)
	481	{
a5fb4114	482	if(n==0) isolated = kTRUE ;
a5fb4114	483	}
b960c7eb	484	else if( fICMethod == kSumPtIC)
b960c7eb	485	{
a5fb4114	486	if(coneptsum < fSumPtThreshold)
	487	isolated = kTRUE ;
	488	}
b960c7eb	489	else if( fICMethod == kPtFracIC)
b960c7eb	490	{
a5fb4114	491	if(nfrac==0) isolated = kTRUE ;
a5fb4114	492	}
b960c7eb	493	else if( fICMethod == kSumPtFracIC)
b960c7eb	494	{
a5fb4114	495	//when the fPtFraction*ptC < fSumPtThreshold then consider the later case
22c7f802	496	if(fFracIsThresh ){
	497	if(fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold) isolated = kTRUE ;
	498	if( fPtFractionptC > fSumPtThreshold && coneptsum < fPtFractionptC) isolated = kTRUE ;
	499	}
	500	else
	501	{
	502	if(coneptsum < fPtFraction*ptC) isolated = kTRUE ;
	503	}
a5fb4114	504	}
22c7f802	505	else if( fICMethod == kSumDensityIC)
f8d07abf	506	{
	507	// Get good cell density (number of active cells over all cells in cone)
	508	// and correct energy in cone
	509	Float_t cellDensity = GetCellDensity(pCandidate,reader);
	510	if(coneptsum < fSumPtThreshold*cellDensity)
	511	isolated = kTRUE;
	512	}
a5fb4114	513
	514	}
	515
c5693f62	516	//_____________________________________________________
a5fb4114	517	void AliIsolationCut::Print(const Option_t * opt) const
	518	{
	519
	520	//Print some relevant parameters set for the analysis
	521	if(! opt)
	522	return;
	523
	524	printf("** Print %s %s ** \n", GetName(), GetTitle() ) ;
	525
c5693f62	526	printf("IC method = %d\n", fICMethod ) ;
c5693f62	527	printf("Cone Size = %1.2f\n", fConeSize ) ;
a5fb4114	528	printf("pT threshold = %2.1f\n", fPtThreshold) ;
c5693f62	529	printf("pT fraction = %3.1f\n", fPtFraction ) ;
c5693f62	530	printf("particle type in cone = %d\n", fPartInCone ) ;
22c7f802	531	printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
a5fb4114	532	printf(" \n") ;
	533
	534	}
f18cb329	535
	536	//___________________________________________________________________________
	537	Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC,
	538	const Float_t eta , const Float_t phi) const
	539	{
	540	// Calculate the distance to trigger from any particle
	541
	542	Float_t dEta = etaC-eta;
	543	Float_t dPhi = phiC-phi;
	544
	545	if(TMath::Abs(dPhi) >= TMath::Pi())
	546	dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
	547
	548	return TMath::Sqrt( dEtadEta + dPhidPhi );
	549
	550	}
	551
	552
	553