1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
23 //*-- Author: Gustavo Conesa (LNF-INFN)
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)
28 //////////////////////////////////////////////////////////////////////////////
31 // --- ROOT system ---
32 #include <TLorentzVector.h>
33 #include <TObjArray.h>
35 // --- AliRoot system ---
36 #include "AliIsolationCut.h"
37 #include "AliAODPWG4ParticleCorrelation.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliEMCALGeoParams.h"
40 #include "AliCalorimeterUtils.h"
41 #include "AliAODTrack.h"
42 #include "AliVCluster.h"
43 #include "AliCaloTrackReader.h"
44 #include "AliMixedEvent.h"
45 #include "AliCaloPID.h"
47 ClassImp(AliIsolationCut)
49 //____________________________________
50 AliIsolationCut::AliIsolationCut() :
54 fPtThresholdMax(10000.),
64 //Initialize parameters
70 //__________________________________________________________________________________
71 Float_t AliIsolationCut::GetCellDensity(const AliAODPWG4ParticleCorrelation * pCandidate,
72 const AliCaloTrackReader * reader) const
74 // Get good cell density (number of active cells over all cells in cone)
76 Double_t coneCells = 0.; //number of cells in cone with radius fConeSize
77 Double_t coneCellsBad = 0.; //number of bad cells in cone with radius fConeSize
78 Double_t cellDensity = 1.;
80 Float_t phiC = pCandidate->Phi() ;
81 if(phiC<0) phiC+=TMath::TwoPi();
82 Float_t etaC = pCandidate->Eta() ;
84 if(pCandidate->GetDetector()=="EMCAL")
86 AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
87 AliCalorimeterUtils *cu = reader->GetCaloUtils();
90 if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
92 //Get absolute (col,row) of candidate
93 Int_t iEta=-1, iPhi=-1, iRCU = -1;
94 Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
97 if (nSupMod % 2) colC = AliEMCALGeoParams::fgkEMCALCols + iEta ;
98 Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
100 Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
101 //loop on cells in a square of side fConeSize to check cells in cone
102 for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
104 for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
106 if (Radius(colC, rowC, icol, irow) < sqrSize)
111 Int_t cellEta = -999;
112 Int_t cellPhi = -999;
113 if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
115 cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
116 cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
117 cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
119 if(icol < AliEMCALGeoParams::fgkEMCALCols)
121 cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
123 cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
126 //Count as bad "cells" out of EMCAL acceptance
127 if(icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2 ||
128 irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*16./3) //5*nRows+1/3*nRows
132 //Count as bad "cells" marked as bad in the DataBase
133 else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1)
142 else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
146 cellDensity = (coneCells-coneCellsBad)/coneCells;
147 //printf("Energy density = %f\n", cellDensity);
155 //____________________________________________
156 TString AliIsolationCut::GetICParametersList()
158 //Put data member values in string to keep in output container
160 TString parList ; //this will be list of parameters used for this analysis.
161 const Int_t buffersize = 255;
162 char onePar[buffersize] ;
164 snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
166 snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
168 snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
170 snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
172 snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
174 snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
176 snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
182 //____________________________________
183 void AliIsolationCut::InitParameters()
185 //Initialize the parameters of the analysis.
189 fPtThresholdMax = 10000. ;
190 fSumPtThreshold = 0.5 ;
192 fPartInCone = kOnlyCharged;
193 fICMethod = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
197 //________________________________________________________________________________
198 void AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS,
199 const TObjArray * plNe,
200 const AliCaloTrackReader * reader,
201 const AliCaloPID * pid,
202 const Bool_t bFillAOD,
203 AliAODPWG4ParticleCorrelation *pCandidate,
204 const TString & aodArrayRefName,
208 Bool_t & isolated) const
210 //Search in cone around a candidate particle if it is isolated
211 Float_t ptC = pCandidate->Pt() ;
212 Float_t phiC = pCandidate->Phi() ;
213 if(phiC<0) phiC+=TMath::TwoPi();
214 Float_t etaC = pCandidate->Eta() ;
216 Float_t eta = -100. ;
217 Float_t phi = -100. ;
218 Float_t rad = -100. ;
227 printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
228 pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
229 if(plCTS) printf(", nTracks %d" ,plCTS->GetEntriesFast());
230 if(plNe) printf(", nClusters %d",plNe ->GetEntriesFast());
235 //Initialize the array with refrences
236 TObjArray * refclusters = 0x0;
237 TObjArray * reftracks = 0x0;
238 Int_t ntrackrefs = 0;
239 Int_t nclusterrefs = 0;
241 //Check charged particles in cone.
243 (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged))
246 for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
248 AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ;
252 //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
253 if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
254 track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3) ) continue ;
256 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
262 {// Mixed event stored in AliAODPWG4Particles
263 AliAODPWG4Particle * trackmix = dynamic_cast<AliAODPWG4Particle*>(plCTS->At(ipr)) ;
266 printf("AliIsolationCut::MakeIsolationCut() - Wrong track data type, continue\n");
271 eta = trackmix->Eta();
272 phi = trackmix->Phi() ;
275 if( phi < 0 ) phi+=TMath::TwoPi();
277 // Only loop the particle at the same side of candidate
278 if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
280 // If at the same side has particle larger than candidate,
281 // then candidate can not be the leading, skip such events
289 pCandidate->SetLeadingParticle(kFALSE);
291 if(bFillAOD && reftracks)
300 //Check if there is any particle inside cone with pt larger than fPtThreshold
302 rad = Radius(etaC, phiC, eta, phi);
305 printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
309 if(fDebug>0) printf(" - inside candidate cone");
316 reftracks = new TObjArray(0);
317 //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
318 TString tempo(aodArrayRefName) ;
320 reftracks->SetName(tempo);
321 reftracks->SetOwner(kFALSE);
323 reftracks->Add(track);
328 if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
329 if(pt > fPtFraction*ptC ) nfrac++;
333 if(fDebug>0) printf("\n");
335 }// charged particle loop
341 //Check neutral particles in cone.
343 (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged))
347 for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
349 AliVCluster * calo = dynamic_cast<AliVCluster *>(plNe->At(ipr)) ;
353 //Get the index where the cluster comes, to retrieve the corresponding vertex
355 if (reader->GetMixedEvent())
356 evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
359 //Do not count the candidate (photon or pi0) or the daughters of the candidate
360 if(calo->GetID() == pCandidate->GetCaloLabel(0) ||
361 calo->GetID() == pCandidate->GetCaloLabel(1) ) continue ;
363 //Skip matched clusters with tracks in case of neutral+charged analysis
364 if( fPartInCone == kNeutralAndCharged &&
365 pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
367 //Assume that come from vertex in straight line
368 calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
375 {// Mixed event stored in AliAODPWG4Particles
376 AliAODPWG4Particle * calomix = dynamic_cast<AliAODPWG4Particle*>(plNe->At(ipr)) ;
379 printf("AliIsolationCut::MakeIsolationCut() - Wrong calo data type, continue\n");
384 eta = calomix->Eta();
385 phi = calomix->Phi() ;
388 if( phi < 0 ) phi+=TMath::TwoPi();
391 // Only loop the particle at the same side of candidate
392 if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
394 // If at the same side has particle larger than candidate,
395 // then candidate can not be the leading, skip such events
403 pCandidate->SetLeadingParticle(kFALSE);
415 refclusters->Clear();
422 //Check if there is any particle inside cone with pt larger than fPtThreshold
424 rad = Radius(etaC, phiC, eta, phi);
427 printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
431 if(fDebug>0) printf(" - inside candidate cone");
438 refclusters = new TObjArray(0);
439 //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
440 TString tempo(aodArrayRefName) ;
441 tempo += "Clusters" ;
442 refclusters->SetName(tempo);
443 refclusters->SetOwner(kFALSE);
445 refclusters->Add(calo);
449 if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
450 //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
452 if( fPtFraction*ptC<fPtThreshold)
454 if(pt>fPtThreshold) nfrac++ ;
458 if(pt>fPtFraction*ptC) nfrac++;
462 if(pt>fPtFraction*ptC) nfrac++;
467 if(fDebug>0) printf("\n");
469 }// neutral particle loop
474 //Add reference arrays to AOD when filling AODs only
477 if(refclusters) pCandidate->AddObjArray(refclusters);
478 if(reftracks) pCandidate->AddObjArray(reftracks);
481 //Check isolation, depending on selected isolation criteria
482 if( fICMethod == kPtThresIC)
484 if(n==0) isolated = kTRUE ;
486 else if( fICMethod == kSumPtIC)
488 if(coneptsum < fSumPtThreshold)
491 else if( fICMethod == kPtFracIC)
493 if(nfrac==0) isolated = kTRUE ;
495 else if( fICMethod == kSumPtFracIC)
497 //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
499 if(fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold) isolated = kTRUE ;
500 if( fPtFraction*ptC > fSumPtThreshold && coneptsum < fPtFraction*ptC) isolated = kTRUE ;
504 if(coneptsum < fPtFraction*ptC) isolated = kTRUE ;
507 else if( fICMethod == kSumDensityIC)
509 // Get good cell density (number of active cells over all cells in cone)
510 // and correct energy in cone
511 Float_t cellDensity = GetCellDensity(pCandidate,reader);
512 if(coneptsum < fSumPtThreshold*cellDensity)
518 //_____________________________________________________
519 void AliIsolationCut::Print(const Option_t * opt) const
522 //Print some relevant parameters set for the analysis
526 printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
528 printf("IC method = %d\n", fICMethod ) ;
529 printf("Cone Size = %1.2f\n", fConeSize ) ;
530 printf("pT threshold = %2.1f\n", fPtThreshold) ;
531 printf("pT fraction = %3.1f\n", fPtFraction ) ;
532 printf("particle type in cone = %d\n", fPartInCone ) ;
533 printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
538 //___________________________________________________________________________
539 Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC,
540 const Float_t eta , const Float_t phi) const
542 // Calculate the distance to trigger from any particle
544 Float_t dEta = etaC-eta;
545 Float_t dPhi = phiC-phi;
547 if(TMath::Abs(dPhi) >= TMath::Pi())
548 dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
550 return TMath::Sqrt( dEta*dEta + dPhi*dPhi );