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