AliPhysics  41af4b0 (41af4b0)
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  Int_t status = 0;
201  // Loop on cells in a square of side fConeSize to check cells in cone
202  for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
203  {
204  for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
205  {
206  if (Radius(colC, rowC, icol, irow) < sqrSize)
207  {
208  coneCells += 1.;
209 
210  Int_t cellSM = -999;
211  Int_t cellEta = -999;
212  Int_t cellPhi = -999;
213  if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
214  {
215  cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
216  cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
217  cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
218  }
219  if(icol < AliEMCALGeoParams::fgkEMCALCols)
220  {
221  cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
222  cellEta = icol;
223  cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
224  }
225 
226  //Count as bad "cells" out of EMCAL acceptance
227  if(icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2 ||
228  irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*16./3) //5*nRows+1/3*nRows
229  {
230  coneCellsBad += 1.;
231  }
232  // Count as bad "cells" marked as bad in the DataBase
233  else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi,status)==1)
234  {
235  coneCellsBad += 1. ;
236  }
237  }
238  }
239  }//end of cells loop
240  }
241  else AliWarning("Cluster with bad (eta,phi) in EMCal for energy density calculation");
242 
243  if (coneCells > 0.)
244  {
245  cellDensity = (coneCells-coneCellsBad)/coneCells;
246  //printf("Energy density = %f\n", cellDensity);
247  }
248  }
249 
250  return cellDensity;
251 }
252 
253 //___________________________________________________________________________________
255 //___________________________________________________________________________________
257  AliCaloTrackReader * reader,
258  Float_t & coneBadCellsCoeff,
259  Float_t & etaBandBadCellsCoeff,
260  Float_t & phiBandBadCellsCoeff)
261 {
262  Double_t coneCells = 0.; //number of cells in cone with radius fConeSize
263  Double_t phiBandCells = 0.; //number of cells in band phi
264  Double_t etaBandCells = 0.; //number of cells in band eta
265 
266  Float_t phiC = pCandidate->Phi() ;
267  if ( phiC < 0 ) phiC+=TMath::TwoPi();
268  Float_t etaC = pCandidate->Eta() ;
269 
270  if(pCandidate->GetDetectorTag() == AliCaloTrackReader::kEMCAL)
271  {
272  AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
273  AliCalorimeterUtils *cu = reader->GetCaloUtils();
274 
275  Int_t absId = -999;
276  if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
277  {
278  //Get absolute (col,row) of candidate
279  Int_t iEta=-1, iPhi=-1, iRCU = -1;
280  Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetectorTag(),
281  iEta, iPhi, iRCU);
282 
283  Int_t colC = iEta;
284  if (nSupMod % 2) colC = AliEMCALGeoParams::fgkEMCALCols + iEta ;
285  Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
286 
287  Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
288  Int_t status = 0;
289  for(Int_t icol = 0; icol < 2*AliEMCALGeoParams::fgkEMCALCols-1;icol++)
290  {
291  for(Int_t irow = 0; irow < 5*AliEMCALGeoParams::fgkEMCALRows -1; irow++)
292  {
293  //loop on cells in a square of side fConeSize to check cells in cone
294  if ( Radius(colC, rowC, icol, irow) < sqrSize ) { coneCells += 1.; }
295  else if( icol>colC-sqrSize && icol<colC+sqrSize ) { phiBandCells += 1 ; }
296  else if( irow>rowC-sqrSize && irow<rowC+sqrSize ) { etaBandCells += 1 ; }
297 
298  Int_t cellSM = -999;
299  Int_t cellEta = -999;
300  Int_t cellPhi = -999;
301  if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
302  {
303  cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
304  cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
305  cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
306  }
307  if(icol < AliEMCALGeoParams::fgkEMCALCols)
308  {
309  cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
310  cellEta = icol;
311  cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
312  }
313 
314  if( (icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2-1 ||
315  irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*5 - 1) //5*nRows+1/3*nRows //Count as bad "cells" out of EMCAL acceptance
316  || (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi,status)==1)) //Count as bad "cells" marked as bad in the DataBase
317  {
318  if ( Radius(colC, rowC, icol, irow) < sqrSize ) coneBadCellsCoeff += 1.;
319  else if( icol>colC-sqrSize && icol<colC+sqrSize ) phiBandBadCellsCoeff += 1 ;
320  else if( irow>rowC-sqrSize && irow<rowC+sqrSize ) etaBandBadCellsCoeff += 1 ;
321  }
322  }
323  }//end of cells loop
324  }
325  else AliWarning("Cluster with bad (eta,phi) in EMCal for energy density coeff calculation");
326 
327  if (coneCells > 0.)
328  {
329  // printf("Energy density coneBadCellsCoeff= %.2f coneCells%.2f\n", coneBadCellsCoeff,coneCells);
330  coneBadCellsCoeff = (coneCells-coneBadCellsCoeff)/coneCells;
331  // printf("coneBadCellsCoeff= %.2f\n", coneBadCellsCoeff);
332  }
333 
334  if (phiBandCells > 0.)
335  {
336  // printf("Energy density phiBandBadCellsCoeff = %.2f phiBandCells%.2f\n", phiBandBadCellsCoeff,phiBandCells);
337  phiBandBadCellsCoeff = (phiBandCells-phiBandBadCellsCoeff)/phiBandCells;
338  // printf("phiBandBadCellsCoeff = %.2f\n", phiBandBadCellsCoeff);
339  }
340 
341  if (etaBandCells > 0.)
342  {
343  //printf("Energy density etaBandBadCellsCoeff = %.2f etaBandCells%.2f\n", etaBandBadCellsCoeff,etaBandCells);
344  etaBandBadCellsCoeff = (etaBandCells-etaBandBadCellsCoeff)/etaBandCells;
345  // printf("etaBandBadCellsCoeff = %.2f\n",etaBandBadCellsCoeff);
346  }
347  }
348 }
349 
350 //____________________________________________
351 // Put data member values in string to keep
352 // in output container.
353 //____________________________________________
355 {
356  TString parList ; //this will be list of parameters used for this analysis.
357  const Int_t buffersize = 255;
358  char onePar[buffersize] ;
359 
360  snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
361  parList+=onePar ;
362  snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
363  parList+=onePar ;
364  snprintf(onePar,buffersize,"fPtThreshold >%2.2f;<%2.2f (isolation pt threshold) \n",fPtThreshold,fPtThresholdMax) ;
365  parList+=onePar ;
366  snprintf(onePar,buffersize,"fSumPtThreshold >%2.2f;<%2.2f (isolation sum pt threshold) \n",fSumPtThreshold,fSumPtThresholdMax) ;
367  parList+=onePar ;
368  snprintf(onePar,buffersize,"fPtFraction=%2.2f (isolation pt threshold fraction) \n",fPtFraction) ;
369  parList+=onePar ;
370  snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
371  parList+=onePar ;
372  snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
373  parList+=onePar ;
374  snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
375  parList+=onePar ;
376  snprintf(onePar,buffersize,"fDistMinToTrigger=%1.2f \n",fDistMinToTrigger) ;
377  parList+=onePar ;
378 
379  return parList;
380 }
381 
382 //____________________________________
383 // Initialize the parameters of the analysis.
384 //____________________________________
386 {
387  fConeSize = 0.4 ;
388  fPtThreshold = 0.5 ;
389  fPtThresholdMax = 10000. ;
390  fSumPtThreshold = 1.0 ;
391  fSumPtThresholdMax = 10000. ;
392  fPtFraction = 0.1 ;
394  fICMethod = kSumPtIC; // 0 pt threshol method, 1 cone pt sum method
395  fFracIsThresh = 1;
396  fDistMinToTrigger = -1.; // no effect
397 }
398 
399 //________________________________________________________________________________
416 //________________________________________________________________________________
418  TObjArray * plNe,
419  AliCaloTrackReader * reader,
420  AliCaloPID * pid,
421  Bool_t bFillAOD,
423  TString aodArrayRefName,
424  Int_t & n,
425  Int_t & nfrac,
426  Float_t & coneptsum, Float_t & ptLead,
427  Bool_t & isolated)
428 {
429  Float_t ptC = pCandidate->Pt() ;
430  Float_t phiC = pCandidate->Phi() ;
431  if ( phiC < 0 ) phiC+=TMath::TwoPi();
432  Float_t etaC = pCandidate->Eta() ;
433 
434  Float_t pt = -100. ;
435  Float_t eta = -100. ;
436  Float_t phi = -100. ;
437  Float_t rad = -100. ;
438 
439  Float_t coneptsumCluster = 0;
440  Float_t coneptsumTrack = 0;
441 
442  Float_t etaBandPtSumTrack = 0;
443  Float_t phiBandPtSumTrack = 0;
444  Float_t etaBandPtSumCluster = 0;
445  Float_t phiBandPtSumCluster = 0;
446 
447  n = 0 ;
448  nfrac = 0 ;
449  isolated = kFALSE;
450 
451  AliDebug(1,Form("Candidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
452  pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD));
453 
454  //Initialize the array with refrences
455  TObjArray * refclusters = 0x0;
456  TObjArray * reftracks = 0x0;
457  Int_t ntrackrefs = 0;
458  Int_t nclusterrefs = 0;
459 
460  // --------------------------------
461  // Check charged tracks in cone.
462  // --------------------------------
463 
464  if(plCTS &&
466  {
467  for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
468  {
469  AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ;
470 
471  if(track)
472  {
473  // In case of isolation of single tracks or conversion photon (2 tracks) or pi0 (4 tracks),
474  // do not count the candidate or the daughters of the candidate
475  // in the isolation conte
476  if ( pCandidate->GetDetectorTag() == AliFiducialCut::kCTS ) // make sure conversions are tagged as kCTS!!!
477  {
478  Int_t trackID = reader->GetTrackID(track) ; // needed instead of track->GetID() since AOD needs some manipulations
479  Bool_t contained = kFALSE;
480 
481  for(Int_t i = 0; i < 4; i++)
482  {
483  if( trackID == pCandidate->GetTrackLabel(i) ) contained = kTRUE;
484  }
485 
486  if ( contained ) continue ;
487  }
488 
489  fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
490  pt = fTrackVector.Pt();
491  eta = fTrackVector.Eta();
492  phi = fTrackVector.Phi() ;
493  }
494  else
495  {// Mixed event stored in AliCaloTrackParticles
496  AliCaloTrackParticle * trackmix = dynamic_cast<AliCaloTrackParticle*>(plCTS->At(ipr)) ;
497  if(!trackmix)
498  {
499  AliWarning("Wrong track data type, continue");
500  continue;
501  }
502 
503  pt = trackmix->Pt();
504  eta = trackmix->Eta();
505  phi = trackmix->Phi() ;
506  }
507 
508  // ** Calculate distance between candidate and tracks **
509 
510  if ( phi < 0 ) phi+=TMath::TwoPi();
511 
512  rad = Radius(etaC, phiC, eta, phi);
513 
514  // ** Exclude tracks too close to the candidate, inactive by default **
515 
516  if(rad < fDistMinToTrigger) continue ;
517 
518  // ** For the background out of cone **
519 
520  if(rad > fConeSize)
521  {
522  if(eta > (etaC-fConeSize) && eta < (etaC+fConeSize)) phiBandPtSumTrack += pt;
523  if(phi > (phiC-fConeSize) && phi < (phiC+fConeSize)) etaBandPtSumTrack += pt;
524  }
525 
526  // ** For the isolated particle **
527 
528  // Only loop the particle at the same side of candidate
529  if(TMath::Abs(phi-phiC) > TMath::PiOver2()) continue ;
530 
531 // // If at the same side has particle larger than candidate,
532 // // then candidate can not be the leading, skip such events
533 // if(pt > ptC)
534 // {
535 // n = -1;
536 // nfrac = -1;
537 // coneptsumTrack = -1;
538 // isolated = kFALSE;
539 //
540 // pCandidate->SetLeadingParticle(kFALSE);
541 //
542 // if(bFillAOD && reftracks)
543 // {
544 // reftracks->Clear();
545 // delete reftracks;
546 // }
547 //
548 // return ;
549 // }
550 
551  AliDebug(2,Form("\t Track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad));
552 
553  //
554  // Select tracks inside the isolation radius
555  //
556  if(rad < fConeSize)
557  {
558  AliDebug(2,"Inside candidate cone");
559 
560  if(bFillAOD)
561  {
562  ntrackrefs++;
563  if(ntrackrefs == 1)
564  {
565  reftracks = new TObjArray(0);
566  //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
567  TString tempo(aodArrayRefName) ;
568  tempo += "Tracks" ;
569  reftracks->SetName(tempo);
570  reftracks->SetOwner(kFALSE);
571  }
572  reftracks->Add(track);
573  }
574 
575  coneptsumTrack+=pt;
576 
577  if( ptLead < pt ) ptLead = pt;
578 
579 // // *Before*, count particles in cone
580 // if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
581 //
582 // //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
583 // if(fFracIsThresh)
584 // {
585 // if( fPtFraction*ptC < fPtThreshold )
586 // {
587 // if( pt > fPtThreshold ) nfrac++ ;
588 // }
589 // else
590 // {
591 // if( pt > fPtFraction*ptC ) nfrac++;
592 // }
593 // }
594 // else
595 // {
596 // if( pt > fPtFraction*ptC ) nfrac++;
597 // }
598 
599  } // Inside cone
600 
601  }// charged particle loop
602 
603  }//Tracks
604 
605  // --------------------------------
606  // Check calorimeter clusters in cone.
607  // --------------------------------
608 
609  if(plNe &&
611  {
612 
613  for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
614  {
615  AliVCluster * calo = dynamic_cast<AliVCluster *>(plNe->At(ipr)) ;
616 
617  if(calo)
618  {
619  // Get the index where the cluster comes, to retrieve the corresponding vertex
620  Int_t evtIndex = 0 ;
621  if (reader->GetMixedEvent())
622  evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
623 
624 
625  // Do not count the candidate (photon or pi0) or the daughters of the candidate
626  if(calo->GetID() == pCandidate->GetCaloLabel(0) ||
627  calo->GetID() == pCandidate->GetCaloLabel(1) ) continue ;
628 
629  // Skip matched clusters with tracks in case of neutral+charged analysis
631  {
633  pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
634  }
635 
636  // Assume that come from vertex in straight line
637  calo->GetMomentum(fMomentum,reader->GetVertex(evtIndex)) ;
638 
639  pt = fMomentum.Pt() ;
640  eta = fMomentum.Eta() ;
641  phi = fMomentum.Phi() ;
642  }
643  else
644  {// Mixed event stored in AliCaloTrackParticles
645  AliCaloTrackParticle * calomix = dynamic_cast<AliCaloTrackParticle*>(plNe->At(ipr)) ;
646  if(!calomix)
647  {
648  AliWarning("Wrong calo data type, continue");
649  continue;
650  }
651 
652  pt = calomix->Pt();
653  eta = calomix->Eta();
654  phi = calomix->Phi() ;
655  }
656 
657  // ** Calculate distance between candidate and tracks **
658 
659  if( phi < 0 ) phi+=TMath::TwoPi();
660 
661  rad = Radius(etaC, phiC, eta, phi);
662 
663  // ** Exclude clusters too close to the candidate, inactive by default **
664 
665  if(rad < fDistMinToTrigger) continue ;
666 
667  // ** For the background out of cone **
668 
669  if(rad > fConeSize)
670  {
671  if(eta > (etaC-fConeSize) && eta < (etaC+fConeSize)) phiBandPtSumCluster += pt;
672  if(phi > (phiC-fConeSize) && phi < (phiC+fConeSize)) etaBandPtSumCluster += pt;
673  }
674 
675  // ** For the isolated particle **
676 
677  // Only loop the particle at the same side of candidate
678  if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
679 
680 // // If at the same side has particle larger than candidate,
681 // // then candidate can not be the leading, skip such events
682 // if(pt > ptC)
683 // {
684 // n = -1;
685 // nfrac = -1;
686 // coneptsumCluster = -1;
687 // isolated = kFALSE;
688 //
689 // pCandidate->SetLeadingParticle(kFALSE);
690 //
691 // if(bFillAOD)
692 // {
693 // if(reftracks)
694 // {
695 // reftracks ->Clear();
696 // delete reftracks;
697 // }
698 //
699 // if(refclusters)
700 // {
701 // refclusters->Clear();
702 // delete refclusters;
703 // }
704 // }
705 // return ;
706 // }
707 
708  AliDebug(2,Form("\t Cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad));
709 
710  //
711  // Select calorimeter clusters inside the isolation radius
712  //
713  if(rad < fConeSize)
714  {
715  AliDebug(2,"Inside candidate cone");
716 
717  if(bFillAOD)
718  {
719  nclusterrefs++;
720  if(nclusterrefs==1)
721  {
722  refclusters = new TObjArray(0);
723  //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
724  TString tempo(aodArrayRefName) ;
725  tempo += "Clusters" ;
726  refclusters->SetName(tempo);
727  refclusters->SetOwner(kFALSE);
728  }
729  refclusters->Add(calo);
730  }
731 
732  coneptsumCluster+=pt;
733 
734  if( ptLead < pt ) ptLead = pt;
735 
736 // // *Before*, count particles in cone
737 // if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
738 //
739 // //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
740 // if(fFracIsThresh)
741 // {
742 // if( fPtFraction*ptC < fPtThreshold )
743 // {
744 // if( pt > fPtThreshold ) nfrac++ ;
745 // }
746 // else
747 // {
748 // if( pt > fPtFraction*ptC ) nfrac++;
749 // }
750 // }
751 // else
752 // {
753 // if( pt > fPtFraction*ptC ) nfrac++;
754 // }
755 
756  }//in cone
757 
758  }// neutral particle loop
759 
760  }//neutrals
761 
762  //Add reference arrays to AOD when filling AODs only
763  if(bFillAOD)
764  {
765  if(refclusters) pCandidate->AddObjArray(refclusters);
766  if(reftracks) pCandidate->AddObjArray(reftracks);
767  }
768 
769  coneptsum = coneptsumCluster + coneptsumTrack;
770 
771  // *Now*, just check the leading particle in the cone if the threshold is passed
772  if(ptLead > fPtThreshold && ptLead < fPtThresholdMax) n = 1;
773 
774  //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
775  if(fFracIsThresh)
776  {
777  if( fPtFraction*ptC < fPtThreshold )
778  {
779  if( ptLead > fPtThreshold ) nfrac = 1 ;
780  }
781  else
782  {
783  if( ptLead > fPtFraction*ptC ) nfrac = 1;
784  }
785  }
786  else
787  {
788  if( ptLead > fPtFraction*ptC ) nfrac = 1;
789  }
790 
791  //-------------------------------------------------------------------
792  // Check isolation, depending on selected isolation criteria requested
793 
794  if( fICMethod == kPtThresIC)
795  {
796  if( n == 0 ) isolated = kTRUE ;
797 
798  AliDebug(1,Form("pT Cand %2.2f, pT Lead %2.2f, %2.2f<pT Lead< %2.2f, isolated %d",
799  ptC,ptLead,fPtThreshold,fPtThresholdMax,isolated));
800  }
801  else if( fICMethod == kSumPtIC )
802  {
803  if( coneptsum > fSumPtThreshold &&
804  coneptsum < fSumPtThresholdMax )
805  isolated = kFALSE ;
806  else
807  isolated = kTRUE ;
808 
809  AliDebug(1,Form("pT Cand %2.2f, SumPt %2.2f, %2.2f<Sum pT< %2.2f, isolated %d",
810  ptC,ptLead,fSumPtThreshold,fSumPtThresholdMax,isolated));
811  }
812  else if( fICMethod == kPtFracIC )
813  {
814  if(nfrac == 0 ) isolated = kTRUE ;
815  }
816  else if( fICMethod == kSumPtFracIC )
817  {
818  //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
819  // printf("photon analysis IsDataMC() ?%i\n",IsDataMC());
820  if( fFracIsThresh )
821  {
822  if( fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold ) isolated = kTRUE ;
823  if( fPtFraction*ptC > fSumPtThreshold && coneptsum < fPtFraction*ptC ) isolated = kTRUE ;
824  }
825  else
826  {
827  if( coneptsum < fPtFraction*ptC ) isolated = kTRUE ;
828  }
829  }
830  else if( fICMethod == kSumDensityIC )
831  {
832  // Get good cell density (number of active cells over all cells in cone)
833  // and correct energy in cone
834 
835  Float_t cellDensity = GetCellDensity(pCandidate,reader);
836 
837  if( coneptsum < fSumPtThreshold*cellDensity )
838  isolated = kTRUE;
839  }
840  else if( fICMethod == kSumBkgSubIC )
841  {
842  Double_t coneptsumBkg = 0.;
843  Float_t etaBandPtSumTrackNorm = 0;
844  Float_t phiBandPtSumTrackNorm = 0;
845  Float_t etaBandPtSumClusterNorm = 0;
846  Float_t phiBandPtSumClusterNorm = 0;
847 
848  Float_t excessFracEtaTrack = 1;
849  Float_t excessFracPhiTrack = 1;
850  Float_t excessFracEtaCluster = 1;
851  Float_t excessFracPhiCluster = 1;
852 
853  // Normalize background to cone area
854  if (fPartInCone != kOnlyCharged )
855  CalculateUEBandClusterNormalization(reader, etaC, phiC,
856  phiBandPtSumCluster , etaBandPtSumCluster,
857  phiBandPtSumClusterNorm, etaBandPtSumClusterNorm,
858  excessFracEtaCluster , excessFracPhiCluster );
859 
860  if (fPartInCone != kOnlyNeutral )
861  CalculateUEBandTrackNormalization(reader, etaC, phiC,
862  phiBandPtSumTrack , etaBandPtSumTrack ,
863  phiBandPtSumTrackNorm, etaBandPtSumTrackNorm,
864  excessFracEtaTrack , excessFracPhiTrack );
865 
866  if (fPartInCone == kOnlyCharged ) coneptsumBkg = etaBandPtSumTrackNorm;
867  else if(fPartInCone == kOnlyNeutral ) coneptsumBkg = etaBandPtSumClusterNorm;
868  else if(fPartInCone == kNeutralAndCharged ) coneptsumBkg = etaBandPtSumClusterNorm + etaBandPtSumTrackNorm;
869 
870  //coneptsumCluster*=(coneBadCellsCoeff*excessFracEtaCluster*excessFracPhiCluster) ; // apply this correction earlier???
871  // line commented out in last modif!!!
872 
873  coneptsum = coneptsumCluster+coneptsumTrack;
874 
875  coneptsum -= coneptsumBkg;
876 
877  if( coneptsum > fSumPtThreshold && coneptsum < fSumPtThresholdMax )
878  isolated = kFALSE ;
879  else
880  isolated = kTRUE ;
881 
882  }
883 }
884 
885 //_____________________________________________________
887 //_____________________________________________________
889 {
890  if(! opt)
891  return;
892 
893  printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
894 
895  printf("IC method = %d\n", fICMethod ) ;
896  printf("Cone Size = %1.2f\n", fConeSize ) ;
897  printf("pT threshold = >%2.1f;<%2.1f\n", fPtThreshold , fPtThresholdMax) ;
898  printf("Sum pT threshold = >%2.1f;<%2.1f\n", fSumPtThreshold,fSumPtThresholdMax) ;
899  printf("pT fraction = %3.1f\n", fPtFraction ) ;
900  printf("particle type in cone = %d\n", fPartInCone ) ;
901  printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
902  printf("minimum distance to candidate, R>%1.2f\n",fDistMinToTrigger);
903  printf(" \n") ;
904 }
905 
906 //______________________________________________________________
912 //______________________________________________________________
914  Float_t eta , Float_t phi) const
915 {
916  Float_t dEta = etaC-eta;
917  Float_t dPhi = phiC-phi;
918 
919  if(TMath::Abs(dPhi) >= TMath::Pi())
920  dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
921 
922  return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
923 }
924 
925 
926 
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
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
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow, Int_t status) 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...