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