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