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