[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / 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.                  *
 **************************************************************************/
/* $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) 
//////////////////////////////////////////////////////////////////////////////
  
  
// --- ROOT system --- 
//#include <Riostream.h>
#include <TLorentzVector.h>
#include <TObjArray.h>

// --- AliRoot system --- 
#include "AliIsolationCut.h" 
#include "AliAODPWG4ParticleCorrelation.h"
#include "AliAODTrack.h"
#include "AliAODCaloCluster.h"
#include "AliCaloTrackReader.h"

ClassImp(AliIsolationCut)
  
//____________________________________________________________________________
  AliIsolationCut::AliIsolationCut() : 
    TObject(),
    fConeSize(0.),fPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)
 
{
  //default ctor
  
  //Initialize parameters
  InitParameters();

}
/*
//____________________________________________________________________________
AliIsolationCut::AliIsolationCut(const AliIsolationCut & g) : 
  TObject(g),
  fConeSize(g.fConeSize),
  fPtThreshold(g.fPtThreshold),
  fPtFraction(g.fPtFraction), 
  fICMethod(g.fICMethod)
{
  // cpy ctor
  
}

//_________________________________________________________________________
AliIsolationCut & AliIsolationCut::operator = (const AliIsolationCut & source)
{
  // assignment operator
  
  if(&source == this) return *this;
  
  fConeSize = source.fConeSize ;
  fPtThreshold = source.fPtThreshold ; 
  fICMethod = source.fICMethod ;
  fPtFraction = source.fPtFraction ;

  return *this;
  
}
*/
//____________________________________________________________________________
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.
  char onePar[255] ;
  
  sprintf(onePar,"--- AliIsolationCut ---\n") ;
  parList+=onePar ;	
  sprintf(onePar,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
  parList+=onePar ;
  sprintf(onePar,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
  parList+=onePar ;
  sprintf(onePar,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
  parList+=onePar ;
  sprintf(onePar,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
  parList+=onePar ;
  sprintf(onePar,"fPartInCone=%d \n",fPartInCone) ;
  parList+=onePar ;

  return parList; 
}

//____________________________________________________________________________
void AliIsolationCut::InitParameters()
{
  //Initialize the parameters of the analysis.
  
  fConeSize    = 0.4 ; 
  fPtThreshold = 1. ; 
  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 fillAOD, 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 ;
  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();
      
      //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(fillAOD) {
	  ntrackrefs++;
	  if(ntrackrefs == 1){
	    reftracks = new TObjArray(0);
	    reftracks->SetName(aodArrayRefName+"Tracks");
	    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++;  
      }
    }// charged particle loop
  }//Tracks
  
  //Check neutral particles in cone.  
  if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){
	  
	//Get vertex for photon momentum calculation
	Double_t vertex[]  = {0,0,0} ; //vertex ;
	Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;
	if(!reader->GetDataType()== AliCaloTrackReader::kMC) 
	{
		reader->GetVertex(vertex);
		//if(reader->GetSecondInputAODTree()) reader->GetSecondInputAODVertex(vertex2);
	}
    TLorentzVector mom ;
    for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
      AliAODCaloCluster * calo = (AliAODCaloCluster *)(plNe->At(ipr)) ;
      
      //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,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  
      
      pt   = mom.Pt();
      eta  = mom.Eta();
      phi  = mom.Phi() ;
      if(phi<0) phi+=TMath::TwoPi();
      
      //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(fillAOD) {
	  nclusterrefs++;
	  if(nclusterrefs==1){
	    refclusters = new TObjArray(0);
	    refclusters->SetName(aodArrayRefName+"Clusters");
	    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(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(fillAOD) {
    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 < fPtThreshold)
      isolated  =  kTRUE ;
  }
  else if( fICMethod == kPtFracIC){
    if(nfrac==0) isolated = kTRUE ;
  }
  else if( fICMethod == kSumPtFracIC){
    if(coneptsum < fPtFraction*ptC)
      isolated  =  kTRUE ;
  }

  //if(refclusters) delete refclusters;
  //if(reftracks)   delete reftracks;

}

//__________________________________________________________________
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") ;
  
}
Commit	Line	Data
1c5acb87	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	**************************************************************************/
	15	/* $Id: $ */
	16
	17	//_________________________________________________________________________
	18	// Class containing methods for the isolation cut.
591cc579	19	// An AOD candidate (AliAODPWG4ParticleCorrelation type)
	20	// is passed. Look in a cone around the candidate and study
	21	// the hadronic activity inside to decide if the candidate is isolated
1c5acb87	22	//
	23	//
	24	//*-- Author: Gustavo Conesa (LNF-INFN)
	25	//////////////////////////////////////////////////////////////////////////////
	26
	27
	28	// --- ROOT system ---
	29	//#include <Riostream.h>
	30	#include <TLorentzVector.h>
591cc579	31	#include <TObjArray.h>
1c5acb87	32
	33	// --- AliRoot system ---
	34	#include "AliIsolationCut.h"
1c5acb87	35	#include "AliAODPWG4ParticleCorrelation.h"
	36	#include "AliAODTrack.h"
	37	#include "AliAODCaloCluster.h"
233e0df8	38	#include "AliCaloTrackReader.h"
1c5acb87	39
	40	ClassImp(AliIsolationCut)
	41
	42	//____________________________________________________________________________
	43	AliIsolationCut::AliIsolationCut() :
	44	TObject(),
b6991fc4	45	fConeSize(0.),fPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)
1c5acb87	46
	47	{
	48	//default ctor
	49
	50	//Initialize parameters
	51	InitParameters();
	52
	53	}
78219bac	54	/*
1c5acb87	55	//____________________________________________________________________________
	56	AliIsolationCut::AliIsolationCut(const AliIsolationCut & g) :
	57	TObject(g),
	58	fConeSize(g.fConeSize),
	59	fPtThreshold(g.fPtThreshold),
	60	fPtFraction(g.fPtFraction),
	61	fICMethod(g.fICMethod)
	62	{
	63	// cpy ctor
477d6cee	64
1c5acb87	65	}
	66
	67	//_________________________________________________________________________
	68	AliIsolationCut & AliIsolationCut::operator = (const AliIsolationCut & source)
	69	{
	70	// assignment operator
	71
	72	if(&source == this) return *this;
477d6cee	73
1c5acb87	74	fConeSize = source.fConeSize ;
	75	fPtThreshold = source.fPtThreshold ;
	76	fICMethod = source.fICMethod ;
	77	fPtFraction = source.fPtFraction ;
	78
	79	return *this;
	80
	81	}
78219bac	82	*/
1c5acb87	83	//____________________________________________________________________________
	84	TString AliIsolationCut::GetICParametersList()
	85	{
477d6cee	86	//Put data member values in string to keep in output container
	87
	88	TString parList ; //this will be list of parameters used for this analysis.
	89	char onePar[255] ;
	90
	91	sprintf(onePar,"--- AliIsolationCut ---\n") ;
	92	parList+=onePar ;
	93	sprintf(onePar,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
	94	parList+=onePar ;
	95	sprintf(onePar,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
	96	parList+=onePar ;
	97	sprintf(onePar,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
	98	parList+=onePar ;
	99	sprintf(onePar,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
	100	parList+=onePar ;
b6991fc4	101	sprintf(onePar,"fPartInCone=%d \n",fPartInCone) ;
	102	parList+=onePar ;
	103
477d6cee	104	return parList;
1c5acb87	105	}
	106
	107	//____________________________________________________________________________
	108	void AliIsolationCut::InitParameters()
	109	{
	110	//Initialize the parameters of the analysis.
477d6cee	111
b6991fc4	112	fConeSize = 0.4 ;
	113	fPtThreshold = 1. ;
	114	fPtFraction = 0.1 ;
	115	fPartInCone = kNeutralAndCharged;
	116	fICMethod = kPtThresIC; // 0 pt threshol method, 1 cone pt sum method
477d6cee	117
1c5acb87	118	}
	119
	120	//__________________________________________________________________
233e0df8	121	void AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS, TObjArray * const plNe, AliCaloTrackReader * const reader,
1c5acb87	122	const Bool_t fillAOD, AliAODPWG4ParticleCorrelation *pCandidate,
a3aebfff	123	const TString aodArrayRefName,
1c5acb87	124	Int_t & n, Int_t & nfrac, Float_t &coneptsum, Bool_t &isolated) const
1c5acb87	125	{
477d6cee	126	//Search in cone around a candidate particle if it is isolated
477d6cee	127	Float_t phiC = pCandidate->Phi() ;
9ccc5a0b	128	Float_t etaC = pCandidate->Eta() ;
	129	Float_t ptC = pCandidate->Pt() ;
	130	Float_t pt = -100. ;
477d6cee	131	Float_t eta = -100. ;
9ccc5a0b	132	Float_t phi = -100. ;
477d6cee	133	Float_t rad = -100 ;
477d6cee	134	n = 0 ;
	135	coneptsum = 0.;
	136	isolated = kFALSE;
a3aebfff	137
a3aebfff	138	//Initialize the array with refrences
591cc579	139	TObjArray * refclusters = 0x0;
9ccc5a0b	140	TObjArray * reftracks = 0x0;
9ccc5a0b	141	Int_t ntrackrefs = 0;
4ba526c9	142	Int_t nclusterrefs = 0;
4ba526c9	143
477d6cee	144	//Check charged particles in cone.
b6991fc4	145	if(plCTS && (fPartInCone==kOnlyCharged \|\| fPartInCone==kNeutralAndCharged)){
477d6cee	146	TVector3 p3;
	147	for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
	148	AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ;
	149	//Do not count the candidate (pion, conversion photon) or the daughters of the candidate
	150	if(track->GetID() == pCandidate->GetTrackLabel(0) \|\| track->GetID() == pCandidate->GetTrackLabel(1)) continue ;
	151	p3.SetXYZ(track->Px(),track->Py(),track->Pz());
	152	pt = p3.Pt();
	153	eta = p3.Eta();
	154	phi = p3.Phi() ;
	155	if(phi<0) phi+=TMath::TwoPi();
	156
	157	//Check if there is any particle inside cone with pt larger than fPtThreshold
	158	rad = TMath::Sqrt((eta-etaC)(eta-etaC)+ (phi-phiC)(phi-phiC));
	159
	160	if(rad < fConeSize){
	161	if(fillAOD) {
4ba526c9	162	ntrackrefs++;
9ccc5a0b	163	if(ntrackrefs == 1){
9ccc5a0b	164	reftracks = new TObjArray(0);
4ba526c9	165	reftracks->SetName(aodArrayRefName+"Tracks");
	166	reftracks->SetOwner(kFALSE);
	167	}
a3aebfff	168	reftracks->Add(track);
1c5acb87	169	}
477d6cee	170	//printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
	171	coneptsum+=pt;
	172	if(pt > fPtThreshold ) n++;
	173	if(pt > fPtFraction*ptC ) nfrac++;
	174	}
	175	}// charged particle loop
	176	}//Tracks
	177
	178	//Check neutral particles in cone.
b6991fc4	179	if(plNe && (fPartInCone==kOnlyNeutral \|\| fPartInCone==kNeutralAndCharged)){
233e0df8	180
	181	//Get vertex for photon momentum calculation
	182	Double_t vertex[] = {0,0,0} ; //vertex ;
	183	Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;
	184	if(!reader->GetDataType()== AliCaloTrackReader::kMC)
	185	{
	186	reader->GetVertex(vertex);
1e68a3f4	187	//if(reader->GetSecondInputAODTree()) reader->GetSecondInputAODVertex(vertex2);
233e0df8	188	}
477d6cee	189	TLorentzVector mom ;
	190	for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
	191	AliAODCaloCluster * calo = (AliAODCaloCluster *)(plNe->At(ipr)) ;
	192
	193	//Do not count the candidate (photon or pi0) or the daughters of the candidate
	194	if(calo->GetID() == pCandidate->GetCaloLabel(0) \|\| calo->GetID() == pCandidate->GetCaloLabel(1)) continue ; //Skip matched clusters with tracks
	195
	196	if(calo->GetNTracksMatched() > 0) continue ;
4ba526c9	197	//Input from second AOD?
4ba526c9	198	Int_t input = 0;
1e68a3f4	199	// if (pCandidate->GetDetector() == "EMCAL" && reader->GetAODEMCALNormalInputEntries() <= ipr) input = 1 ;
1e68a3f4	200	// else if(pCandidate->GetDetector() == "PHOS" && reader->GetAODPHOSNormalInputEntries() <= ipr) input = 1;
4ba526c9	201
	202	//Get Momentum vector,
	203	if (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
	204	else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
	205
	206	pt = mom.Pt();
477d6cee	207	eta = mom.Eta();
	208	phi = mom.Phi() ;
	209	if(phi<0) phi+=TMath::TwoPi();
	210
	211	//Check if there is any particle inside cone with pt larger than fPtThreshold
	212	rad = TMath::Sqrt((eta-etaC)(eta-etaC)+ (phi-phiC)(phi-phiC));
	213	if(rad < fConeSize){
	214	if(fillAOD) {
4ba526c9	215	nclusterrefs++;
4ba526c9	216	if(nclusterrefs==1){
9ccc5a0b	217	refclusters = new TObjArray(0);
4ba526c9	218	refclusters->SetName(aodArrayRefName+"Clusters");
9ccc5a0b	219	refclusters->SetOwner(kFALSE);
4ba526c9	220	}
a3aebfff	221	refclusters->Add(calo);
1c5acb87	222	}
477d6cee	223	//printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
	224	coneptsum+=pt;
	225	if(pt > fPtThreshold ) n++;
	226	if(pt > fPtFraction*ptC ) nfrac++;
	227	}//in cone
	228	}// neutral particle loop
	229	}//neutrals
	230
	231	//printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
	232
a3aebfff	233	//Add reference arrays to AOD when filling AODs only
a3aebfff	234	if(fillAOD) {
4ba526c9	235	if(refclusters) pCandidate->AddObjArray(refclusters);
4ba526c9	236	if(reftracks) pCandidate->AddObjArray(reftracks);
a3aebfff	237	}
a3aebfff	238
477d6cee	239	//Check isolation, depending on method.
	240	if( fICMethod == kPtThresIC){
	241	if(n==0) isolated = kTRUE ;
	242	}
	243	else if( fICMethod == kSumPtIC){
	244	if(coneptsum < fPtThreshold)
	245	isolated = kTRUE ;
	246	}
	247	else if( fICMethod == kPtFracIC){
	248	if(nfrac==0) isolated = kTRUE ;
	249	}
	250	else if( fICMethod == kSumPtFracIC){
	251	if(coneptsum < fPtFraction*ptC)
	252	isolated = kTRUE ;
	253	}
4ba526c9	254
9ccc5a0b	255	//if(refclusters) delete refclusters;
4ba526c9	256	//if(reftracks) delete reftracks;
4ba526c9	257
1c5acb87	258	}
	259
	260	//__________________________________________________________________
	261	void AliIsolationCut::Print(const Option_t * opt) const
	262	{
	263
	264	//Print some relevant parameters set for the analysis
	265	if(! opt)
	266	return;
	267
	268	printf("** Print %s %s ** \n", GetName(), GetTitle() ) ;
	269
	270	printf("IC method = %d\n", fICMethod) ;
	271	printf("Cone Size = %1.2f\n", fConeSize) ;
	272	printf("pT threshold = %2.1f\n", fPtThreshold) ;
	273	printf("pT fraction = %3.1f\n", fPtFraction) ;
b6991fc4	274	printf("particle type in cone = %d\n",fPartInCone);
1c5acb87	275	printf(" \n") ;
	276
	277	}