]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/CaloTrackCorrBase/AliIsolationCut.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
From Marta:
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliIsolationCut.cxx
CommitLineData
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
47ClassImp(AliIsolationCut)
48
c5693f62 49//____________________________________
50AliIsolationCut::AliIsolationCut() :
51TObject(),
52fConeSize(0.),
043c8ec2 53fPtThreshold(0.),
54fPtThresholdMax(10000.),
c5693f62 55fSumPtThreshold(0.),
56fPtFraction(0.),
57fICMethod(0),
03bae431 58fPartInCone(0),
22c7f802 59fDebug(-1),
60fFracIsThresh(1)
a5fb4114 61{
62 //default ctor
63
64 //Initialize parameters
65 InitParameters();
c5693f62 66
a5fb4114 67}
68
70561f53 69//__________________________________________________________________________________________________________________________________________
70void AliIsolationCut::CalculateUEBandClusterNormalization( AliCaloTrackReader * /*reader*/,const Float_t etaC, const Float_t /*phiC*/,
71 const Float_t phiUEptsumCluster, const Float_t etaUEptsumCluster,
72 Float_t & phiUEptsumClusterNorm, Float_t & etaUEptsumClusterNorm,
73 Float_t & excessFracEta, Float_t & excessFracPhi ) const
74{
75 // Normalize cluster background band
76
77 Float_t coneA = fConeSize*fConeSize*TMath::Pi(); // A = pi R^2, isolation cone area
78
79 //Careful here if EMCal limits changed .. 2010 (4 SM) to 2011-12 (10 SM), for the moment consider 100 deg in phi
80 Float_t emcEtaSize = 0.7*2; // TO FIX
81 Float_t emcPhiSize = TMath::DegToRad()*100.; // TO FIX
82
83 /* //Catherine code
84 if(((((2*fConeSize*emcPhiSize)-coneA))*phiBandBadCellsCoeff)!=0)phiUEptsumClusterNorm = phiUEptsumCluster*(coneA*coneBadCellsCoeff / (((2*fConeSize*emcPhiSize)-coneA))*phiBandBadCellsCoeff); // pi * R^2 / (2 R * 2 100 deg) - trigger cone
85 if(((((2*(fConeSize-excess)*emcPhiSize)-(coneA-excessFracEta))*etaBandBadCellsCoeff))!=0)phiUEptsumClusterNorm = phiUEptsumCluster*(coneA *coneBadCellsCoeff/ (((2*(fConeSize-excess)*emcPhiSize)-(coneA/excessFracEta))*etaBandBadCellsCoeff));
86 if(((2*(fConeSize-excess)*emcEtaSize)-(coneA-excessFracPhi))*phiBandBadCellsCoeff!=0) etaUEptsumClusterNorm = etaUEptsumCluster*(coneA*coneBadCellsCoeff / (((2*(fConeSize-excess)*emcEtaSize)-(coneA/excessFracPhi))*phiBandBadCellsCoeff));
87 */
88
89 if((2*fConeSize*emcPhiSize-coneA)!=0) phiUEptsumClusterNorm = phiUEptsumCluster*(coneA / (((2*fConeSize*emcPhiSize)-coneA))); // pi * R^2 / (2 R * 2 100 deg) - trigger cone
90 if((2*fConeSize*emcEtaSize-coneA)!=0) etaUEptsumClusterNorm = etaUEptsumCluster*(coneA / (((2*fConeSize*emcEtaSize)-coneA))); // pi * R^2 / (2 R * 2*0.7) - trigger cone
91
92 //out of eta acceptance
93 excessFracEta = 1;
94 excessFracPhi = 1;
f8d07abf 95
70561f53 96 if(TMath::Abs(etaC)+fConeSize > emcEtaSize/2.)
97 {
98 Float_t excess = TMath::Abs(etaC) + fConeSize - emcEtaSize/2.;
99 excessFracEta = CalculateExcessAreaFraction(excess);
100
101 if ( excessFracEta != 0) coneA /= excessFracEta;
102
103 //UE band is also out of acceptance, need to estimate corrected area
104 if(((2*fConeSize-excess)*emcPhiSize-coneA) != 0 ) phiUEptsumClusterNorm = phiUEptsumCluster*(coneA / ((((2*fConeSize-excess)*emcPhiSize)-coneA)));
105 if(( 2*fConeSize *emcEtaSize-coneA) != 0 ) etaUEptsumClusterNorm = etaUEptsumCluster*(coneA / ((( 2*fConeSize *emcEtaSize)-coneA)));
106 }
107
108}
109
110//________________________________________________________________________________________________________________________________________
111void AliIsolationCut::CalculateUEBandTrackNormalization ( AliCaloTrackReader * reader, const Float_t etaC, const Float_t /*phiC*/,
112 const Float_t phiUEptsumTrack, const Float_t etaUEptsumTrack,
113 Float_t & phiUEptsumTrackNorm, Float_t & etaUEptsumTrackNorm,
114 Float_t & excessFracEta, Float_t & excessFracPhi ) const
115{
116 // Normalize track background band
117
118 Float_t coneA = fConeSize*fConeSize*TMath::Pi(); // A = pi R^2, isolation cone area
119
120 // Get the cut used for the TPC tracks in the reader, +-0.8, +-0.9 ...
121 // Only valid in simple fidutial cut case and if the cut is applied, careful!
122 Float_t tpcEtaSize = reader->GetFiducialCut()->GetCTSFidCutMaxEtaArray()->At(0) -
123 reader->GetFiducialCut()->GetCTSFidCutMinEtaArray()->At(0) ;
124 Float_t tpcPhiSize = TMath::TwoPi();
125
126 /*//Catherine code
127 //phiUEptsumTrackNorm = phiUEptsumTrack*(coneA*coneBadCellsCoeff / (((2*fConeSize*tpcPhiSize)-coneA))*phiBandBadCellsCoeff); // pi * R^2 / (2 R * 2 pi) - trigger cone
128 //etaUEptsumTrackNorm = etaUEptsumTrack*(coneA*coneBadCellsCoeff / (((2*fConeSize*tpcEtaSize)-coneA))*etaBandBadCellsCoeff); // pi * R^2 / (2 R * 1.6) - trigger cone
129 if((2*fConeSize*tpcPhiSize-coneA)!=0)phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / (((2*fConeSize*tpcPhiSize)-coneA))); // pi * R^2 / (2 R * 2 pi) - trigger cone
130 if((2*fConeSize*tpcEtaSize-coneA)!=0)etaUEptsumTrackNorm = etaUEptsumTrack*(coneA / (((2*fConeSize*tpcEtaSize)-coneA))); // pi * R^2 / (2 R * 1.6) - trigger cone
131 if((2*(fConeSize-excess)*tpcPhiSize)-(coneA-excessFracEta)!=0)phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / (((2*(fConeSize-excess)*tpcPhiSize)-(coneA/excessFracEta))));
132 */ //end Catherine code
133
134 //correct out of eta acceptance
135 excessFracEta = 1;
136 excessFracPhi = 1;
137
138 if((2*fConeSize*tpcPhiSize-coneA)!=0) phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / (((2*fConeSize*tpcPhiSize)-coneA))); // pi * R^2 / (2 R * 2 pi) - trigger cone
139 if((2*fConeSize*tpcEtaSize-coneA)!=0) etaUEptsumTrackNorm = etaUEptsumTrack*(coneA / (((2*fConeSize*tpcEtaSize)-coneA))); // pi * R^2 / (2 R * 1.6) - trigger cone
140
141 if(TMath::Abs(etaC)+fConeSize > tpcEtaSize/2.)
142 {
143 Float_t excess = TMath::Abs(etaC) + fConeSize - tpcEtaSize/2.;
144 excessFracEta = CalculateExcessAreaFraction(excess);
145 if (excessFracEta != 0) coneA /= excessFracEta;
146
147 //UE band is also out of acceptance, need to estimate corrected area
148 if(((2*fConeSize-excess)*tpcPhiSize - coneA) !=0 ) phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / ((((2*fConeSize-excess)*tpcPhiSize)-coneA)));
149 if(( 2*fConeSize *tpcEtaSize - coneA) !=0 ) etaUEptsumTrackNorm = etaUEptsumTrack*(coneA / ((( 2*fConeSize *tpcEtaSize)-coneA)));
150 }
151
152}
153
154//______________________________________________________________________________
155Float_t AliIsolationCut::CalculateExcessAreaFraction(const Float_t excess) const
156{
157 // Area of a circunference segment segment 1/2 R^2 (angle-sin(angle)), angle = 2*ACos((R-excess)/R)
158
159
160 Float_t angle = 2*TMath::ACos( (fConeSize-excess) / fConeSize );
161
162 Float_t coneA = fConeSize*fConeSize*TMath::Pi(); // A = pi R^2, isolation cone area
163
164 Float_t excessA = fConeSize*fConeSize / 2 * (angle-TMath::Sin(angle));
165
166 if(coneA > excessA) return coneA / (coneA-excessA);
167 else
168 {
169 printf("AliIsolationCut::CalculateExcessAreaFraction() - Please Check : Excess Track %2.3f, coneA %2.2f, excessA %2.2f, angle %2.2f,factor %2.2f\n",
170 excess,coneA, excessA, angle*TMath::RadToDeg(), coneA / (coneA-excessA));
171 return 1;
172 }
173}
174
175//_______________________________________________________________________________________
176Float_t AliIsolationCut::GetCellDensity(AliAODPWG4ParticleCorrelation * pCandidate,
177 AliCaloTrackReader * reader) const
f8d07abf 178{
179 // Get good cell density (number of active cells over all cells in cone)
180
181 Double_t coneCells = 0.; //number of cells in cone with radius fConeSize
182 Double_t coneCellsBad = 0.; //number of bad cells in cone with radius fConeSize
183 Double_t cellDensity = 1.;
184
185 Float_t phiC = pCandidate->Phi() ;
186 if(phiC<0) phiC+=TMath::TwoPi();
187 Float_t etaC = pCandidate->Eta() ;
188
189 if(pCandidate->GetDetector()=="EMCAL")
190 {
191 AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
192 AliCalorimeterUtils *cu = reader->GetCaloUtils();
193
194 Int_t absId = -999;
195 if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
196 {
197 //Get absolute (col,row) of candidate
198 Int_t iEta=-1, iPhi=-1, iRCU = -1;
199 Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
200
201 Int_t colC = iEta;
202 if (nSupMod % 2) colC = AliEMCALGeoParams::fgkEMCALCols + iEta ;
203 Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
204
205 Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
206 //loop on cells in a square of side fConeSize to check cells in cone
207 for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
208 {
209 for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
210 {
211 if (Radius(colC, rowC, icol, irow) < sqrSize)
212 {
213 coneCells += 1.;
214
215 Int_t cellSM = -999;
216 Int_t cellEta = -999;
217 Int_t cellPhi = -999;
218 if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
219 {
220 cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
221 cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
222 cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
223 }
224 if(icol < AliEMCALGeoParams::fgkEMCALCols)
225 {
226 cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
227 cellEta = icol;
228 cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
229 }
230
231 //Count as bad "cells" out of EMCAL acceptance
232 if(icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2 ||
233 irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*16./3) //5*nRows+1/3*nRows
234 {
235 coneCellsBad += 1.;
236 }
237 //Count as bad "cells" marked as bad in the DataBase
238 else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1)
239 {
240 coneCellsBad += 1. ;
241 }
242 }
243 }
244 }//end of cells loop
245 }
246
247 else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
248
249 if (coneCells > 0.)
250 {
251 cellDensity = (coneCells-coneCellsBad)/coneCells;
8c00ddb0 252 //printf("Energy density = %f\n", cellDensity);
f8d07abf 253 }
254 }
255
256 return cellDensity;
257
258}
259
70561f53 260//__________________________________________________________________________________
261void AliIsolationCut::GetCoeffNormBadCell(AliAODPWG4ParticleCorrelation * pCandidate,
262 AliCaloTrackReader * reader,
263 Float_t & coneBadCellsCoeff,
264 Float_t & etaBandBadCellsCoeff,
265 Float_t & phiBandBadCellsCoeff)
266{
267 // Get good cell density (number of active cells over all cells in cone)
268
269 Double_t coneCells = 0.; //number of cells in cone with radius fConeSize
270 Double_t phiBandCells = 0.; //number of cells in band phi
271 Double_t etaBandCells = 0.; //number of cells in band eta
272
273 Float_t phiC = pCandidate->Phi() ;
274 if(phiC<0) phiC+=TMath::TwoPi();
275 Float_t etaC = pCandidate->Eta() ;
276
277 if(pCandidate->GetDetector()=="EMCAL")
278 {
279 AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
280 AliCalorimeterUtils *cu = reader->GetCaloUtils();
281
282 Int_t absId = -999;
283 if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
284 {
285 //Get absolute (col,row) of candidate
286 Int_t iEta=-1, iPhi=-1, iRCU = -1;
287 Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(),
288 iEta, iPhi, iRCU);
289
290 Int_t colC = iEta;
291 if (nSupMod % 2) colC = AliEMCALGeoParams::fgkEMCALCols + iEta ;
292 Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
293
294 Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
295 for(Int_t icol = 0; icol < 2*AliEMCALGeoParams::fgkEMCALCols-1;icol++)
296 {
297 for(Int_t irow = 0; irow < 5*AliEMCALGeoParams::fgkEMCALRows -1; irow++)
298 {
299 //loop on cells in a square of side fConeSize to check cells in cone
300 if ( Radius(colC, rowC, icol, irow) < sqrSize ) { coneCells += 1.; }
301 else if( icol>colC-sqrSize && icol<colC+sqrSize ) { phiBandCells += 1 ; }
302 else if( irow>rowC-sqrSize && irow<rowC+sqrSize ) { etaBandCells += 1 ; }
303
304 Int_t cellSM = -999;
305 Int_t cellEta = -999;
306 Int_t cellPhi = -999;
307 if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
308 {
309 cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
310 cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
311 cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
312 }
313 if(icol < AliEMCALGeoParams::fgkEMCALCols)
314 {
315 cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
316 cellEta = icol;
317 cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
318 }
319
320 if( (icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2-1 ||
321 irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*5 - 1) //5*nRows+1/3*nRows //Count as bad "cells" out of EMCAL acceptance
322 || (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1)) //Count as bad "cells" marked as bad in the DataBase
323 {
324 if ( Radius(colC, rowC, icol, irow) < sqrSize ) coneBadCellsCoeff += 1.;
325 else if( icol>colC-sqrSize && icol<colC+sqrSize ) phiBandBadCellsCoeff += 1 ;
326 else if( irow>rowC-sqrSize && irow<rowC+sqrSize ) etaBandBadCellsCoeff += 1 ;
327 }
328 }
329 }//end of cells loop
330 }
331
332 else if(fDebug > 0) printf("cluster with bad (eta,phi) in EMCal for energy density coeff calculation\n");
333
334 if (coneCells > 0.)
335 {
336 // printf("Energy density coneBadCellsCoeff= %.2f coneCells%.2f\n", coneBadCellsCoeff,coneCells);
337 coneBadCellsCoeff = (coneCells-coneBadCellsCoeff)/coneCells;
338 // printf("coneBadCellsCoeff= %.2f\n", coneBadCellsCoeff);
339 }
340 if (phiBandCells > 0.)
341 {
342 // printf("Energy density phiBandBadCellsCoeff = %.2f phiBandCells%.2f\n", phiBandBadCellsCoeff,phiBandCells);
343 phiBandBadCellsCoeff = (phiBandCells-phiBandBadCellsCoeff)/phiBandCells;
344 // printf("phiBandBadCellsCoeff = %.2f\n", phiBandBadCellsCoeff);
345 }
346 if (etaBandCells > 0.)
347 {
348 //printf("Energy density etaBandBadCellsCoeff = %.2f etaBandCells%.2f\n", etaBandBadCellsCoeff,etaBandCells);
349 etaBandBadCellsCoeff = (etaBandCells-etaBandBadCellsCoeff)/etaBandCells;
350 // printf("etaBandBadCellsCoeff = %.2f\n",etaBandBadCellsCoeff);
351 }
352
353 }
354
355}
356
c5693f62 357//____________________________________________
a5fb4114 358TString AliIsolationCut::GetICParametersList()
359{
360 //Put data member values in string to keep in output container
361
362 TString parList ; //this will be list of parameters used for this analysis.
363 const Int_t buffersize = 255;
364 char onePar[buffersize] ;
365
366 snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
367 parList+=onePar ;
368 snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
369 parList+=onePar ;
370 snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
371 parList+=onePar ;
372 snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
373 parList+=onePar ;
374 snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
375 parList+=onePar ;
376 snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
377 parList+=onePar ;
22c7f802 378 snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
379 parList+=onePar ;
380
a5fb4114 381 return parList;
382}
383
f18cb329 384//____________________________________
a5fb4114 385void AliIsolationCut::InitParameters()
386{
387 //Initialize the parameters of the analysis.
388
389 fConeSize = 0.4 ;
043c8ec2 390 fPtThreshold = 1. ;
391 fPtThresholdMax = 10000. ;
a5fb4114 392 fSumPtThreshold = 0.5 ;
393 fPtFraction = 0.1 ;
394 fPartInCone = kOnlyCharged;
395 fICMethod = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
22c7f802 396 fFracIsThresh = 1;
a5fb4114 397}
398
c5693f62 399//________________________________________________________________________________
70561f53 400void AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,
401 TObjArray * plNe,
402 AliCaloTrackReader * reader,
403 AliCaloPID * pid,
c5693f62 404 const Bool_t bFillAOD,
405 AliAODPWG4ParticleCorrelation *pCandidate,
406 const TString & aodArrayRefName,
b960c7eb 407 Int_t & n,
408 Int_t & nfrac,
409 Float_t & coneptsum,
410 Bool_t & isolated) const
70561f53 411{
a5fb4114 412 //Search in cone around a candidate particle if it is isolated
03bae431 413 Float_t ptC = pCandidate->Pt() ;
a5fb4114 414 Float_t phiC = pCandidate->Phi() ;
415 if(phiC<0) phiC+=TMath::TwoPi();
416 Float_t etaC = pCandidate->Eta() ;
70561f53 417
a5fb4114 418 Float_t pt = -100. ;
419 Float_t eta = -100. ;
420 Float_t phi = -100. ;
421 Float_t rad = -100. ;
422
70561f53 423 Float_t coneptsumCluster = 0;
424 Float_t coneptsumTrack = 0;
425
426 Float_t etaBandPtSumTrack = 0;
427 Float_t phiBandPtSumTrack = 0;
428 Float_t etaBandPtSumCluster = 0;
429 Float_t phiBandPtSumCluster = 0;
430
a5fb4114 431 n = 0 ;
432 nfrac = 0 ;
a5fb4114 433 isolated = kFALSE;
70561f53 434
03bae431 435 if(fDebug>0)
436 {
437 printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
438 pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
439 if(plCTS) printf(", nTracks %d" ,plCTS->GetEntriesFast());
440 if(plNe) printf(", nClusters %d",plNe ->GetEntriesFast());
441
442 printf("\n");
443 }
1a31a9ab 444
a5fb4114 445 //Initialize the array with refrences
a258315d 446 TObjArray * refclusters = 0x0;
447 TObjArray * reftracks = 0x0;
448 Int_t ntrackrefs = 0;
449 Int_t nclusterrefs = 0;
ac5111f9 450
a5fb4114 451 //Check charged particles in cone.
b960c7eb 452 if(plCTS &&
453 (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged))
454 {
a5fb4114 455 TVector3 p3;
b960c7eb 456 for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
457 {
a258315d 458 AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ;
b960c7eb 459
a258315d 460 if(track)
461 {
462 //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
463 if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
464 track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3) ) continue ;
465
466 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
467 pt = p3.Pt();
468 eta = p3.Eta();
469 phi = p3.Phi() ;
470 }
471 else
472 {// Mixed event stored in AliAODPWG4Particles
473 AliAODPWG4Particle * trackmix = dynamic_cast<AliAODPWG4Particle*>(plCTS->At(ipr)) ;
474 if(!trackmix)
475 {
476 printf("AliIsolationCut::MakeIsolationCut() - Wrong track data type, continue\n");
477 continue;
478 }
479
480 pt = trackmix->Pt();
481 eta = trackmix->Eta();
482 phi = trackmix->Phi() ;
483 }
b960c7eb 484
a258315d 485 if( phi < 0 ) phi+=TMath::TwoPi();
a5fb4114 486
70561f53 487 rad = Radius(etaC, phiC, eta, phi);
488
489 // ** For the background out of cone **
490
491 if(rad > fConeSize)
492 {
493 if(eta > (etaC-fConeSize) && eta < (etaC+fConeSize)) phiBandPtSumTrack += pt;
494 if(phi > (phiC-fConeSize) && phi < (phiC+fConeSize)) etaBandPtSumTrack += pt;
495 }
496
497 // ** For the isolated particle **
498
b960c7eb 499 // Only loop the particle at the same side of candidate
70561f53 500 if(TMath::Abs(phi-phiC) > TMath::PiOver2()) continue ;
501
b960c7eb 502 // If at the same side has particle larger than candidate,
503 // then candidate can not be the leading, skip such events
504 if(pt > ptC)
505 {
a5fb4114 506 n = -1;
507 nfrac = -1;
70561f53 508 coneptsumTrack = -1;
a5fb4114 509 isolated = kFALSE;
70561f53 510
3f150b4b 511 pCandidate->SetLeadingParticle(kFALSE);
512
b960c7eb 513 if(bFillAOD && reftracks)
514 {
a5fb4114 515 reftracks->Clear();
516 delete reftracks;
517 }
b960c7eb 518
a5fb4114 519 return ;
520 }
f18cb329 521
a5fb4114 522 //Check if there is any particle inside cone with pt larger than fPtThreshold
a5fb4114 523
70561f53 524 if( fDebug > 0 )
03bae431 525 printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
70561f53 526
b960c7eb 527 if(rad < fConeSize)
528 {
70561f53 529 if(fDebug > 0) printf(" - inside candidate cone");
530
b960c7eb 531 if(bFillAOD)
532 {
a5fb4114 533 ntrackrefs++;
b960c7eb 534 if(ntrackrefs == 1)
535 {
a5fb4114 536 reftracks = new TObjArray(0);
537 //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
538 TString tempo(aodArrayRefName) ;
539 tempo += "Tracks" ;
540 reftracks->SetName(tempo);
541 reftracks->SetOwner(kFALSE);
542 }
543 reftracks->Add(track);
544 }
b960c7eb 545
70561f53 546 coneptsumTrack+=pt;
547 if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
a5fb4114 548 if(pt > fPtFraction*ptC ) nfrac++;
b960c7eb 549
a5fb4114 550 } // Inside cone
70561f53 551
03bae431 552 if(fDebug>0) printf("\n");
70561f53 553
a5fb4114 554 }// charged particle loop
03bae431 555
556
a5fb4114 557 }//Tracks
558
70561f53 559
a5fb4114 560 //Check neutral particles in cone.
70561f53 561 if(plNe &&
b960c7eb 562 (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged))
563 {
a5fb4114 564 TLorentzVector mom ;
b960c7eb 565
566 for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
567 {
a258315d 568 AliVCluster * calo = dynamic_cast<AliVCluster *>(plNe->At(ipr)) ;
b960c7eb 569
a258315d 570 if(calo)
571 {
572 //Get the index where the cluster comes, to retrieve the corresponding vertex
573 Int_t evtIndex = 0 ;
574 if (reader->GetMixedEvent())
575 evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
576
577
578 //Do not count the candidate (photon or pi0) or the daughters of the candidate
579 if(calo->GetID() == pCandidate->GetCaloLabel(0) ||
580 calo->GetID() == pCandidate->GetCaloLabel(1) ) continue ;
581
582 //Skip matched clusters with tracks in case of neutral+charged analysis
583 if( fPartInCone == kNeutralAndCharged &&
584 pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
585
586 //Assume that come from vertex in straight line
587 calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
588
589 pt = mom.Pt() ;
590 eta = mom.Eta() ;
591 phi = mom.Phi() ;
592 }
593 else
594 {// Mixed event stored in AliAODPWG4Particles
595 AliAODPWG4Particle * calomix = dynamic_cast<AliAODPWG4Particle*>(plNe->At(ipr)) ;
596 if(!calomix)
597 {
598 printf("AliIsolationCut::MakeIsolationCut() - Wrong calo data type, continue\n");
599 continue;
600 }
601
602 pt = calomix->Pt();
603 eta = calomix->Eta();
604 phi = calomix->Phi() ;
605 }
a5fb4114 606
a258315d 607 if( phi < 0 ) phi+=TMath::TwoPi();
a5fb4114 608
70561f53 609 rad = Radius(etaC, phiC, eta, phi);
610
611 // ** For the background out of cone **
612
613 if(rad > fConeSize)
614 {
615 if(eta > (etaC-fConeSize) && eta < (etaC+fConeSize)) phiBandPtSumCluster += pt;
616 if(phi > (phiC-fConeSize) && phi < (phiC+fConeSize)) etaBandPtSumCluster += pt;
617 }
618
619 // ** For the isolated particle **
a5fb4114 620
b960c7eb 621 // Only loop the particle at the same side of candidate
ac5111f9 622 if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
623
b960c7eb 624 // If at the same side has particle larger than candidate,
625 // then candidate can not be the leading, skip such events
626 if(pt > ptC)
627 {
a5fb4114 628 n = -1;
629 nfrac = -1;
70561f53 630 coneptsumCluster = -1;
a5fb4114 631 isolated = kFALSE;
b960c7eb 632
3f150b4b 633 pCandidate->SetLeadingParticle(kFALSE);
634
b960c7eb 635 if(bFillAOD)
636 {
637 if(reftracks)
638 {
a5fb4114 639 reftracks ->Clear();
640 delete reftracks;
641 }
b960c7eb 642
643 if(refclusters)
644 {
a5fb4114 645 refclusters->Clear();
646 delete refclusters;
647 }
648 }
649 return ;
650 }
651
652 //Check if there is any particle inside cone with pt larger than fPtThreshold
f18cb329 653
70561f53 654 if(fDebug > 0 )
03bae431 655 printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
656
b960c7eb 657 if(rad < fConeSize)
658 {
70561f53 659 if(fDebug > 0 ) printf(" - inside candidate cone");
660
b960c7eb 661 if(bFillAOD)
662 {
a5fb4114 663 nclusterrefs++;
b960c7eb 664 if(nclusterrefs==1)
665 {
a5fb4114 666 refclusters = new TObjArray(0);
667 //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
668 TString tempo(aodArrayRefName) ;
669 tempo += "Clusters" ;
670 refclusters->SetName(tempo);
671 refclusters->SetOwner(kFALSE);
672 }
673 refclusters->Add(calo);
674 }
b960c7eb 675
70561f53 676 coneptsumCluster+=pt;
677 if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
a5fb4114 678 //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
70561f53 679 if(fFracIsThresh)
680 {
681 if( fPtFraction*ptC<fPtThreshold)
682 {
683 if(pt>fPtThreshold) nfrac++ ;
684 }
685 else
686 {
687 if(pt>fPtFraction*ptC) nfrac++;
688 }
689 }
690 else
691 {
692 if(pt>fPtFraction*ptC) nfrac++;
693 }
b960c7eb 694
a5fb4114 695 }//in cone
b960c7eb 696
03bae431 697 if(fDebug>0) printf("\n");
70561f53 698
a5fb4114 699 }// neutral particle loop
70561f53 700
a5fb4114 701 }//neutrals
1a31a9ab 702
a5fb4114 703
704 //Add reference arrays to AOD when filling AODs only
b960c7eb 705 if(bFillAOD)
706 {
a5fb4114 707 if(refclusters) pCandidate->AddObjArray(refclusters);
708 if(reftracks) pCandidate->AddObjArray(reftracks);
709 }
b960c7eb 710
70561f53 711 coneptsum = coneptsumCluster+coneptsumTrack;
712
b960c7eb 713 //Check isolation, depending on selected isolation criteria
714 if( fICMethod == kPtThresIC)
715 {
a5fb4114 716 if(n==0) isolated = kTRUE ;
717 }
b960c7eb 718 else if( fICMethod == kSumPtIC)
719 {
a5fb4114 720 if(coneptsum < fSumPtThreshold)
721 isolated = kTRUE ;
722 }
b960c7eb 723 else if( fICMethod == kPtFracIC)
724 {
a5fb4114 725 if(nfrac==0) isolated = kTRUE ;
726 }
b960c7eb 727 else if( fICMethod == kSumPtFracIC)
728 {
a5fb4114 729 //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
70561f53 730 // printf("photon analysis IsDataMC() ?%i\n",IsDataMC());
731 if(fFracIsThresh )
732 {
733 if( fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold) isolated = kTRUE ;
22c7f802 734 if( fPtFraction*ptC > fSumPtThreshold && coneptsum < fPtFraction*ptC) isolated = kTRUE ;
735 }
736 else
70561f53 737 {
738 if(coneptsum < fPtFraction*ptC) isolated = kTRUE ;
739 }
a5fb4114 740 }
70561f53 741 else if( fICMethod == kSumDensityIC)
f8d07abf 742 {
743 // Get good cell density (number of active cells over all cells in cone)
744 // and correct energy in cone
70561f53 745
f8d07abf 746 Float_t cellDensity = GetCellDensity(pCandidate,reader);
70561f53 747
f8d07abf 748 if(coneptsum < fSumPtThreshold*cellDensity)
749 isolated = kTRUE;
750 }
70561f53 751 else if( fICMethod == kSumBkgSubIC)
752 {
753 Double_t coneptsumBkg = 0.;
754 Float_t etaBandPtSumTrackNorm = 0;
755 Float_t phiBandPtSumTrackNorm = 0;
756 Float_t etaBandPtSumClusterNorm = 0;
757 Float_t phiBandPtSumClusterNorm = 0;
758
759 Float_t excessFracEtaTrack = 1;
760 Float_t excessFracPhiTrack = 1;
761 Float_t excessFracEtaCluster = 1;
762 Float_t excessFracPhiCluster = 1;
763
764 // Normalize background to cone area
765 if (fPartInCone != kOnlyCharged )
766 CalculateUEBandClusterNormalization(reader, etaC, phiC,
767 phiBandPtSumCluster , etaBandPtSumCluster,
768 phiBandPtSumClusterNorm, etaBandPtSumClusterNorm,
769 excessFracEtaCluster , excessFracPhiCluster );
770 if (fPartInCone != kOnlyNeutral )
771 CalculateUEBandTrackNormalization(reader, etaC, phiC,
772 phiBandPtSumTrack , etaBandPtSumTrack ,
773 phiBandPtSumTrackNorm, etaBandPtSumTrackNorm,
774 excessFracEtaTrack , excessFracPhiTrack );
775
776 if (fPartInCone == kOnlyCharged ) coneptsumBkg = etaBandPtSumTrackNorm;
777 else if(fPartInCone == kOnlyNeutral ) coneptsumBkg = etaBandPtSumClusterNorm;
778 else if(fPartInCone == kNeutralAndCharged ) coneptsumBkg = etaBandPtSumClusterNorm + etaBandPtSumTrackNorm;
779
780 //coneptsumCluster*=(coneBadCellsCoeff*excessFracEtaCluster*excessFracPhiCluster) ; // apply this correction earlier???
781 // line commented out in last modif!!!
782
783 coneptsum = coneptsumCluster+coneptsumTrack;
784
785 coneptsum -= coneptsumBkg;
786 if(coneptsum < fSumPtThreshold)
787 isolated = kTRUE ;
788 }
a5fb4114 789
790}
791
c5693f62 792//_____________________________________________________
a5fb4114 793void AliIsolationCut::Print(const Option_t * opt) const
794{
795
796 //Print some relevant parameters set for the analysis
797 if(! opt)
798 return;
799
800 printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
801
c5693f62 802 printf("IC method = %d\n", fICMethod ) ;
803 printf("Cone Size = %1.2f\n", fConeSize ) ;
a5fb4114 804 printf("pT threshold = %2.1f\n", fPtThreshold) ;
c5693f62 805 printf("pT fraction = %3.1f\n", fPtFraction ) ;
806 printf("particle type in cone = %d\n", fPartInCone ) ;
22c7f802 807 printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
a5fb4114 808 printf(" \n") ;
809
810}
f18cb329 811
812//___________________________________________________________________________
813Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC,
814 const Float_t eta , const Float_t phi) const
815{
816 // Calculate the distance to trigger from any particle
817
818 Float_t dEta = etaC-eta;
819 Float_t dPhi = phiC-phi;
820
821 if(TMath::Abs(dPhi) >= TMath::Pi())
822 dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
823
824 return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
825
826}
827
828
829