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