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