AliPhysics  a76316e (a76316e)
AliEMCALRecoUtils.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 includes
17 #include <TH2F.h>
18 #include <TArrayI.h>
19 #include <TArrayF.h>
20 #include <TObjArray.h>
21 
22 // STEER includes
23 #include "AliVCluster.h"
24 #include "AliVCaloCells.h"
25 #include "AliLog.h"
26 #include "AliPID.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliESDtrack.h"
30 #include "AliAODTrack.h"
31 #include "AliExternalTrackParam.h"
32 #include "AliESDfriendTrack.h"
33 #include "AliTrackerBase.h"
34 #include "AliMCEvent.h"
35 
36 // EMCAL includes
37 #include "AliEMCALRecoUtils.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliTrackerBase.h"
40 #include "AliEMCALPIDUtils.h"
41 
43 ClassImp(AliEMCALRecoUtils) ;
45 
51 //_____________________________________
52 AliEMCALRecoUtils::AliEMCALRecoUtils():
53  fParticleType(0), fPosAlgo(0),
54  fW0(0), fShowerShapeCellLocationType(0),
55  fNonLinearityFunction(0), fNonLinearThreshold(0),
56  fSmearClusterEnergy(kFALSE), fRandom(),
57  fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
58  fConstantTimeShift(0), fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fLowGain(kFALSE),
59  fUseL1PhaseInTimeRecalibration(kFALSE), fEMCALL1PhaseInTimeRecalibration(),
60  fUseRunCorrectionFactors(kFALSE),
61  fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
62  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
63  fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
64  fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
65  fPIDUtils(), fAODFilterMask(0),
66  fAODHybridTracks(0), fAODTPCOnlyTracks(0),
67  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
68  fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
69  fCutR(0), fCutEta(0), fCutPhi(0),
70  fClusterWindow(0), fMass(0),
71  fStepSurface(0), fStepCluster(0),
72  fITSTrackSA(kFALSE), fEMCalSurfaceDistance(440.),
73  fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
74  fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
75  fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
76  fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
77  fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE),
78  fNMCGenerToAccept(0), fMCGenerToAcceptForTrack(1)
79 {
80  // Init parameters
81  InitParameters();
82 
83  for(Int_t j = 0; j < 5; j++) fMCGenerToAccept[j] = "";
84 
85  // Track matching arrays init
86  fMatchedTrackIndex = new TArrayI();
87  fMatchedClusterIndex = new TArrayI();
88  fResidualPhi = new TArrayF();
89  fResidualEta = new TArrayF();
90  fPIDUtils = new AliEMCALPIDUtils();
91 
92  fBadStatusSelection[0] = kTRUE;
93  fBadStatusSelection[1] = kTRUE;
94  fBadStatusSelection[2] = kTRUE;
95  fBadStatusSelection[3] = kTRUE;
96 }
97 
98 //
99 // Copy constructor.
100 //
101 //______________________________________________________________________
102 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
103 : AliEMCALRecoUtilsBase(reco),
104  fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo),
105  fW0(reco.fW0), fShowerShapeCellLocationType(reco.fShowerShapeCellLocationType),
106  fNonLinearityFunction(reco.fNonLinearityFunction), fNonLinearThreshold(reco.fNonLinearThreshold),
107  fSmearClusterEnergy(reco.fSmearClusterEnergy), fRandom(),
108  fCellsRecalibrated(reco.fCellsRecalibrated),
109  fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(NULL),
110  fConstantTimeShift(reco.fConstantTimeShift),
111  fTimeRecalibration(reco.fTimeRecalibration), fEMCALTimeRecalibrationFactors(NULL),
112  fLowGain(reco.fLowGain),
113  fUseL1PhaseInTimeRecalibration(reco.fUseL1PhaseInTimeRecalibration),
114  fEMCALL1PhaseInTimeRecalibration(reco.fEMCALL1PhaseInTimeRecalibration),
115  fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),
116  fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
117  fEMCALBadChannelMap(NULL),
118  fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
119  fRejectExoticCluster(reco.fRejectExoticCluster), fRejectExoticCells(reco.fRejectExoticCells),
120  fExoticCellFraction(reco.fExoticCellFraction), fExoticCellDiffTime(reco.fExoticCellDiffTime),
121  fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
122  fPIDUtils(reco.fPIDUtils), fAODFilterMask(reco.fAODFilterMask),
123  fAODHybridTracks(reco.fAODHybridTracks), fAODTPCOnlyTracks(reco.fAODTPCOnlyTracks),
124  fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
125  fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
126  fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
127  fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
128  fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
129  fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
130  fClusterWindow(reco.fClusterWindow),
131  fMass(reco.fMass), fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
132  fITSTrackSA(reco.fITSTrackSA), fEMCalSurfaceDistance(440.),
133  fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
134  fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
135  fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
136  fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
137  fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
138  fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
139  fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone), fCutRequireITSpureSA(reco.fCutRequireITSpureSA),
140  fNMCGenerToAccept(reco.fNMCGenerToAccept), fMCGenerToAcceptForTrack(reco.fMCGenerToAcceptForTrack)
141 {
142  for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
143  fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
144  for (Int_t i = 0; i < 10 ; i++){ fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
145  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
146  for (Int_t j = 0; j < 5 ; j++) { fMCGenerToAccept[j] = reco.fMCGenerToAccept[j] ; }
147  for (Int_t j = 0; j < 4 ; j++) { fBadStatusSelection[j] = reco.fBadStatusSelection[j] ; }
148 
149  if(reco.fEMCALBadChannelMap) {
150  // Copy constructor - not taking ownership over calibration histograms
151  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
152  fEMCALBadChannelMap->SetOwner(false);
153  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
154  }
155 
156  if(reco.fEMCALRecalibrationFactors) {
157  // Copy constructor - not taking ownership over calibration histograms
158  fEMCALRecalibrationFactors = new TObjArray(reco.fEMCALRecalibrationFactors->GetEntries());
159  fEMCALRecalibrationFactors->SetOwner(false);
160  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
161  }
162 
163  if(reco.fEMCALTimeRecalibrationFactors) {
164  // Copy constructor - not taking ownership over calibration histograms
165  fEMCALTimeRecalibrationFactors = new TObjArray(reco.fEMCALTimeRecalibrationFactors->GetEntries());
166  fEMCALTimeRecalibrationFactors->SetOwner(false);
167  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
168  }
169 }
170 
174 //______________________________________________________________________
175 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
176 {
177  if (this == &reco)return *this;
178  ((TNamed *)this)->operator=(reco);
179 
180  for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
181  fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
182  for (Int_t i = 0; i < 10 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
183  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
184 
185  fParticleType = reco.fParticleType;
186  fPosAlgo = reco.fPosAlgo;
187  fW0 = reco.fW0;
188  fShowerShapeCellLocationType = reco.fShowerShapeCellLocationType;
189 
190  fNonLinearityFunction = reco.fNonLinearityFunction;
191  fNonLinearThreshold = reco.fNonLinearThreshold;
192  fSmearClusterEnergy = reco.fSmearClusterEnergy;
193 
194  fCellsRecalibrated = reco.fCellsRecalibrated;
195  fRecalibration = reco.fRecalibration;
196 
197  fConstantTimeShift = reco.fConstantTimeShift;
198  fTimeRecalibration = reco.fTimeRecalibration;
199  fLowGain = reco.fLowGain;
200 
201  fUseL1PhaseInTimeRecalibration = reco.fUseL1PhaseInTimeRecalibration;
202  fEMCALL1PhaseInTimeRecalibration = reco.fEMCALL1PhaseInTimeRecalibration;
203 
204  fUseRunCorrectionFactors = reco.fUseRunCorrectionFactors;
205 
206  fRemoveBadChannels = reco.fRemoveBadChannels;
207  fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
208 
209  fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
210  fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
211 
212  fRejectExoticCluster = reco.fRejectExoticCluster;
213  fRejectExoticCells = reco.fRejectExoticCells;
214  fExoticCellFraction = reco.fExoticCellFraction;
215  fExoticCellDiffTime = reco.fExoticCellDiffTime;
216  fExoticCellMinAmplitude = reco.fExoticCellMinAmplitude;
217 
218  fPIDUtils = reco.fPIDUtils;
219 
220  fAODFilterMask = reco.fAODFilterMask;
221  fAODHybridTracks = reco.fAODHybridTracks;
222  fAODTPCOnlyTracks = reco.fAODTPCOnlyTracks;
223 
224  fCutEtaPhiSum = reco.fCutEtaPhiSum;
225  fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
226  fCutR = reco.fCutR;
227  fCutEta = reco.fCutEta;
228  fCutPhi = reco.fCutPhi;
229  fClusterWindow = reco.fClusterWindow;
230  fMass = reco.fMass;
231  fStepSurface = reco.fStepSurface;
232  fStepCluster = reco.fStepCluster;
233  fITSTrackSA = reco.fITSTrackSA;
234  fEMCalSurfaceDistance = reco.fEMCalSurfaceDistance;
235 
236  fTrackCutsType = reco.fTrackCutsType;
237  fCutMinTrackPt = reco.fCutMinTrackPt;
238  fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
239  fCutMinNClusterITS = reco.fCutMinNClusterITS;
240  fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
241  fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
242  fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
243  fCutRequireITSRefit = reco.fCutRequireITSRefit;
244  fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
245  fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
246  fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
247  fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
248  fCutRequireITSStandAlone = reco.fCutRequireITSStandAlone;
249  fCutRequireITSpureSA = reco.fCutRequireITSpureSA;
250 
251  fNMCGenerToAccept = reco.fNMCGenerToAccept;
252  fMCGenerToAcceptForTrack = reco.fMCGenerToAcceptForTrack;
253  for (Int_t j = 0; j < 5 ; j++)
254  fMCGenerToAccept[j] = reco.fMCGenerToAccept[j];
255 
256  //
257  // Assign or copy construct the different TArrays
258  //
259  if (reco.fResidualEta)
260  {
261  if (fResidualEta)
262  *fResidualEta = *reco.fResidualEta;
263  else
264  fResidualEta = new TArrayF(*reco.fResidualEta);
265  }
266  else
267  {
268  if (fResidualEta) delete fResidualEta;
269  fResidualEta = 0;
270  }
271 
272  if (reco.fResidualPhi)
273  {
274  if (fResidualPhi)
275  *fResidualPhi = *reco.fResidualPhi;
276  else
277  fResidualPhi = new TArrayF(*reco.fResidualPhi);
278  }
279  else
280  {
281  if (fResidualPhi) delete fResidualPhi;
282  fResidualPhi = 0;
283  }
284 
285  if (reco.fMatchedTrackIndex)
286  {
287  if (fMatchedTrackIndex)
288  *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
289  else
290  fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
291  }
292  else
293  {
294  if (fMatchedTrackIndex) delete fMatchedTrackIndex;
295  fMatchedTrackIndex = 0;
296  }
297 
298  if (reco.fMatchedClusterIndex)
299  {
300  if (fMatchedClusterIndex)
301  *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
302  else
303  fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
304  }
305  else
306  {
307  if (fMatchedClusterIndex) delete fMatchedClusterIndex;
308  fMatchedClusterIndex = 0;
309  }
310 
311  for (Int_t j = 0; j < 4 ; j++)
312  fBadStatusSelection[j] = reco.fBadStatusSelection[j] ;
313 
314  if(fEMCALBadChannelMap) delete fEMCALBadChannelMap;
315  if(reco.fEMCALBadChannelMap) {
316  // Copy constructor - not taking ownership over calibration histograms
317  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
318  fEMCALBadChannelMap->SetOwner(false);
319  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
320  }
321 
322  if(fEMCALRecalibrationFactors) delete fEMCALRecalibrationFactors;
323  if(reco.fEMCALRecalibrationFactors) {
324  // Copy constructor - not taking ownership over calibration histograms
325  fEMCALRecalibrationFactors = new TObjArray(reco.fEMCALRecalibrationFactors->GetEntries());
326  fEMCALRecalibrationFactors->SetOwner(false);
327  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
328  }
329 
330  if(fEMCALTimeRecalibrationFactors) delete fEMCALTimeRecalibrationFactors;
331  if(reco.fEMCALTimeRecalibrationFactors) {
332  // Copy constructor - not taking ownership over calibration histograms
333  fEMCALTimeRecalibrationFactors = new TObjArray(reco.fEMCALTimeRecalibrationFactors->GetEntries());
334  fEMCALTimeRecalibrationFactors->SetOwner(false);
335  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
336  }
337 
338  if(this->fEMCALL1PhaseInTimeRecalibration) delete fEMCALL1PhaseInTimeRecalibration;
339  if(reco.fEMCALL1PhaseInTimeRecalibration) {
340  // Copy constructor - not taking ownership over calibration histograms
341  fEMCALL1PhaseInTimeRecalibration = new TObjArray(reco.fEMCALL1PhaseInTimeRecalibration->GetEntries());
342  fEMCALL1PhaseInTimeRecalibration->SetOwner(false);
343  for(int ism = 0; ism < reco.fEMCALL1PhaseInTimeRecalibration->GetEntries(); ism++) fEMCALL1PhaseInTimeRecalibration->AddAt(reco.fEMCALL1PhaseInTimeRecalibration->At(ism), ism);
344  }
345 
346  return *this;
347 }
348 
352 //_____________________________________
353 AliEMCALRecoUtils::~AliEMCALRecoUtils()
354 {
355  if (fEMCALRecalibrationFactors)
356  {
357  delete fEMCALRecalibrationFactors;
358  }
359 
360  if (fEMCALTimeRecalibrationFactors)
361  {
362  delete fEMCALTimeRecalibrationFactors;
363  }
364 
365  if(fEMCALL1PhaseInTimeRecalibration)
366  {
367  delete fEMCALL1PhaseInTimeRecalibration;
368  }
369 
370  if (fEMCALBadChannelMap)
371  {
372  delete fEMCALBadChannelMap;
373  }
374 
375  delete fMatchedTrackIndex ;
376  delete fMatchedClusterIndex ;
377  delete fResidualEta ;
378  delete fResidualPhi ;
379  delete fPIDUtils ;
380 
381  InitTrackCuts();
382 }
383 
396 //_______________________________________________________________________________
397 Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(Int_t absID, Int_t bc,
398  Float_t & amp, Double_t & time,
399  AliVCaloCells* cells)
400 {
401  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
402 
403  if(!geom)
404  {
405  AliError("No instance of the geometry is available");
406  return kFALSE;
407  }
408 
409  if ( absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules() )
410  return kFALSE;
411 
412  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1, status=0;
413 
414  if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
415  {
416  // cell absID does not exist
417  amp=0; time = 1.e9;
418  return kFALSE;
419  }
420 
421  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
422 
423  // Do not include bad channels found in analysis,
424  if ( IsBadChannelsRemovalSwitchedOn() )
425  {
426  Bool_t bad = GetEMCALChannelStatus(imod, ieta, iphi,status);
427 
428  if ( status > 0 )
429  AliDebug(1,Form("Channel absId %d, status %d, set as bad %d",absID, status, bad));
430 
431  if ( bad ) return kFALSE;
432  }
433 
434  //Recalibrate energy
435  amp = cells->GetCellAmplitude(absID);
436  if (!fCellsRecalibrated && IsRecalibrationOn())
437  amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
438 
439  // Recalibrate time
440  time = cells->GetCellTime(absID);
441  time-=fConstantTimeShift*1e-9; // only in case of old Run1 simulation
442  Bool_t isLowGain = !(cells->GetCellHighGain(absID));//HG = false -> LG = true
443 
444  RecalibrateCellTime(absID,bc,time,isLowGain);
445 
446  //Recalibrate time with L1 phase
447  RecalibrateCellTimeL1Phase(imod, bc, time);
448 
449  return kTRUE;
450 }
451 
461 //_____________________________________________________________________________
462 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom,
463  const AliVCluster* cluster,
464  AliVCaloCells* cells)
465 {
466  if (!cluster)
467  {
468  AliInfo("Cluster pointer null!");
469  return kFALSE;
470  }
471 
472  // If the distance to the border is 0 or negative just exit accept all clusters
473  if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
474  return kTRUE;
475 
476  Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
477  Bool_t shared = kFALSE;
478  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
479 
480  AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
481  absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
482 
483  if (absIdMax==-1) return kFALSE;
484 
485  // Check if the cell is close to the borders:
486  Bool_t okrow = kFALSE;
487  Bool_t okcol = kFALSE;
488 
489  if (iSM < 0 || iphi < 0 || ieta < 0 )
490  {
491  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
492  iSM,ieta,iphi));
493  return kFALSE; // trick coverity
494  }
495 
496  // Check rows/phi
497  Int_t iPhiLast = 24;
498  if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
499  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
500  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Ext ) iPhiLast /= 3;// 1/3 sm case
501 
502  if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
503 
504  // Check columns/eta
505  Int_t iEtaLast = 48;
506  if ( !fNoEMCALBorderAtEta0 || geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard )
507  {
508  // consider inner border
509  if ( geom->IsDCALSM(iSM) ) iEtaLast = iEtaLast*2/3;
510 
511  if ( ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder ) okcol = kTRUE;
512  }
513  else
514  {
515  if (iSM%2==0)
516  {
517  if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
518  }
519  else
520  {
521  if (ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
522  }
523  }//eta 0 not checked
524 
525  AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
526  fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
527 
528  if (okcol && okrow)
529  {
530  return kTRUE;
531  }
532  else
533  {
534  AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
535 
536  return kFALSE;
537  }
538 }
539 
550 //_______________________________________________________________________________
551 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom,
552  const UShort_t* cellList,
553  Int_t nCells)
554 {
555  if (!fRemoveBadChannels) return kFALSE;
556  if (!fEMCALBadChannelMap) return kFALSE;
557 
558  Int_t icol = -1;
559  Int_t irow = -1;
560  Int_t imod = -1;
561  for (Int_t iCell = 0; iCell<nCells; iCell++)
562  {
563  //Get the column and row
564  Int_t iTower = -1, iIphi = -1, iIeta = -1;
565  geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
566 
567  if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
568 
569  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
570 
571  Int_t status = 0;
572  if (GetEMCALChannelStatus(imod, icol, irow, status))
573  {
574  AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d, status %d\n",imod, icol, irow, status));
575  return kTRUE;
576  }
577  }// cell cluster loop
578 
579  return kFALSE;
580 }
581 
593 //___________________________________________________________________________
594 Float_t AliEMCALRecoUtils::GetECross(Int_t absID, Double_t tcell,
595  AliVCaloCells* cells, Int_t bc)
596 {
597  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
598 
599  if(!geom)
600  {
601  AliError("No instance of the geometry is available");
602  return -1;
603  }
604 
605  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
606  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
607  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
608 
609  // Get close cells index, energy and time, not in corners
610 
611  Int_t absID1 = -1;
612  Int_t absID2 = -1;
613 
614  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
615  if ( iphi > 0 ) absID2 = geom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
616 
617  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
618 
619  Int_t absID3 = -1;
620  Int_t absID4 = -1;
621 
622  if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
623  {
624  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
625  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
626  }
627  else if ( ieta == 0 && imod%2 )
628  {
629  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
630  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
631  }
632  else
633  {
634  if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
635  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
636  if ( ieta > 0 )
637  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
638  }
639 
640  //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
641 
642  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
643  Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
644 
645  AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
646  AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
647  AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
648  AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
649 
650  if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
651  if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
652  if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
653  if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
654 
655  return ecell1+ecell2+ecell3+ecell4;
656 }
657 
668 //_____________________________________________________________________________________________
669 Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
670 {
671  if (!fRejectExoticCells) return kFALSE;
672 
673  Float_t ecell = 0;
674  Double_t tcell = 0;
675  Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
676 
677  if (!accept) return kTRUE; // reject this cell
678 
679  if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
680 
681  Float_t eCross = GetECross(absID,tcell,cells,bc);
682 
683  if (1-eCross/ecell > fExoticCellFraction)
684  {
685  AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
686  absID,ecell,eCross,1-eCross/ecell));
687  return kTRUE;
688  }
689 
690  return kFALSE;
691 }
692 
702 //___________________________________________________________________
703 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
704  AliVCaloCells *cells,
705  Int_t bc)
706 {
707  if (!cluster)
708  {
709  AliInfo("Cluster pointer null!");
710  return kFALSE;
711  }
712 
713  if (!fRejectExoticCluster) return kFALSE;
714 
715  // Get highest energy tower
716  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
717 
718  if(!geom)
719  {
720  AliError("No instance of the geometry is available");
721  return kFALSE;
722  }
723 
724  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
725  Bool_t shared = kFALSE;
726  GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
727 
728  return IsExoticCell(absId,cells,bc);
729 }
730 
739 //_______________________________________________________________________
740 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
741 {
742  if (!cluster)
743  {
744  AliInfo("Cluster pointer null!");
745  return 0;
746  }
747 
748  Float_t energy = cluster->E() ;
749  Float_t rdmEnergy = energy ;
750 
751  if (fSmearClusterEnergy)
752  {
753  rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
754  fSmearClusterParam[1] * energy +
755  fSmearClusterParam[2] );
756  AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
757  }
758 
759  return rdmEnergy;
760 }
761 
770 //____________________________________________________________________________
771 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
772 {
773  if (!cluster)
774  {
775  AliInfo("Cluster pointer null!");
776  return 0;
777  }
778 
779  Float_t energy = cluster->E();
780  if (energy < 0.100)
781  {
782  // Clusters with less than 100 MeV or negative are not possible
783  AliInfo(Form("Too low cluster energy!, E = %f < 0.100 GeV",energy));
784  return 0;
785  }
786  else if(energy > 300.)
787  {
788  AliInfo(Form("Too high cluster energy!, E = %f GeV, do not apply non linearity",energy));
789  return energy;
790  }
791 
792  switch (fNonLinearityFunction)
793  {
794  case kPi0MC:
795  {
796  //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
797  //fNonLinearityParams[0] = 1.014;
798  //fNonLinearityParams[1] =-0.03329;
799  //fNonLinearityParams[2] =-0.3853;
800  //fNonLinearityParams[3] = 0.5423;
801  //fNonLinearityParams[4] =-0.4335;
802  energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
803  ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
804  exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
805  break;
806  }
807 
808  case kPi0MCv2:
809  {
810  //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
811  //fNonLinearityParams[0] = 3.11111e-02;
812  //fNonLinearityParams[1] =-5.71666e-02;
813  //fNonLinearityParams[2] = 5.67995e-01;
814 
815  energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
816  break;
817  }
818 
819  case kPi0MCv3:
820  {
821  //Same as beam test corrected, change parameters
822  //fNonLinearityParams[0] = 9.81039e-01
823  //fNonLinearityParams[1] = 1.13508e-01;
824  //fNonLinearityParams[2] = 1.00173e+00;
825  //fNonLinearityParams[3] = 9.67998e-02;
826  //fNonLinearityParams[4] = 2.19381e+02;
827  //fNonLinearityParams[5] = 6.31604e+01;
828  //fNonLinearityParams[6] = 1;
829  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
830 
831  break;
832  }
833 
834 
835  case kPi0GammaGamma:
836  {
837  //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
838  //fNonLinearityParams[0] = 1.04;
839  //fNonLinearityParams[1] = -0.1445;
840  //fNonLinearityParams[2] = 1.046;
841  energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
842  break;
843  }
844 
845  case kPi0GammaConversion:
846  {
847  //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
848  //fNonLinearityParams[0] = 0.139393/0.1349766;
849  //fNonLinearityParams[1] = 0.0566186;
850  //fNonLinearityParams[2] = 0.982133;
851  energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
852 
853  break;
854  }
855 
856  case kBeamTest:
857  {
858  //From beam test, Alexei's results, for different ZS thresholds
859  // th=30 MeV; th = 45 MeV; th = 75 MeV
860  //fNonLinearityParams[0] = 1.007; 1.003; 1.002
861  //fNonLinearityParams[1] = 0.894; 0.719; 0.797
862  //fNonLinearityParams[2] = 0.246; 0.334; 0.358
863  //Rescale the param[0] with 1.03
864  energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
865 
866  break;
867  }
868 
869  case kBeamTestCorrected:
870  {
871  // From beam test, corrected for material between beam and EMCAL
872  //fNonLinearityParams[0] = 0.99078
873  //fNonLinearityParams[1] = 0.161499;
874  //fNonLinearityParams[2] = 0.655166;
875  //fNonLinearityParams[3] = 0.134101;
876  //fNonLinearityParams[4] = 163.282;
877  //fNonLinearityParams[5] = 23.6904;
878  //fNonLinearityParams[6] = 0.978;
879  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
880 
881  break;
882  }
883 
884  case kBeamTestCorrectedv2:
885  {
886  // From beam test, corrected for material between beam and EMCAL
887  // Different function to kBeamTestCorrected
888  //fNonLinearityParams[0] = 0.983504;
889  //fNonLinearityParams[1] = 0.210106;
890  //fNonLinearityParams[2] = 0.897274;
891  //fNonLinearityParams[3] = 0.0829064;
892  //fNonLinearityParams[4] = 152.299;
893  //fNonLinearityParams[5] = 31.5028;
894  //fNonLinearityParams[6] = 0.968;
895  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
896 
897  break;
898  }
899 
900  case kBeamTestCorrectedv3:
901  {
902  // Same function as kBeamTestCorrectedv2, different default parametrization.
903  //fNonLinearityParams[0] = 0.976941;
904  //fNonLinearityParams[1] = 0.162310;
905  //fNonLinearityParams[2] = 1.08689;
906  //fNonLinearityParams[3] = 0.0819592;
907  //fNonLinearityParams[4] = 152.338;
908  //fNonLinearityParams[5] = 30.9594;
909  //fNonLinearityParams[6] = 0.9615;
910  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
911 
912  break;
913  }
914 
915  case kSDMv5:
916  {
917  //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
918  //fNonLinearityParams[0] = 1.0;
919  //fNonLinearityParams[1] = 6.64778e-02;
920  //fNonLinearityParams[2] = 1.570;
921  //fNonLinearityParams[3] = 9.67998e-02;
922  //fNonLinearityParams[4] = 2.19381e+02;
923  //fNonLinearityParams[5] = 6.31604e+01;
924  //fNonLinearityParams[6] = 1.01286;
925  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5])))) * (0.964 + exp(-3.132-0.435*energy*2.0));
926 
927  break;
928  }
929 
930  case kPi0MCv5:
931  {
932  //Based on comparing MC truth information to the reconstructed energy of clusters.
933  //fNonLinearityParams[0] = 1.0;
934  //fNonLinearityParams[1] = 6.64778e-02;
935  //fNonLinearityParams[2] = 1.570;
936  //fNonLinearityParams[3] = 9.67998e-02;
937  //fNonLinearityParams[4] = 2.19381e+02;
938  //fNonLinearityParams[5] = 6.31604e+01;
939  //fNonLinearityParams[6] = 1.01286;
940  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
941 
942  break;
943  }
944 
945  case kSDMv6:
946  {
947  //Based on fit to the MC/data using kNoCorrection on the data
948  // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
949  // - parameters constrained by the test beam data as well
950  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
951  //fNonLinearityParams[0] = 1.0;
952  //fNonLinearityParams[1] = 0.237767;
953  //fNonLinearityParams[2] = 0.651203;
954  //fNonLinearityParams[3] = 0.183741;
955  //fNonLinearityParams[4] = 155.427;
956  //fNonLinearityParams[5] = 17.0335;
957  //fNonLinearityParams[6] = 0.987054;
958  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
959 
960  break;
961  }
962 
963  case kPi0MCv6:
964  {
965  //Based on comparing MC truth information to the reconstructed energy of clusters.
966  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
967  //fNonLinearityParams[0] = 1.0;
968  //fNonLinearityParams[1] = 0.0797873;
969  //fNonLinearityParams[2] = 1.68322;
970  //fNonLinearityParams[3] = 0.0806098;
971  //fNonLinearityParams[4] = 244.586;
972  //fNonLinearityParams[5] = 116.938;
973  //fNonLinearityParams[6] = 1.00437;
974  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
975 
976  break;
977  }
978 
979  case kPCMv1:
980  {
981  //based on symmetric decays of pi0 meson
982  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
983  // parameters vary from MC to MC
984  //fNonLinearityParams[0] = 0.984876;
985  //fNonLinearityParams[1] = -9.999609;
986  //fNonLinearityParams[2] = -4.999891;
987  //fNonLinearityParams[3] = 0.;
988  //fNonLinearityParams[4] = 0.;
989  //fNonLinearityParams[5] = 0.;
990  //fNonLinearityParams[6] = 0.;
991  energy /= TMath::Power(fNonLinearityParams[0] + exp(fNonLinearityParams[1] + fNonLinearityParams[2]*energy),2);//result coming from calo-conv method
992 
993  break;
994  }
995 
996  case kPCMplusBTCv1:
997  {
998  //convolution of TestBeamCorrectedv3 with PCM method
999  //Based on comparing MC truth information to the reconstructed energy of clusters.
1000  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1001  // parameters vary from MC to MC
1002  //fNonLinearityParams[0] = 0.976941;
1003  //fNonLinearityParams[1] = 0.162310;
1004  //fNonLinearityParams[2] = 1.08689;
1005  //fNonLinearityParams[3] = 0.0819592;
1006  //fNonLinearityParams[4] = 152.338;
1007  //fNonLinearityParams[5] = 30.9594;
1008  //fNonLinearityParams[6] = 0.9615;
1009  //fNonLinearityParams[7] = 0.984876;
1010  //fNonLinearityParams[8] = -9.999609;
1011  //fNonLinearityParams[9] = -4.999891;
1012  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1013  energy /= TMath::Power(fNonLinearityParams[7] + exp(fNonLinearityParams[8] + fNonLinearityParams[9]*energy),2);//result coming from calo-conv method
1014 
1015  break;
1016  }
1017 
1018  case kPCMsysv1:
1019  {
1020  // Systematic variation of kPCMv1
1021  //Based on comparing MC truth information to the reconstructed energy of clusters.
1022  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1023  // parameters vary from MC to MC
1024  //fNonLinearityParams[0] = 0.0;
1025  //fNonLinearityParams[1] = 1.0;
1026  //fNonLinearityParams[2] = 1.0;
1027  //fNonLinearityParams[3] = 0.0;
1028  //fNonLinearityParams[4] = 1.0;
1029  //fNonLinearityParams[5] = 0.0;
1030  //fNonLinearityParams[6] = 0.0;
1031  energy /= TMath::Power( (fNonLinearityParams[0] + fNonLinearityParams[1] * TMath::Power(energy,fNonLinearityParams[2]) ) /
1032  (fNonLinearityParams[3] + fNonLinearityParams[4] * TMath::Power(energy,fNonLinearityParams[5]) ) + fNonLinearityParams[6], 2);//result coming from calo-conv method
1033 
1034  break;
1035  }
1036 
1037 
1038  case kNoCorrection:
1039  AliDebug(2,"No correction on the energy\n");
1040  break;
1041 
1042  }
1043 
1044  return energy;
1045 }
1046 
1051 //__________________________________________________
1052 void AliEMCALRecoUtils::InitNonLinearityParam()
1053 {
1054  if (fNonLinearityFunction == kPi0MC) {
1055  fNonLinearityParams[0] = 1.014;
1056  fNonLinearityParams[1] = -0.03329;
1057  fNonLinearityParams[2] = -0.3853;
1058  fNonLinearityParams[3] = 0.5423;
1059  fNonLinearityParams[4] = -0.4335;
1060  }
1061 
1062  if (fNonLinearityFunction == kPi0MCv2) {
1063  fNonLinearityParams[0] = 3.11111e-02;
1064  fNonLinearityParams[1] =-5.71666e-02;
1065  fNonLinearityParams[2] = 5.67995e-01;
1066  }
1067 
1068  if (fNonLinearityFunction == kPi0MCv3) {
1069  fNonLinearityParams[0] = 9.81039e-01;
1070  fNonLinearityParams[1] = 1.13508e-01;
1071  fNonLinearityParams[2] = 1.00173e+00;
1072  fNonLinearityParams[3] = 9.67998e-02;
1073  fNonLinearityParams[4] = 2.19381e+02;
1074  fNonLinearityParams[5] = 6.31604e+01;
1075  fNonLinearityParams[6] = 1;
1076  }
1077 
1078  if (fNonLinearityFunction == kPi0GammaGamma) {
1079  fNonLinearityParams[0] = 1.04;
1080  fNonLinearityParams[1] = -0.1445;
1081  fNonLinearityParams[2] = 1.046;
1082  }
1083 
1084  if (fNonLinearityFunction == kPi0GammaConversion) {
1085  fNonLinearityParams[0] = 0.139393;
1086  fNonLinearityParams[1] = 0.0566186;
1087  fNonLinearityParams[2] = 0.982133;
1088  }
1089 
1090  if (fNonLinearityFunction == kBeamTest) {
1091  if (fNonLinearThreshold == 30) {
1092  fNonLinearityParams[0] = 1.007;
1093  fNonLinearityParams[1] = 0.894;
1094  fNonLinearityParams[2] = 0.246;
1095  }
1096  if (fNonLinearThreshold == 45) {
1097  fNonLinearityParams[0] = 1.003;
1098  fNonLinearityParams[1] = 0.719;
1099  fNonLinearityParams[2] = 0.334;
1100  }
1101  if (fNonLinearThreshold == 75) {
1102  fNonLinearityParams[0] = 1.002;
1103  fNonLinearityParams[1] = 0.797;
1104  fNonLinearityParams[2] = 0.358;
1105  }
1106  }
1107 
1108  if (fNonLinearityFunction == kBeamTestCorrected) {
1109  fNonLinearityParams[0] = 0.99078;
1110  fNonLinearityParams[1] = 0.161499;
1111  fNonLinearityParams[2] = 0.655166;
1112  fNonLinearityParams[3] = 0.134101;
1113  fNonLinearityParams[4] = 163.282;
1114  fNonLinearityParams[5] = 23.6904;
1115  fNonLinearityParams[6] = 0.978;
1116  }
1117 
1118  if (fNonLinearityFunction == kBeamTestCorrectedv2) {
1119  // Parameters until November 2015, use now kBeamTestCorrectedv3
1120  fNonLinearityParams[0] = 0.983504;
1121  fNonLinearityParams[1] = 0.210106;
1122  fNonLinearityParams[2] = 0.897274;
1123  fNonLinearityParams[3] = 0.0829064;
1124  fNonLinearityParams[4] = 152.299;
1125  fNonLinearityParams[5] = 31.5028;
1126  fNonLinearityParams[6] = 0.968;
1127  }
1128 
1129  if (fNonLinearityFunction == kBeamTestCorrectedv3) {
1130 
1131  // New parametrization of kBeamTestCorrectedv2
1132  // excluding point at 0.5 GeV from Beam Test Data
1133  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1134 
1135  fNonLinearityParams[0] = 0.976941;
1136  fNonLinearityParams[1] = 0.162310;
1137  fNonLinearityParams[2] = 1.08689;
1138  fNonLinearityParams[3] = 0.0819592;
1139  fNonLinearityParams[4] = 152.338;
1140  fNonLinearityParams[5] = 30.9594;
1141  fNonLinearityParams[6] = 0.9615;
1142  }
1143 
1144  if (fNonLinearityFunction == kSDMv5) {
1145  fNonLinearityParams[0] = 1.0;
1146  fNonLinearityParams[1] = 6.64778e-02;
1147  fNonLinearityParams[2] = 1.570;
1148  fNonLinearityParams[3] = 9.67998e-02;
1149  fNonLinearityParams[4] = 2.19381e+02;
1150  fNonLinearityParams[5] = 6.31604e+01;
1151  fNonLinearityParams[6] = 1.01286;
1152  }
1153 
1154  if (fNonLinearityFunction == kPi0MCv5) {
1155  fNonLinearityParams[0] = 1.0;
1156  fNonLinearityParams[1] = 6.64778e-02;
1157  fNonLinearityParams[2] = 1.570;
1158  fNonLinearityParams[3] = 9.67998e-02;
1159  fNonLinearityParams[4] = 2.19381e+02;
1160  fNonLinearityParams[5] = 6.31604e+01;
1161  fNonLinearityParams[6] = 1.01286;
1162  }
1163 
1164  if (fNonLinearityFunction == kSDMv6) {
1165  fNonLinearityParams[0] = 1.0;
1166  fNonLinearityParams[1] = 0.237767;
1167  fNonLinearityParams[2] = 0.651203;
1168  fNonLinearityParams[3] = 0.183741;
1169  fNonLinearityParams[4] = 155.427;
1170  fNonLinearityParams[5] = 17.0335;
1171  fNonLinearityParams[6] = 0.987054;
1172  }
1173 
1174  if (fNonLinearityFunction == kPi0MCv6) {
1175  fNonLinearityParams[0] = 1.0;
1176  fNonLinearityParams[1] = 0.0797873;
1177  fNonLinearityParams[2] = 1.68322;
1178  fNonLinearityParams[3] = 0.0806098;
1179  fNonLinearityParams[4] = 244.586;
1180  fNonLinearityParams[5] = 116.938;
1181  fNonLinearityParams[6] = 1.00437;
1182  }
1183 
1184 if (fNonLinearityFunction == kPCMv1) {
1185  //parameters change from MC production to MC production, they need to set for each period
1186  fNonLinearityParams[0] = 0.984876;
1187  fNonLinearityParams[1] = -9.999609;
1188  fNonLinearityParams[2] = -4.999891;
1189  fNonLinearityParams[3] = 0.;
1190  fNonLinearityParams[4] = 0.;
1191  fNonLinearityParams[5] = 0.;
1192  fNonLinearityParams[6] = 0.;
1193  }
1194 
1195  if (fNonLinearityFunction == kPCMplusBTCv1) {
1196  // test beam corrected values convoluted with symmetric meson decays values
1197  // for test beam:
1198  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1199  // for PCM method:
1200  // https://aliceinfo.cern.ch/Notes/node/211
1201  fNonLinearityParams[0] = 0.976941;
1202  fNonLinearityParams[1] = 0.162310;
1203  fNonLinearityParams[2] = 1.08689;
1204  fNonLinearityParams[3] = 0.0819592;
1205  fNonLinearityParams[4] = 152.338;
1206  fNonLinearityParams[5] = 30.9594;
1207  fNonLinearityParams[6] = 0.9615;
1208  fNonLinearityParams[7] = 0.984876;
1209  fNonLinearityParams[8] = -9.999609;
1210  fNonLinearityParams[9] = -4.999891;
1211  }
1212 
1213  if (fNonLinearityFunction == kPCMsysv1) {
1214  //systematics for kPCMv1
1215  // for PCM method:
1216  // https://aliceinfo.cern.ch/Notes/node/211
1217  fNonLinearityParams[0] = 0.0;
1218  fNonLinearityParams[1] = 1.0;
1219  fNonLinearityParams[2] = 1.0;
1220  fNonLinearityParams[3] = 0.0;
1221  fNonLinearityParams[4] = 1.0;
1222  fNonLinearityParams[5] = 0.0;
1223  fNonLinearityParams[6] = 0.0;
1224  }
1225 }
1226 
1234 //____________________________________________________________________
1235 void AliEMCALRecoUtils::SetEMCALBadChannelStatusSelection(Bool_t all, Bool_t dead, Bool_t hot, Bool_t warm)
1236 {
1237  fBadStatusSelection[0] = all; // declare all as bad if true, never mind the other settings
1238  fBadStatusSelection[1] = dead;
1239  fBadStatusSelection[2] = hot;
1240  fBadStatusSelection[3] = warm;
1241 }
1242 
1253 //____________________________________________________________________
1254 Bool_t AliEMCALRecoUtils::GetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow, Int_t & status) const
1255 {
1256  if(fEMCALBadChannelMap)
1257  status = (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
1258  else
1259  status = 0; // Channel is ok by default
1260 
1261  if ( status == AliCaloCalibPedestal::kAlive )
1262  {
1263  return kFALSE; // Good channel
1264  }
1265  else
1266  {
1267  if ( fBadStatusSelection[0] == kTRUE )
1268  {
1269  return kTRUE; // consider bad hot, dead and warm
1270  }
1271  else
1272  {
1273  if ( fBadStatusSelection[AliCaloCalibPedestal::kDead] == kTRUE &&
1274  status == AliCaloCalibPedestal::kDead )
1275  return kTRUE; // consider bad dead
1276  else if ( fBadStatusSelection[AliCaloCalibPedestal::kHot] == kTRUE &&
1277  status == AliCaloCalibPedestal::kHot )
1278  return kTRUE; // consider bad hot
1279  else if ( fBadStatusSelection[AliCaloCalibPedestal::kWarning] == kTRUE &&
1280  status == AliCaloCalibPedestal::kWarning )
1281  return kTRUE; // consider bad warm
1282  }
1283  }
1284 
1285  AliWarning(Form("Careful, bad channel selection not properly done: ism %d, icol %d, irow %d, status %d,\n"
1286  " fBadAll %d, fBadHot %d, fBadWarm %d, fBadDead %d",
1287  iSM, iCol, iRow, status,
1288  fBadStatusSelection[0], fBadStatusSelection[1],
1289  fBadStatusSelection[2], fBadStatusSelection[3]));
1290 
1291  return kFALSE; // if everything fails, accept it.
1292 }
1293 
1307 //____________________________________________________________________
1308 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1309  AliVCaloCells* cells,
1310  const AliVCluster* clu,
1311  Int_t & absId,
1312  Int_t & iSupMod,
1313  Int_t & ieta,
1314  Int_t & iphi,
1315  Bool_t & shared)
1316 {
1317  Double_t eMax = -1.;
1318  Double_t eCell = -1.;
1319  Float_t fraction = 1.;
1320  Float_t recalFactor = 1.;
1321  Int_t cellAbsId = -1 ;
1322 
1323  Int_t iTower = -1;
1324  Int_t iIphi = -1;
1325  Int_t iIeta = -1;
1326  Int_t iSupMod0= -1;
1327 
1328  if (!clu)
1329  {
1330  AliInfo("Cluster pointer null!");
1331  absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1332  return;
1333  }
1334 
1335  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1336  {
1337  cellAbsId = clu->GetCellAbsId(iDig);
1338  fraction = clu->GetCellAmplitudeFraction(iDig);
1339  //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1340  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1341 
1342  geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1343  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1344 
1345  if (iDig==0)
1346  {
1347  iSupMod0=iSupMod;
1348  } else if (iSupMod0!=iSupMod)
1349  {
1350  shared = kTRUE;
1351  //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1352  }
1353 
1354  if (!fCellsRecalibrated && IsRecalibrationOn())
1355  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1356 
1357  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1358  //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1359  if (eCell > eMax)
1360  {
1361  eMax = eCell;
1362  absId = cellAbsId;
1363  //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1364  }
1365  }// cell loop
1366 
1367  //Get from the absid the supermodule, tower and eta/phi numbers
1368  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1369  //Gives SuperModule and Tower numbers
1370  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1371  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1372  //printf("Max end---\n");
1373 }
1374 
1383 //_________________________________________________________
1384 Float_t AliEMCALRecoUtils::GetCellWeight(Float_t eCell, Float_t eCluster) const
1385 {
1386  if (eCell > 0 && eCluster > 0)
1387  {
1388  if ( fW0 > 0 ) return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ) ) ;
1389  else return TMath::Log( eCluster / eCell ) ;
1390  }
1391  else
1392  return 0. ;
1393 }
1394 
1398 //______________________________________
1399 void AliEMCALRecoUtils::InitParameters()
1400 {
1401  fParticleType = kPhoton;
1402  fPosAlgo = kUnchanged;
1403  fW0 = 4.5;
1404 
1405  fNonLinearityFunction = kNoCorrection;
1406  fNonLinearThreshold = 30;
1407 
1408  fExoticCellFraction = 0.97;
1409  fExoticCellDiffTime = 1e6;
1410  fExoticCellMinAmplitude = 4.0;
1411 
1412  fAODFilterMask = 128;
1413  fAODHybridTracks = kFALSE;
1414  fAODTPCOnlyTracks = kTRUE;
1415 
1416  fCutEtaPhiSum = kTRUE;
1417  fCutEtaPhiSeparate = kFALSE;
1418 
1419  fCutR = 0.05;
1420  fCutEta = 0.025;
1421  fCutPhi = 0.05;
1422 
1423  fClusterWindow = 100;
1424  fMass = 0.139;
1425 
1426  fStepSurface = 20.;
1427  fStepCluster = 5.;
1428  fTrackCutsType = kLooseCut;
1429 
1430  fCutMinTrackPt = 0;
1431  fCutMinNClusterTPC = -1;
1432  fCutMinNClusterITS = -1;
1433 
1434  fCutMaxChi2PerClusterTPC = 1e10;
1435  fCutMaxChi2PerClusterITS = 1e10;
1436 
1437  fCutRequireTPCRefit = kFALSE;
1438  fCutRequireITSRefit = kFALSE;
1439  fCutAcceptKinkDaughters = kFALSE;
1440 
1441  fCutMaxDCAToVertexXY = 1e10;
1442  fCutMaxDCAToVertexZ = 1e10;
1443  fCutDCAToVertex2D = kFALSE;
1444 
1445  fCutRequireITSStandAlone = kFALSE; //MARCEL
1446  fCutRequireITSpureSA = kFALSE; //Marcel
1447 
1448  // Misalignment matrices
1449  for (Int_t i = 0; i < 15 ; i++)
1450  {
1451  fMisalTransShift[i] = 0.;
1452  fMisalRotShift[i] = 0.;
1453  }
1454 
1455  // Non linearity
1456  for (Int_t i = 0; i < 10 ; i++) fNonLinearityParams[i] = 0.;
1457 
1458  // For kBeamTestCorrectedv2 case, but default is no correction
1459  fNonLinearityParams[0] = 0.983504;
1460  fNonLinearityParams[1] = 0.210106;
1461  fNonLinearityParams[2] = 0.897274;
1462  fNonLinearityParams[3] = 0.0829064;
1463  fNonLinearityParams[4] = 152.299;
1464  fNonLinearityParams[5] = 31.5028;
1465  fNonLinearityParams[6] = 0.968;
1466 
1467  // Cluster energy smearing
1468  fSmearClusterEnergy = kFALSE;
1469  fSmearClusterParam[0] = 0.07; // * sqrt E term
1470  fSmearClusterParam[1] = 0.00; // * E term
1471  fSmearClusterParam[2] = 0.00; // constant
1472 }
1473 
1477 //_____________________________________________________
1478 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
1479 {
1480  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1481 
1482  // In order to avoid rewriting the same histograms
1483  Bool_t oldStatus = TH1::AddDirectoryStatus();
1484  TH1::AddDirectory(kFALSE);
1485 
1486  fEMCALRecalibrationFactors = new TObjArray(22);
1487  for (int i = 0; i < 22; i++)
1488  fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1489  Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1490  //Init the histograms with 1
1491  for (Int_t sm = 0; sm < 22; sm++)
1492  {
1493  for (Int_t i = 0; i < 48; i++)
1494  {
1495  for (Int_t j = 0; j < 24; j++)
1496  {
1497  SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1498  }
1499  }
1500  }
1501 
1502  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1503  fEMCALRecalibrationFactors->Compress();
1504 
1505  // In order to avoid rewriting the same histograms
1506  TH1::AddDirectory(oldStatus);
1507 }
1508 
1512 //_________________________________________________________
1513 void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1514 {
1515  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1516 
1517  // In order to avoid rewriting the same histograms
1518  Bool_t oldStatus = TH1::AddDirectoryStatus();
1519  TH1::AddDirectory(kFALSE);
1520 
1521  if(fLowGain) fEMCALTimeRecalibrationFactors = new TObjArray(8);
1522  else fEMCALTimeRecalibrationFactors = new TObjArray(4);
1523 
1524  for (int i = 0; i < 4; i++)
1525  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1526  Form("hAllTimeAvBC%d",i),
1527  48*24*22,0.,48*24*22) );
1528  // Init the histograms with 0
1529  for (Int_t iBC = 0; iBC < 4; iBC++)
1530  {
1531  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1532  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1533  }
1534 
1535  if(fLowGain) {
1536  for (int iBC = 0; iBC < 4; iBC++) {
1537  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1538  Form("hAllTimeAvLGBC%d",iBC),
1539  48*24*22,0.,48*24*22) );
1540  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1541  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1542  }
1543  }
1544 
1545 
1546 
1547  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1548  fEMCALTimeRecalibrationFactors->Compress();
1549 
1550  // In order to avoid rewriting the same histograms
1551  TH1::AddDirectory(oldStatus);
1552 }
1553 
1557 //____________________________________________________
1558 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1559 {
1560  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1561 
1562  // In order to avoid rewriting the same histograms
1563  Bool_t oldStatus = TH1::AddDirectoryStatus();
1564  TH1::AddDirectory(kFALSE);
1565 
1566  fEMCALBadChannelMap = new TObjArray(22);
1567  //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1568 
1569  for (int i = 0; i < 22; i++)
1570  fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1571 
1572  fEMCALBadChannelMap->SetOwner(kTRUE);
1573  fEMCALBadChannelMap->Compress();
1574 
1575  // In order to avoid rewriting the same histograms
1576  TH1::AddDirectory(oldStatus);
1577 }
1578 
1582 //___________________________________________________________
1583 void AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibration()
1584 {
1585  AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
1586 
1587  // In order to avoid rewriting the same histograms
1588  Bool_t oldStatus = TH1::AddDirectoryStatus();
1589  TH1::AddDirectory(kFALSE);
1590 
1591  fEMCALL1PhaseInTimeRecalibration = new TObjArray(1);
1592 
1593  fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
1594 
1595  for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
1596  SetEMCALL1PhaseInTimeRecalibrationForSM(i,0);
1597 
1598  fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
1599  fEMCALL1PhaseInTimeRecalibration->Compress();
1600 
1601  // In order to avoid rewriting the same histograms
1602  TH1::AddDirectory(oldStatus);
1603 }
1604 
1614 //____________________________________________________________________________
1615 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1616  AliVCluster * cluster,
1617  AliVCaloCells * cells,
1618  Int_t bc)
1619 {
1620  if (!cluster)
1621  {
1622  AliInfo("Cluster pointer null!");
1623  return;
1624  }
1625 
1626  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1627  UShort_t * index = cluster->GetCellsAbsId() ;
1628  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1629  Int_t ncells = cluster->GetNCells();
1630 
1631  // Initialize some used variables
1632  Float_t energy = 0;
1633  Int_t absId =-1;
1634  Int_t icol =-1, irow =-1, imod=1;
1635  Float_t factor = 1, frac = 0;
1636  Int_t absIdMax = -1;
1637  Float_t emax = 0;
1638 
1639  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1640  for (Int_t icell = 0; icell < ncells; icell++)
1641  {
1642  absId = index[icell];
1643  frac = fraction[icell];
1644  if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1645 
1646  if (!fCellsRecalibrated && IsRecalibrationOn())
1647  {
1648  // Energy
1649  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1650  geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1651  if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1652  continue;
1653  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1654  factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1655 
1656  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1657  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1658 
1659  }
1660 
1661  energy += cells->GetCellAmplitude(absId)*factor*frac;
1662 
1663  if (emax < cells->GetCellAmplitude(absId)*factor*frac)
1664  {
1665  emax = cells->GetCellAmplitude(absId)*factor*frac;
1666  absIdMax = absId;
1667  }
1668  }
1669 
1670  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1671 
1672  cluster->SetE(energy);
1673 
1674  // Recalculate time of cluster
1675  Double_t timeorg = cluster->GetTOF();
1676  Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
1677 
1678  Double_t time = cells->GetCellTime(absIdMax);
1679  if (!fCellsRecalibrated && IsTimeRecalibrationOn())
1680  RecalibrateCellTime(absIdMax,bc,time,isLowGain);
1681  time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
1682 
1683  // Recalibrate time with L1 phase
1684  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn())
1685  RecalibrateCellTimeL1Phase(imod, bc, time);
1686 
1687  cluster->SetTOF(time);
1688 
1689  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1690 }
1691 
1699 //_______________________________________________________________________
1700 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
1701 {
1702  if (!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn())
1703  return;
1704 
1705  if (!cells)
1706  {
1707  AliInfo("Cells pointer null!");
1708  return;
1709  }
1710 
1711  Short_t absId =-1;
1712  Bool_t accept = kFALSE;
1713  Float_t ecell = 0;
1714  Double_t tcell = 0;
1715  Double_t ecellin = 0;
1716  Double_t tcellin = 0;
1717  Int_t mclabel = -1;
1718  Double_t efrac = 0;
1719 
1720  Int_t nEMcell = cells->GetNumberOfCells() ;
1721  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1722  {
1723  cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1724 
1725  accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1726  if (!accept)
1727  {
1728  ecell = 0;
1729  tcell = -1;
1730  }
1731 
1732  // Set new values
1733  cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1734  }
1735 
1736  fCellsRecalibrated = kTRUE;
1737 }
1738 
1747 //_______________________________________________________________________________________________________
1748 void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
1749 {
1750  if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
1751  if(fLowGain)
1752  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
1753  else
1754  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
1755  }
1756 }
1757 
1765 //_______________________________________________________________________________________________________
1766 void AliEMCALRecoUtils::RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t & celltime) const
1767 {
1768  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
1769  {
1770  bc=bc%4;
1771 
1772  Float_t offsetPerSM=0.;
1773  Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM);
1774  Int_t l1Phase=l1PhaseShift & 3; //bit operation
1775 
1776  if(bc >= l1Phase)
1777  offsetPerSM = (bc - l1Phase)*25;
1778  else
1779  offsetPerSM = (bc - l1Phase + 4)*25;
1780 
1781  Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
1782  l1shiftOffset*=25;
1783 
1784  celltime -= offsetPerSM*1.e-9;
1785  celltime -= l1shiftOffset*1.e-9;
1786  }
1787 }
1788 
1804 //______________________________________________________________________________
1805 void AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator( Int_t absID, AliVCluster* clus, AliMCEvent* mc,
1806  Float_t & amp, TArrayI & labeArr, TArrayF & eDepArr) const
1807 {
1808  TString genName ;
1809  Float_t eDepFrac[4];
1810 
1811  Float_t edepTotFrac = 1;
1812  Bool_t found = kFALSE;
1813  Float_t ampOrg = amp;
1814 
1815  //
1816  // Get the energy deposition fraction from cluster.
1817  //
1818  for(Int_t icluscell = 0; icluscell < clus->GetNCells(); icluscell++ )
1819  {
1820  if(absID == clus->GetCellAbsId(icluscell))
1821  {
1822  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
1823 
1824  found = kTRUE;
1825 
1826  break;
1827  }
1828  }
1829 
1830  if ( !found )
1831  {
1832  AliWarning(Form("Cell abs ID %d NOT found in cluster",absID));
1833  return;
1834  }
1835 
1836  //
1837  // Check if there is a particle contribution from a given generator name.
1838  // If it is not one of the selected generators,
1839  // remove the constribution from the cell.
1840  //
1841  if ( mc && fNMCGenerToAccept > 0 )
1842  {
1843  //printf("Accept contribution from generator? \n");
1844  for(UInt_t imc = 0; imc < 4; imc++)
1845  {
1846  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc )
1847  {
1848  mc->GetCocktailGenerator(clus->GetLabelAt(imc),genName);
1849 
1850  Bool_t generOK = kFALSE;
1851  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
1852  {
1853  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
1854  }
1855 
1856  if ( !generOK )
1857  {
1858  amp-=ampOrg*eDepFrac[imc];
1859 
1860  edepTotFrac-=eDepFrac[imc];
1861 
1862  eDepFrac[imc] = 0;
1863  }
1864 
1865  } // eDep > 0
1866  } // 4 possible loop
1867 
1868  } // accept at least one kind of generator
1869 
1870  //
1871  // Add MC label and Edep to corresponding array (to be used later in digits)
1872  //
1873  Int_t nLabels = 0;
1874  for(UInt_t imc = 0; imc < 4; imc++)
1875  {
1876  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc && edepTotFrac > 0 )
1877  {
1878  nLabels++;
1879 
1880  labeArr.Set(nLabels);
1881  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
1882 
1883  eDepArr.Set(nLabels);
1884  eDepArr.AddAt( (eDepFrac[imc]/edepTotFrac) * amp, nLabels-1);
1885  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
1886  } // edep > 0
1887 
1888  } // mc cell label loop
1889 
1890  //
1891  // If no label found, reject cell
1892  // It can happen to have this case (4 MC labels per cell is not enough for some cases)
1893  // better to remove. To be treated carefully.
1894  //
1895  if ( nLabels == 0 ) amp = 0;
1896 
1897 }
1898 
1906 //______________________________________________________________________________
1907 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1908  AliVCaloCells* cells,
1909  AliVCluster* clu)
1910 {
1911  if (!clu)
1912  {
1913  AliInfo("Cluster pointer null!");
1914  return;
1915  }
1916 
1917  if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1918  else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1919  else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1920 }
1921 
1930 //_____________________________________________________________________________________________
1931 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom,
1932  AliVCaloCells* cells,
1933  AliVCluster* clu)
1934 {
1935  Double_t eCell = 0.;
1936  Float_t fraction = 1.;
1937  Float_t recalFactor = 1.;
1938 
1939  Int_t absId = -1;
1940  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1941  Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1942  Float_t weight = 0., totalWeight=0.;
1943  Float_t newPos[3] = {-1.,-1.,-1.};
1944  Double_t pLocal[3], pGlobal[3];
1945  Bool_t shared = kFALSE;
1946 
1947  Float_t clEnergy = clu->E(); //Energy already recalibrated previously
1948  if (clEnergy <= 0) return;
1949 
1950  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1951  Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1952 
1953  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1954 
1955  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1956  {
1957  absId = clu->GetCellAbsId(iDig);
1958  fraction = clu->GetCellAmplitudeFraction(iDig);
1959  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1960 
1961  if (!fCellsRecalibrated)
1962  {
1963  geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1964  geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1965  if (IsRecalibrationOn()) {
1966  recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1967  }
1968  }
1969 
1970  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1971 
1972  weight = GetCellWeight(eCell,clEnergy);
1973  totalWeight += weight;
1974 
1975  geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1976  //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1977  geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1978  //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1979 
1980  for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1981  }// cell loop
1982 
1983  if (totalWeight>0)
1984  {
1985  for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1986  }
1987 
1988  //Float_t pos[]={0,0,0};
1989  //clu->GetPosition(pos);
1990  //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1991  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1992 
1993  if (iSupModMax > 1) { //sector 1
1994  newPos[0] +=fMisalTransShift[3];//-=3.093;
1995  newPos[1] +=fMisalTransShift[4];//+=6.82;
1996  newPos[2] +=fMisalTransShift[5];//+=1.635;
1997  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1998  } else { //sector 0
1999  newPos[0] +=fMisalTransShift[0];//+=1.134;
2000  newPos[1] +=fMisalTransShift[1];//+=8.2;
2001  newPos[2] +=fMisalTransShift[2];//+=1.197;
2002  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
2003  }
2004  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2005 
2006  clu->SetPosition(newPos);
2007 }
2008 
2017 //____________________________________________________________________________________________
2018 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom,
2019  AliVCaloCells* cells,
2020  AliVCluster* clu)
2021 {
2022  Double_t eCell = 1.;
2023  Float_t fraction = 1.;
2024  Float_t recalFactor = 1.;
2025 
2026  Int_t absId = -1;
2027  Int_t iTower = -1;
2028  Int_t iIphi = -1, iIeta = -1;
2029  Int_t iSupMod = -1, iSupModMax = -1;
2030  Int_t iphi = -1, ieta =-1;
2031  Bool_t shared = kFALSE;
2032 
2033  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
2034 
2035  if (clEnergy <= 0)
2036  return;
2037  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
2038  Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
2039 
2040  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
2041  Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
2042  Int_t startingSM = -1;
2043 
2044  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
2045  {
2046  absId = clu->GetCellAbsId(iDig);
2047  fraction = clu->GetCellAmplitudeFraction(iDig);
2048  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2049 
2050  if (iDig==0) startingSM = iSupMod;
2051  else if (iSupMod != startingSM) areInSameSM = kFALSE;
2052 
2053  eCell = cells->GetCellAmplitude(absId);
2054 
2055  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2056  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2057 
2058  if (!fCellsRecalibrated)
2059  {
2060  if (IsRecalibrationOn()) {
2061  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2062  }
2063  }
2064 
2065  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2066 
2067  weight = GetCellWeight(eCell,clEnergy);
2068  if (weight < 0) weight = 0;
2069  totalWeight += weight;
2070  weightedCol += ieta*weight;
2071  weightedRow += iphi*weight;
2072 
2073  //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
2074  }// cell loop
2075 
2076  Float_t xyzNew[]={-1.,-1.,-1.};
2077  if (areInSameSM == kTRUE)
2078  {
2079  //printf("In Same SM\n");
2080  weightedCol = weightedCol/totalWeight;
2081  weightedRow = weightedRow/totalWeight;
2082  geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2083  }
2084  else
2085  {
2086  //printf("In Different SM\n");
2087  geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2088  }
2089 
2090  clu->SetPosition(xyzNew);
2091 }
2092 
2100 //___________________________________________________________________________________________
2101 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
2102  AliVCaloCells* cells,
2103  AliVCluster * cluster)
2104 {
2105  if (!fRecalDistToBadChannels) return;
2106 
2107  if (!cluster)
2108  {
2109  AliInfo("Cluster pointer null!");
2110  return;
2111  }
2112 
2113  // Get channels map of the supermodule where the cluster is.
2114  Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
2115  Bool_t shared = kFALSE;
2116  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
2117  TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
2118 
2119  Int_t dRrow, dRcol;
2120  Float_t minDist = 10000.;
2121  Float_t dist = 0.;
2122 
2123  // Loop on tower status map
2124  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2125  {
2126  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2127  {
2128  // Check if tower is bad.
2129  if (hMap->GetBinContent(icol,irow)==0) continue;
2130  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
2131  // iSupMod,icol, irow, icolM,irowM);
2132 
2133  dRrow=TMath::Abs(irowM-irow);
2134  dRcol=TMath::Abs(icolM-icol);
2135  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2136  if (dist < minDist)
2137  {
2138  //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
2139  minDist = dist;
2140  }
2141  }
2142  }
2143 
2144  // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
2145  if (shared)
2146  {
2147  TH2D* hMap2 = 0;
2148  Int_t iSupMod2 = -1;
2149 
2150  // The only possible combinations are (0,1), (2,3) ... (8,9)
2151  if (iSupMod%2) iSupMod2 = iSupMod-1;
2152  else iSupMod2 = iSupMod+1;
2153  hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
2154 
2155  // Loop on tower status map of second super module
2156  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2157  {
2158  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2159  {
2160  // Check if tower is bad.
2161  if (hMap2->GetBinContent(icol,irow)==0) continue;
2162 
2163  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
2164  // iSupMod2,icol, irow,iSupMod,icolM,irowM);
2165  dRrow=TMath::Abs(irow-irowM);
2166 
2167  if (iSupMod%2)
2168  dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
2169  else
2170  dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
2171 
2172  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2173  if (dist < minDist) minDist = dist;
2174  }
2175  }
2176  }// shared cluster in 2 SuperModules
2177 
2178  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
2179  cluster->SetDistanceToBadChannel(minDist);
2180 }
2181 
2188 //__________________________________________________________________
2189 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
2190 {
2191  if (!cluster)
2192  {
2193  AliInfo("Cluster pointer null!");
2194  return;
2195  }
2196 
2197  if (cluster->GetM02() != 0)
2198  fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
2199 
2200  Float_t pidlist[AliPID::kSPECIESCN+1];
2201  for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
2202 
2203  cluster->SetPID(pidlist);
2204 }
2205 
2226 //___________________________________________________________________________________________________________________
2227 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry * geom,
2228  AliVCaloCells* cells, AliVCluster * cluster,
2229  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2230  Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
2231  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2232  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2233 {
2234  if (!cluster)
2235  {
2236  AliInfo("Cluster pointer null!");
2237  return;
2238  }
2239 
2240  Double_t eCell = 0.;
2241  Double_t tCell = 0.;
2242  Float_t fraction = 1.;
2243  Float_t recalFactor = 1.;
2244  Bool_t isLowGain = kFALSE;
2245 
2246  Int_t iSupMod = -1;
2247  Int_t iTower = -1;
2248  Int_t iIphi = -1;
2249  Int_t iIeta = -1;
2250  Int_t iphi = -1;
2251  Int_t ieta = -1;
2252  Double_t etai = -1.;
2253  Double_t phii = -1.;
2254 
2255  Int_t nstat = 0 ;
2256  Float_t wtot = 0.;
2257  Double_t w = 0.;
2258  Double_t etaMean = 0.;
2259  Double_t phiMean = 0.;
2260 
2261  Double_t pGlobal[3];
2262 
2263  // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
2264  // or the non linearity correction was applied
2265  // and to check if the cluster is between 2 SM in eta
2266  Int_t iSM0 = -1;
2267  Bool_t shared = kFALSE;
2268  Float_t energy = 0;
2269 
2270  enAfterCuts = 0;
2271  l0 = 0; l1 = 0;
2272  disp = 0; dEta = 0; dPhi = 0;
2273  sEta = 0; sPhi = 0; sEtaPhi = 0;
2274 
2275  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2276  {
2277  // Get from the absid the supermodule, tower and eta/phi numbers
2278 
2279  Int_t absId = cluster->GetCellAbsId(iDigit);
2280 
2281  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2282  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2283 
2284  // Check if there are cells of different SM
2285  if (iDigit == 0 ) iSM0 = iSupMod;
2286  else if (iSupMod!= iSM0) shared = kTRUE;
2287 
2288  // Get the cell energy, if recalibration is on, apply factors
2289  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2290  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2291 
2292  if (IsRecalibrationOn())
2293  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2294 
2295  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2296  tCell = cells->GetCellTime (absId);
2297  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2298 
2299  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2300  tCell*=1e9;
2301  tCell-=fConstantTimeShift;
2302 
2303  if(eCell > 0.05) // at least a minimum 50 MeV cut
2304  energy += eCell;
2305 
2306  if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
2307  enAfterCuts += eCell;
2308  } // cell loop
2309 
2310 
2311  // Loop on cells to calculate weights and shower shape terms parameters
2312  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2313  {
2314  // Get from the absid the supermodule, tower and eta/phi numbers
2315  Int_t absId = cluster->GetCellAbsId(iDigit);
2316 
2317  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2318  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2319 
2320  //Get the cell energy, if recalibration is on, apply factors
2321  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2322  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2323 
2324  if (!fCellsRecalibrated && IsRecalibrationOn())
2325  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2326 
2327  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2328  tCell = cells->GetCellTime (absId);
2329  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2330 
2331  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2332  tCell*=1e9;
2333  tCell-=fConstantTimeShift;
2334 
2335  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2336  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2337  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2338 
2339  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2340  {
2341  w = GetCellWeight(eCell, energy);
2342 
2343  // Cell index
2344  if ( fShowerShapeCellLocationType == 0 )
2345  {
2346  etai=(Double_t)ieta;
2347  phii=(Double_t)iphi;
2348  }
2349  // Cell angle location
2350  else if( fShowerShapeCellLocationType == 1 )
2351  {
2352  geom->EtaPhiFromIndex(absId, etai, phii);
2353  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2354  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2355  }
2356  else
2357  {
2358  geom->GetGlobal(absId,pGlobal);
2359 
2360  // Cell x-z location
2361  if( fShowerShapeCellLocationType == 2 )
2362  {
2363  etai = pGlobal[2];
2364  phii = pGlobal[0];
2365  }
2366  // Cell r-z location
2367  else
2368  {
2369  etai = pGlobal[2];
2370  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2371  }
2372  }
2373 
2374  if (w > 0.0)
2375  {
2376  wtot += w ;
2377  nstat++;
2378 
2379  // Shower shape
2380  sEta += w * etai * etai ;
2381  etaMean += w * etai ;
2382  sPhi += w * phii * phii ;
2383  phiMean += w * phii ;
2384  sEtaPhi += w * etai * phii ;
2385  }
2386  }
2387  else if(eCell > 0.05)
2388  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2389  } // cell loop
2390 
2391  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot);
2392 
2393  // Normalize to the weight
2394  if (wtot > 0)
2395  {
2396  etaMean /= wtot ;
2397  phiMean /= wtot ;
2398  }
2399  else
2400  AliDebug(2,Form("Wrong weight %f\n", wtot));
2401 
2402  // Loop on cells to calculate dispersion
2403  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2404  {
2405  // Get from the absid the supermodule, tower and eta/phi numbers
2406  Int_t absId = cluster->GetCellAbsId(iDigit);
2407 
2408  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2409  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2410 
2411  //Get the cell energy, if recalibration is on, apply factors
2412  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2413  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2414  if (IsRecalibrationOn())
2415  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2416 
2417  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2418  tCell = cells->GetCellTime (absId);
2419  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2420 
2421  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2422  tCell*=1e9;
2423  tCell-=fConstantTimeShift;
2424 
2425  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2426  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2427  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2428 
2429  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2430  {
2431  w = GetCellWeight(eCell,cluster->E());
2432 
2433  // Cell index
2434  if ( fShowerShapeCellLocationType == 0 )
2435  {
2436  etai=(Double_t)ieta;
2437  phii=(Double_t)iphi;
2438  }
2439  // Cell angle location
2440  else if( fShowerShapeCellLocationType == 1 )
2441  {
2442  geom->EtaPhiFromIndex(absId, etai, phii);
2443  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2444  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2445  }
2446  else
2447  {
2448  geom->GetGlobal(absId,pGlobal);
2449 
2450  // Cell x-z location
2451  if( fShowerShapeCellLocationType == 2 )
2452  {
2453  etai = pGlobal[2];
2454  phii = pGlobal[0];
2455  }
2456  // Cell r-z location
2457  else
2458  {
2459  etai = pGlobal[2];
2460  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2461  }
2462  }
2463 
2464  if (w > 0.0)
2465  {
2466  disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2467  dEta += w * (etai-etaMean)*(etai-etaMean) ;
2468  dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2469  }
2470  }
2471  else
2472  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2473  }// cell loop
2474 
2475  // Normalize to the weigth and set shower shape parameters
2476  if (wtot > 0 && nstat > 1)
2477  {
2478  disp /= wtot ;
2479  dEta /= wtot ;
2480  dPhi /= wtot ;
2481  sEta /= wtot ;
2482  sPhi /= wtot ;
2483  sEtaPhi /= wtot ;
2484 
2485  sEta -= etaMean * etaMean ;
2486  sPhi -= phiMean * phiMean ;
2487  sEtaPhi -= etaMean * phiMean ;
2488 
2489  l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2490  l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2491 
2492  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f l0 %f l1 %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot,l0,l1);
2493 
2494  }
2495  else
2496  {
2497  l0 = 0. ;
2498  l1 = 0. ;
2499  dEta = 0. ; dPhi = 0. ; disp = 0. ;
2500  sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2501  }
2502 }
2503 
2522 //___________________________________________________________________________________________________________________
2523 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
2524  AliVCaloCells* cells, AliVCluster * cluster,
2525  Float_t & l0, Float_t & l1,
2526  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2527  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2528 {
2529  Float_t newEnergy = 0;
2530  Float_t cellEmin = 0.05; // 50 MeV
2531  Float_t cellTimeWindow = 1000; // open cut
2532  Int_t bc = 0;
2533  AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(geom, cells, cluster,
2534  cellEmin, cellTimeWindow, bc,
2535  newEnergy, l0, l1, disp,
2536  dEta, dPhi, sEta, sPhi, sEtaPhi);
2537 }
2538 
2547 //____________________________________________________________________________________________
2548 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
2549  AliVCaloCells* cells,
2550  AliVCluster * cluster)
2551 {
2552  Float_t l0 = 0., l1 = 0.;
2553  Float_t disp = 0., dEta = 0., dPhi = 0.;
2554  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2555 
2556  AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
2557  dEta, dPhi, sEta, sPhi, sEtaPhi);
2558 
2559  cluster->SetM02(l0);
2560  cluster->SetM20(l1);
2561  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2562 
2563 }
2564 
2577 //____________________________________________________________________________________________
2578 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry * geom,
2579  AliVCaloCells* cells, AliVCluster * cluster,
2580  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2581  Float_t & enAfterCuts)
2582 {
2583  Float_t l0 = 0., l1 = 0.;
2584  Float_t disp = 0., dEta = 0., dPhi = 0.;
2585  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2586 
2587  AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(geom, cells, cluster,
2588  cellEcut, cellTimeCut, bc,
2589  enAfterCuts, l0, l1, disp,
2590  dEta, dPhi, sEta, sPhi, sEtaPhi);
2591 
2592  cluster->SetM02(l0);
2593  cluster->SetM20(l1);
2594  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2595 
2596 }
2597 
2611 //____________________________________________________________________________
2612 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
2613  TObjArray * clusterArr,
2614  const AliEMCALGeometry *geom,
2615  AliMCEvent * mc)
2616 {
2617  fMatchedTrackIndex ->Reset();
2618  fMatchedClusterIndex->Reset();
2619  fResidualPhi->Reset();
2620  fResidualEta->Reset();
2621 
2622  fMatchedTrackIndex ->Set(1000);
2623  fMatchedClusterIndex->Set(1000);
2624  fResidualPhi->Set(1000);
2625  fResidualEta->Set(1000);
2626 
2627  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
2628  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
2629 
2630  // Init the magnetic field if not already on
2631  if (!TGeoGlobalMagField::Instance()->GetField())
2632  {
2633  if (!event->InitMagneticField())
2634  {
2635  AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
2636  }
2637  } // Init mag field
2638 
2639  if (esdevent)
2640  {
2641  UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
2642  UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
2643  Bool_t desc1 = (mask1 >> 3) & 0x1;
2644  Bool_t desc2 = (mask2 >> 3) & 0x1;
2645  if (desc1==0 || desc2==0) {
2646  // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
2647  // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
2648  // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
2649  fITSTrackSA=kTRUE;
2650  }
2651  }
2652 
2653  TObjArray *clusterArray = 0x0;
2654  if (!clusterArr)
2655  {
2656  clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
2657  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2658  {
2659  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2660  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
2661  continue;
2662  clusterArray->AddAt(cluster,icl);
2663  }
2664  }
2665 
2666  Int_t matched=0;
2667  Double_t cv[21];
2668  TString genName;
2669  for (Int_t i=0; i<21;i++) cv[i]=0;
2670  for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
2671  {
2672  AliExternalTrackParam *trackParam = 0;
2673  Int_t mcLabel = -1;
2674  // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
2675  AliESDtrack *esdTrack = 0;
2676  AliAODTrack *aodTrack = 0;
2677  if (esdevent)
2678  {
2679  esdTrack = esdevent->GetTrack(itr);
2680  if (!esdTrack) continue;
2681  if (!IsAccepted(esdTrack)) continue;
2682  if (esdTrack->Pt()<fCutMinTrackPt) continue;
2683 
2684  if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
2685 
2686  // Save some time and memory in case of no DCal present
2687  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2688  {
2689  Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
2690  if ( phi <= 10 || phi >= 250 ) continue;
2691  }
2692 
2693  if (!fITSTrackSA)
2694  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
2695  else
2696  trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
2697 
2698  mcLabel = TMath::Abs(esdTrack->GetLabel());
2699  }
2700 
2701  // If the input event is AOD, the starting point for extrapolation is at vertex
2702  // AOD tracks are selected according to its filterbit.
2703  else if (aodevent)
2704  {
2705  aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
2706 
2707  if (!aodTrack) AliFatal("Not a standard AOD");
2708 
2709  if (!aodTrack) continue;
2710 
2711  if (fAODTPCOnlyTracks)
2712  { // Match with TPC only tracks, default from May 2013, before filter bit 32
2713  //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
2714  if (!aodTrack->IsTPCConstrained()) continue ;
2715  }
2716  else if (fAODHybridTracks)
2717  { // Match with hybrid tracks
2718  //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
2719  if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
2720  } else
2721  { // Match with tracks on a mask
2722  //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
2723  if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
2724  }
2725 
2726  if (aodTrack->Pt() < fCutMinTrackPt) continue;
2727 
2728  if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
2729 
2730  // Save some time and memory in case of no DCal present
2731  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2732  {
2733  Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
2734  if ( phi <= 10 || phi >= 250 ) continue;
2735  }
2736 
2737  Double_t pos[3],mom[3];
2738  aodTrack->GetXYZ(pos);
2739  aodTrack->GetPxPyPz(mom);
2740  AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
2741  itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
2742 
2743  trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
2744 
2745  mcLabel = TMath::Abs(aodTrack->GetLabel());
2746  }
2747 
2748  //Return if the input data is not "AOD" or "ESD"
2749  else
2750  {
2751  AliWarning("Wrong input data type! Should be \"AOD\" or \"ESD\" ");
2752  if (clusterArray)
2753  {
2754  clusterArray->Clear();
2755  delete clusterArray;
2756  }
2757  return;
2758  }
2759 
2760  if (!trackParam) continue;
2761 
2762  //
2763  // Check if track comes from a particular MC generator, do not include it
2764  // if it is not a selected one
2765  //
2766  if( mc && fMCGenerToAcceptForTrack && fNMCGenerToAccept > 0 )
2767  {
2768  mc->GetCocktailGenerator(mcLabel,genName);
2769 
2770  Bool_t generOK = kFALSE;
2771  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
2772  {
2773  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
2774  }
2775 
2776  if ( !generOK ) continue;
2777  }
2778 
2779  // Extrapolate the track to EMCal surface, see AliEMCALRecoUtilsBase
2780  AliExternalTrackParam emcalParam(*trackParam);
2781  Float_t eta, phi, pt;
2782  if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2783  {
2784  if (aodevent && trackParam) delete trackParam;
2785  if (fITSTrackSA && trackParam) delete trackParam;
2786  continue;
2787  }
2788 
2789  if ( TMath::Abs(eta) > 0.75 )
2790  {
2791  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2792  continue;
2793  }
2794 
2795  // Save some time and memory in case of no DCal present
2796  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2797  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2798  {
2799  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2800  continue;
2801  }
2802 
2803  //Find matched clusters
2804  Int_t index = -1;
2805  Float_t dEta = -999, dPhi = -999;
2806  if (!clusterArr)
2807  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
2808  else
2809  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
2810 
2811 
2812  if (index>-1)
2813  {
2814  fMatchedTrackIndex ->AddAt(itr,matched);
2815  fMatchedClusterIndex ->AddAt(index,matched);
2816  fResidualEta ->AddAt(dEta,matched);
2817  fResidualPhi ->AddAt(dPhi,matched);
2818  matched++;
2819  }
2820 
2821  if (aodevent && trackParam) delete trackParam;
2822  if (fITSTrackSA && trackParam) delete trackParam;
2823  }//track loop
2824 
2825  if (clusterArray)
2826  {
2827  clusterArray->Clear();
2828  delete clusterArray;
2829  }
2830 
2831  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2832 
2833  fMatchedTrackIndex ->Set(matched);
2834  fMatchedClusterIndex ->Set(matched);
2835  fResidualPhi ->Set(matched);
2836  fResidualEta ->Set(matched);
2837 }
2838 
2850 //________________________________________________________________________________
2851 Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
2852  const AliVEvent *event,
2853  const AliEMCALGeometry *geom,
2854  Float_t &dEta, Float_t &dPhi)
2855 {
2856  Int_t index = -1;
2857 
2858  if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
2859 
2860  // Save some time and memory in case of no DCal present
2861  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2862  {
2863  Double_t phiV = track->Phi()*TMath::RadToDeg();
2864  if ( phiV <= 10 || phiV >= 250 ) return index;
2865  }
2866 
2867  AliExternalTrackParam *trackParam = 0;
2868  if (!fITSTrackSA)
2869  trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
2870  else
2871  trackParam = new AliExternalTrackParam(*track);
2872 
2873  if (!trackParam) return index;
2874  AliExternalTrackParam emcalParam(*trackParam);
2875 
2876  Float_t eta, phi, pt;
2877  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2878  {
2879  if (fITSTrackSA) delete trackParam;
2880  return index;
2881  }
2882 
2883  if ( TMath::Abs(eta) > 0.75 )
2884  {
2885  if (fITSTrackSA) delete trackParam;
2886  return index;
2887  }
2888 
2889  // Save some time and memory in case of no DCal present
2890  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2891  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2892  {
2893  if (fITSTrackSA) delete trackParam;
2894  return index;
2895  }
2896 
2897  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2898 
2899  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2900  {
2901  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2902  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2903  clusterArr->AddAt(cluster,icl);
2904  }
2905 
2906  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2907  clusterArr->Clear();
2908  delete clusterArr;
2909  if (fITSTrackSA) delete trackParam;
2910 
2911  return index;
2912 }
2913 
2924 //_______________________________________________________________________________________________
2925 Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2926  AliExternalTrackParam *trkParam,
2927  const TObjArray * clusterArr,
2928  Float_t &dEta, Float_t &dPhi)
2929 {
2930  dEta=-999, dPhi=-999;
2931  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2932  Int_t index = -1;
2933  Float_t tmpEta=-999, tmpPhi=-999;
2934 
2935  Double_t exPos[3] = {0.,0.,0.};
2936  if (!emcalParam->GetXYZ(exPos)) return index;
2937 
2938  Float_t clsPos[3] = {0.,0.,0.};
2939  for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2940  {
2941  AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2942 
2943  if (!cluster || !cluster->IsEMCAL()) continue;
2944 
2945  cluster->GetPosition(clsPos);
2946 
2947  Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
2948  if (dR > fClusterWindow) continue;
2949 
2950  AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
2951 
2952  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2953 
2954  if (fCutEtaPhiSum)
2955  {
2956  Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2957  if (tmpR<dRMax)
2958  {
2959  dRMax=tmpR;
2960  dEtaMax=tmpEta;
2961  dPhiMax=tmpPhi;
2962  index=icl;
2963  }
2964  }
2965  else if (fCutEtaPhiSeparate)
2966  {
2967  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
2968  {
2969  dEtaMax = tmpEta;
2970  dPhiMax = tmpPhi;
2971  index=icl;
2972  }
2973  }
2974  else
2975  {
2976  AliError("Please specify your cut criteria");
2977  AliError("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()");
2978  AliError("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()");
2979  return index;
2980  }
2981  }
2982 
2983  dEta=dEtaMax;
2984  dPhi=dPhiMax;
2985 
2986  return index;
2987 }
2988 
3001 //---------------------------------------------------------------------------------
3002 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
3003  const AliVCluster *cluster,
3004  Float_t &tmpEta,
3005  Float_t &tmpPhi)
3006 {
3007  return AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
3008 }
3009 
3019 //_______________________________________________________________________
3020 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex,
3021  Float_t &dEta, Float_t &dPhi)
3022 {
3023  if (FindMatchedPosForCluster(clsIndex) >= 999)
3024  {
3025  AliDebug(2,"No matched tracks found!\n");
3026  dEta=999.;
3027  dPhi=999.;
3028  return;
3029  }
3030 
3031  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
3032  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
3033 }
3034 
3044 //______________________________________________________________________________________________
3045 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
3046 {
3047  if (FindMatchedPosForTrack(trkIndex) >= 999)
3048  {
3049  AliDebug(2,"No matched cluster found!\n");
3050  dEta=999.;
3051  dPhi=999.;
3052  return;
3053  }
3054 
3055  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
3056  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
3057 }
3058 
3065 //__________________________________________________________
3066 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
3067 {
3068  if (IsClusterMatched(clsIndex))
3069  return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
3070  else
3071  return -1;
3072 }
3073 
3080 //__________________________________________________________
3081 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
3082 {
3083  if (IsTrackMatched(trkIndex))
3084  return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
3085  else
3086  return -1;
3087 }
3088 
3096 //______________________________________________________________
3097 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
3098 {
3099  if (FindMatchedPosForCluster(clsIndex) < 999)
3100  return kTRUE;
3101  else
3102  return kFALSE;
3103 }
3104 
3112 //____________________________________________________________
3113 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
3114 {
3115  if (FindMatchedPosForTrack(trkIndex) < 999)
3116  return kTRUE;
3117  else
3118  return kFALSE;
3119 }
3120 
3128 //______________________________________________________________________
3129 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
3130 {
3131  Float_t tmpR = fCutR;
3132  UInt_t pos = 999;
3133 
3134  for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
3135  {
3136  if (fMatchedClusterIndex->At(i)==clsIndex)
3137  {
3138  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3139  if (r<tmpR)
3140  {
3141  pos=i;
3142  tmpR=r;
3143  AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3144  fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3145  }
3146  }
3147  }
3148  return pos;
3149 }
3150 
3158 //____________________________________________________________________
3159 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
3160 {
3161  Float_t tmpR = fCutR;
3162  UInt_t pos = 999;
3163 
3164  for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3165  {
3166  if (fMatchedTrackIndex->At(i)==trkIndex)
3167  {
3168  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3169  if (r<tmpR)
3170  {
3171  pos=i;
3172  tmpR=r;
3173  AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3174  fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3175  }
3176  }
3177  }
3178  return pos;
3179 }
3180 
3195 //__________________________________________________________________________
3196 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
3197  const AliEMCALGeometry *geom,
3198  AliVCaloCells* cells, Int_t bc)
3199 {
3200  Bool_t isGood=kTRUE;
3201 
3202  if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3203  if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3204  if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3205  if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3206 
3207  return isGood;
3208 }
3209 
3220 //__________________________________________________________
3221 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
3222 {
3223  UInt_t status = esdTrack->GetStatus();
3224 
3225  Int_t nClustersITS = esdTrack->GetITSclusters(0);
3226  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3227 
3228  Float_t chi2PerClusterITS = -1;
3229  Float_t chi2PerClusterTPC = -1;
3230  if (nClustersITS!=0)
3231  chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3232  if (nClustersTPC!=0)
3233  chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3234 
3235  //
3236  // DCA cuts
3237  // Only to be used for primary
3238  //
3239  if ( fTrackCutsType == kGlobalCut )
3240  {
3241  Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3242  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3243 
3244  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3245 
3246  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3247  }
3248  else if( fTrackCutsType == kGlobalCut2011 )
3249  {
3250  Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3251  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3252 
3253  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3254 
3255  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3256  }
3257 
3258  Float_t b[2];
3259  Float_t bCov[3];
3260  esdTrack->GetImpactParameters(b,bCov);
3261  if (bCov[0]<=0 || bCov[2]<=0)
3262  {
3263  AliDebug(1, "Estimated b resolution lower or equal zero!");
3264  bCov[0]=0; bCov[2]=0;
3265  }
3266 
3267  Float_t dcaToVertexXY = b[0];
3268  Float_t dcaToVertexZ = b[1];
3269  Float_t dcaToVertex = -1;
3270 
3271  if (fCutDCAToVertex2D)
3272  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3273  dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3274  else
3275  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3276 
3277  // cut the track?
3278 
3279  Bool_t cuts[kNCuts];
3280  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3281 
3282  // track quality cuts
3283  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3284  cuts[0]=kTRUE;
3285  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3286  cuts[1]=kTRUE;
3287  if (nClustersTPC<fCutMinNClusterTPC)
3288  cuts[2]=kTRUE;
3289  if (nClustersITS<fCutMinNClusterITS)
3290  cuts[3]=kTRUE;
3291  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3292  cuts[4]=kTRUE;
3293  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3294  cuts[5]=kTRUE;
3295  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3296  cuts[6]=kTRUE;
3297  if (fCutDCAToVertex2D && dcaToVertex > 1)
3298  cuts[7] = kTRUE;
3299  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3300  cuts[8] = kTRUE;
3301  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3302  cuts[9] = kTRUE;
3303 
3304  if (fTrackCutsType == kGlobalCut || fTrackCutsType == kGlobalCut2011)
3305  {
3306  //Require at least one SPD point + anything else in ITS
3307  if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3308  cuts[10] = kTRUE;
3309  }
3310 
3311  // ITS
3312  if (fCutRequireITSStandAlone || fCutRequireITSpureSA)
3313  {
3314  if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3315  {
3316  // TPC tracks
3317  cuts[11] = kTRUE;
3318  }
3319  else
3320  {
3321  // ITS standalone tracks
3322  if (fCutRequireITSStandAlone && !fCutRequireITSpureSA)
3323  {
3324  if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3325  }
3326  else if (fCutRequireITSpureSA)
3327  {
3328  if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3329  }
3330  }
3331  }
3332 
3333  Bool_t cut=kFALSE;
3334  for (Int_t i=0; i<kNCuts; i++)
3335  if (cuts[i]) { cut = kTRUE ; }
3336 
3337  // cut the track
3338  if (cut)
3339  return kFALSE;
3340  else
3341  return kTRUE;
3342 }
3343 
3349 //_____________________________________
3350 void AliEMCALRecoUtils::InitTrackCuts()
3351 {
3352  switch (fTrackCutsType)
3353  {
3354  case kTPCOnlyCut:
3355  {
3356  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3357  //TPC
3358  SetMinNClustersTPC(70);
3359  SetMaxChi2PerClusterTPC(4);
3360  SetAcceptKinkDaughters(kFALSE);
3361  SetRequireTPCRefit(kFALSE);
3362 
3363  //ITS
3364  SetRequireITSRefit(kFALSE);
3365  SetMaxDCAToVertexZ(3.2);
3366  SetMaxDCAToVertexXY(2.4);
3367  SetDCAToVertex2D(kTRUE);
3368 
3369  break;
3370  }
3371 
3372  case kGlobalCut:
3373  {
3374  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3375  //TPC
3376  SetMinNClustersTPC(70);
3377  SetMaxChi2PerClusterTPC(4);
3378  SetAcceptKinkDaughters(kFALSE);
3379  SetRequireTPCRefit(kTRUE);
3380 
3381  //ITS
3382  SetRequireITSRefit(kTRUE);
3383  SetMaxDCAToVertexZ(2);
3384  SetMaxDCAToVertexXY();
3385  SetDCAToVertex2D(kFALSE);
3386 
3387  break;
3388  }
3389 
3390  case kLooseCut:
3391  {
3392  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3393  SetMinNClustersTPC(50);
3394  SetAcceptKinkDaughters(kTRUE);
3395 
3396  break;
3397  }
3398 
3399  case kITSStandAlone:
3400  {
3401  AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3402  SetRequireITSRefit(kTRUE);
3403  SetRequireITSStandAlone(kTRUE);
3404  SetITSTrackSA(kTRUE);
3405  break;
3406  }
3407 
3408  case kGlobalCut2011:
3409  {
3410  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3411  //TPC
3412  SetMinNClustersTPC(50);
3413  SetMaxChi2PerClusterTPC(4);
3414  SetAcceptKinkDaughters(kFALSE);
3415  SetRequireTPCRefit(kTRUE);
3416 
3417  //ITS
3418  SetRequireITSRefit(kTRUE);
3419  SetMaxDCAToVertexZ(2);
3420  SetMaxDCAToVertexXY();
3421  SetDCAToVertex2D(kFALSE);
3422 
3423  break;
3424  }
3425 
3426  case kLooseCutWithITSrefit:
3427  {
3428  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3429  SetMinNClustersTPC(50);
3430  SetAcceptKinkDaughters(kTRUE);
3431  SetRequireITSRefit(kTRUE);
3432 
3433  break;
3434  }
3435  }
3436 }
3437 
3443 //________________________________________________________________________
3444 void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
3445 {
3446  Int_t nTracks = event->GetNumberOfTracks();
3447  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3448  {
3449  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3450  if (!track)
3451  {
3452  AliWarning(Form("Could not receive track %d", iTrack));
3453  continue;
3454  }
3455 
3456  Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3457  track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3458 
3459  /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3460  AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3461  if (esdtrack)
3462  {
3463  if (matchClusIndex != -1)
3464  esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3465  else
3466  esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3467  }
3468  else
3469  {
3470  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3471  if (matchClusIndex != -1)
3472  aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3473  else
3474  aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3475  }
3476  }
3477  AliDebug(2,"Track matched to closest cluster");
3478 }
3479 
3486 //_________________________________________________________________________
3487 void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
3488 {
3489  for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
3490  {
3491  AliVCluster *cluster = event->GetCaloCluster(iClus);
3492  if (!cluster->IsEMCAL())
3493  continue;
3494 
3495  //
3496  // Remove old matches in cluster
3497  //
3498  if(cluster->GetNTracksMatched() > 0)
3499  {
3500  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
3501  {
3502  TArrayI arrayTrackMatched(0);
3503  ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
3504  }
3505  else
3506  {
3507  for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
3508  {
3509  ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
3510  }
3511  }
3512  }
3513 
3514  //
3515  // Find new matches and put them in the cluster
3516  //
3517  Int_t nTracks = event->GetNumberOfTracks();
3518  TArrayI arrayTrackMatched(nTracks);
3519 
3520  // Get the closest track matched to the cluster
3521  Int_t nMatched = 0;
3522  Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
3523  if (matchTrackIndex != -1)
3524  {
3525  arrayTrackMatched[nMatched] = matchTrackIndex;
3526  nMatched++;
3527  }
3528 
3529  // Get all other tracks matched to the cluster
3530  for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
3531  {
3532  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
3533 
3534  if( !track ) continue;
3535 
3536  if ( iTrk == matchTrackIndex ) continue;
3537 
3538  if ( track->GetEMCALcluster() == iClus )
3539  {
3540  arrayTrackMatched[nMatched] = iTrk;
3541  ++nMatched;
3542  }
3543  }
3544 
3545  AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
3546 
3547  arrayTrackMatched.Set(nMatched);
3548  AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
3549  if (esdcluster)
3550  esdcluster->AddTracksMatched(arrayTrackMatched);
3551  else if ( nMatched > 0 )
3552  {
3553  AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
3554  if (aodcluster)
3555  {
3556  aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
3557  //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
3558  //printf("Is the closest matching track with ID %d a 128? %d what's its full filter map? %u\n",aodtrack->GetID(), aodtrack->TestFilterBit(128),aodtrack->GetFilterMap());
3559  //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
3560  }
3561  }
3562 
3563  Float_t eta= -999, phi = -999;
3564  if (matchTrackIndex != -1)
3565  GetMatchedResiduals(iClus, eta, phi);
3566 
3567  cluster->SetTrackDistance(phi, eta);
3568  }
3569 
3570  AliDebug(2,"Cluster matched to tracks");
3571 }
3572 
3573 void AliEMCALRecoUtils::SetEMCALChannelRecalibrationFactors(const TObjArray *map) {
3574  if(fEMCALRecalibrationFactors) fEMCALRecalibrationFactors->Clear();
3575  else {
3576  fEMCALRecalibrationFactors = new TObjArray(map->GetEntries());
3577  fEMCALRecalibrationFactors->SetOwner(true);
3578  }
3579  if(!fEMCALRecalibrationFactors->IsOwner()){
3580  // Must claim ownership since the new objects are owend by this instance
3581  fEMCALRecalibrationFactors->SetOwner(kTRUE);
3582  }
3583  for(int i = 0; i < map->GetEntries(); i++){
3584  TH2F *hist = dynamic_cast<TH2F *>(map->At(i));
3585  if(!hist) continue;
3586  this->SetEMCALChannelRecalibrationFactors(i, hist);
3587  }
3588 }
3589 
3590 void AliEMCALRecoUtils::SetEMCALChannelRecalibrationFactors(Int_t iSM , const TH2F* h) {
3591  if(!fEMCALRecalibrationFactors){
3592  fEMCALRecalibrationFactors = new TObjArray(iSM);
3593  fEMCALRecalibrationFactors->SetOwner(true);
3594  }
3595  if(fEMCALRecalibrationFactors->GetEntries() <= iSM) fEMCALRecalibrationFactors->Expand(iSM+1);
3596  if(fEMCALRecalibrationFactors->At(iSM)) fEMCALRecalibrationFactors->RemoveAt(iSM);
3597  TH2F *clone = new TH2F(*h);
3598  clone->SetDirectory(NULL);
3599  fEMCALRecalibrationFactors->AddAt(clone,iSM);
3600 }
3601 
3602 void AliEMCALRecoUtils::SetEMCALChannelStatusMap(const TObjArray *map) {
3603  if(fEMCALBadChannelMap) fEMCALBadChannelMap->Clear();
3604  else {
3605  fEMCALBadChannelMap = new TObjArray(map->GetEntries());
3606  fEMCALBadChannelMap->SetOwner(true);
3607  }
3608  if(!fEMCALBadChannelMap->IsOwner()) {
3609  // Must claim ownership since the new objects are owend by this instance
3610  fEMCALBadChannelMap->SetOwner(true);
3611  }
3612  for(int i = 0; i < map->GetEntries(); i++){
3613  TH2I *hist = dynamic_cast<TH2I *>(map->At(i));
3614  if(!hist) continue;
3615  this->SetEMCALChannelStatusMap(i, hist);
3616  }
3617 }
3618 
3619 void AliEMCALRecoUtils::SetEMCALChannelStatusMap(Int_t iSM , const TH2I* h) {
3620  if(!fEMCALBadChannelMap){
3621  fEMCALBadChannelMap = new TObjArray(iSM);
3622  fEMCALBadChannelMap->SetOwner(true);
3623  }
3624  if(fEMCALBadChannelMap->GetEntries() <= iSM) fEMCALBadChannelMap->Expand(iSM+1);
3625  if(fEMCALBadChannelMap->At(iSM)) fEMCALBadChannelMap->RemoveAt(iSM);
3626  TH2I *clone = new TH2I(*h);
3627  clone->SetDirectory(NULL);
3628  fEMCALBadChannelMap->AddAt(clone,iSM);
3629 }
3630 
3631 void AliEMCALRecoUtils::SetEMCALChannelTimeRecalibrationFactors(const TObjArray *map) {
3632  if(fEMCALTimeRecalibrationFactors) fEMCALTimeRecalibrationFactors->Clear();
3633  else {
3634  fEMCALTimeRecalibrationFactors = new TObjArray(map->GetEntries());
3635  fEMCALTimeRecalibrationFactors->SetOwner(true);
3636  }
3637  if(!fEMCALTimeRecalibrationFactors->IsOwner()) {
3638  // Must claim ownership since the new objects are owend by this instance
3639  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
3640  }
3641  for(int i = 0; i < map->GetEntries(); i++){
3642  TH1F *hist = dynamic_cast<TH1F *>(map->At(i));
3643  if(!hist) continue;
3644  this->SetEMCALChannelTimeRecalibrationFactors(i, hist);
3645  }
3646 }
3647 
3648 void AliEMCALRecoUtils::SetEMCALChannelTimeRecalibrationFactors(Int_t bc, const TH1F* h){
3649  if(!fEMCALTimeRecalibrationFactors){
3650  fEMCALTimeRecalibrationFactors = new TObjArray(bc);
3651  fEMCALTimeRecalibrationFactors->SetOwner(true);
3652  }
3653  if(fEMCALTimeRecalibrationFactors->GetEntries() <= bc) fEMCALTimeRecalibrationFactors->Expand(bc+1);
3654  if(fEMCALTimeRecalibrationFactors->At(bc)) fEMCALTimeRecalibrationFactors->RemoveAt(bc);
3655  TH1F *clone = new TH1F(*h);
3656  clone->SetDirectory(NULL);
3657  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
3658 }
3659 
3660 void AliEMCALRecoUtils::SetEMCALL1PhaseInTimeRecalibrationForAllSM(const TObjArray *map) {
3661  if(fEMCALL1PhaseInTimeRecalibration) fEMCALL1PhaseInTimeRecalibration->Clear();
3662  else {
3663  fEMCALL1PhaseInTimeRecalibration = new TObjArray(map->GetEntries());
3664  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3665  }
3666  if(!fEMCALL1PhaseInTimeRecalibration->IsOwner()) {
3667  // Must claim ownership since the new objects are owend by this instance
3668  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3669  }
3670  for(int i = 0; i < map->GetEntries(); i++) {
3671  TH1C *hist = dynamic_cast<TH1C *>(map->At(i));
3672  if(!hist) continue;
3673  SetEMCALL1PhaseInTimeRecalibrationForAllSM(hist);
3674  }
3675 }
3676 
3677 void AliEMCALRecoUtils::SetEMCALL1PhaseInTimeRecalibrationForAllSM(const TH1C* h) {
3678  if(!fEMCALL1PhaseInTimeRecalibration){
3679  fEMCALL1PhaseInTimeRecalibration = new TObjArray(1);
3680  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3681  }
3682  if(fEMCALL1PhaseInTimeRecalibration->GetEntries()<1) fEMCALL1PhaseInTimeRecalibration->Expand(1);
3683  TH1C *clone = new TH1C(*h);
3684  clone->SetDirectory(NULL);
3685  fEMCALL1PhaseInTimeRecalibration->AddAt(clone,0);
3686 }
3687 
3691 //___________________________________________________
3692 void AliEMCALRecoUtils::Print(const Option_t *) const
3693 {
3694  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3695  printf("AliEMCALRecoUtils Settings: \n");
3696  printf("\tMisalignment shifts\n");
3697  for (Int_t i=0; i<5; i++) printf("\t\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
3698  fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
3699  fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
3700  printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
3701  if (fNonLinearityFunction != 3) // print only if not kNoCorrection
3702  for (Int_t i=0; i<10; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
3703 
3704  printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
3705 
3706  printf("\tMatching criteria: ");
3707  if (fCutEtaPhiSum) {
3708  printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
3709  } else if (fCutEtaPhiSeparate) {
3710  printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
3711  } else {
3712  printf("\t\tError\n");
3713  printf("\t\tplease specify your cut criteria\n");
3714  printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
3715  printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
3716  }
3717 
3718  printf("\tMass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
3719  printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
3720 
3721  printf("\tTrack cuts: \n");
3722  printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
3723  printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
3724  printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
3725  printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
3726  printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
3727  printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
3728  printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
3729  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3730 }
AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
Definition: AliEMCALRecoUtils.cxx:1805
AliEMCALRecoUtils::IsL1PhaseInTimeRecalibrationOn
Bool_t IsL1PhaseInTimeRecalibrationOn() const
Definition: AliEMCALRecoUtils.h:217
AliEMCALRecoUtils::fMass
Double_t fMass
Mass hypothesis of the track.
Definition: AliEMCALRecoUtils.h:524
AliEMCALRecoUtils::fMCGenerToAcceptForTrack
Bool_t fMCGenerToAcceptForTrack
Activate the removal of tracks entering the track matching that come from a particular generator...
Definition: AliEMCALRecoUtils.h:549
Double_t
double Double_t
Definition: External.C:58
AliEMCALRecoUtils::CheckCellFiducialRegion
Bool_t CheckCellFiducialRegion(const AliEMCALGeometry *geom, const AliVCluster *cluster, AliVCaloCells *cells)
Definition: AliEMCALRecoUtils.cxx:462
AliEMCALRecoUtils::fAODTPCOnlyTracks
Bool_t fAODTPCOnlyTracks
Match with TPC only tracks.
Definition: AliEMCALRecoUtils.h:512
AliEMCALRecoUtils.h
AliEMCALRecoUtils::fCutRequireTPCRefit
Bool_t fCutRequireTPCRefit
Require TPC refit.
Definition: AliEMCALRecoUtils.h:537
TH2F
Definition: External.C:236
AliEMCALRecoUtils::SetRequireITSStandAlone
void SetRequireITSStandAlone(Bool_t b=kFALSE)
Definition: AliEMCALRecoUtils.h:379
AliEMCALRecoUtils::IsClusterMatched
Bool_t IsClusterMatched(Int_t clsIndex) const
Definition: AliEMCALRecoUtils.cxx:3097
AliEMCALRecoUtils::GetEMCALChannelTimeRecalibrationFactor
Float_t GetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Bool_t isLGon=kFALSE) const
Definition: AliEMCALRecoUtils.h:199
AliEMCALRecoUtils::fMatchedClusterIndex
TArrayI * fMatchedClusterIndex
Array that stores indexes of matched clusters.
Definition: AliEMCALRecoUtils.h:515
AliEMCALRecoUtils::fEMCALTimeRecalibrationFactors
TObjArray * fEMCALTimeRecalibrationFactors
Array of histograms with map of time recalibration factors, EMCAL.
Definition: AliEMCALRecoUtils.h:475
AliEMCALRecoUtils::SetMaxDCAToVertexXY
void SetMaxDCAToVertexXY(Float_t dist=1e10)
Definition: AliEMCALRecoUtils.h:376
AliEMCALRecoUtils::fRandom
TRandom3 fRandom
Random generator for cluster energy smearing.
Definition: AliEMCALRecoUtils.h:465
AliEMCALRecoUtils::fTimeRecalibration
Bool_t fTimeRecalibration
Switch on or off the time recalibration.
Definition: AliEMCALRecoUtils.h:474
AliEMCALRecoUtils::fCutEtaPhiSum
Bool_t fCutEtaPhiSum
Place cut on sqrt(dEta^2+dPhi^2)
Definition: AliEMCALRecoUtils.h:518
AliEMCALRecoUtils::FindMatchedClusterInClusterArr
Int_t FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, const TObjArray *clusterArr, Float_t &dEta, Float_t &dPhi)
Definition: AliEMCALRecoUtils.cxx:2925
AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibration
void InitEMCALL1PhaseInTimeRecalibration()
Definition: AliEMCALRecoUtils.cxx:1583
AliEMCALRecoUtils::fCutMinNClusterTPC
Int_t fCutMinNClusterTPC
Min number of tpc clusters.
Definition: AliEMCALRecoUtils.h:533
AliEMCALRecoUtils::fCutMaxDCAToVertexZ
Float_t fCutMaxDCAToVertexZ
Track-to-vertex cut in max absolute distance in z-plane.
Definition: AliEMCALRecoUtils.h:541
AliEMCALRecoUtils::SetRequireITSRefit
void SetRequireITSRefit(Bool_t b=kFALSE)
Definition: AliEMCALRecoUtils.h:374
energy
energy
Definition: HFPtSpectrum.C:44
AliEMCALRecoUtils::SetEMCALL1PhaseInTimeRecalibrationForSM
void SetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM, Int_t c=0)
Definition: AliEMCALRecoUtils.h:229
AliEMCALRecoUtils::SetDCAToVertex2D
void SetDCAToVertex2D(Bool_t b=kFALSE)
Definition: AliEMCALRecoUtils.h:378
AliEMCALRecoUtils::SetEMCALChannelRecalibrationFactor
void SetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
Definition: AliEMCALRecoUtils.h:178
AliEMCALRecoUtils::fNonLinearThreshold
Int_t fNonLinearThreshold
Non linearity threshold value for kBeamTest non linearity function.
Definition: AliEMCALRecoUtils.h:460
AliEMCALRecoUtils::fITSTrackSA
Bool_t fITSTrackSA
If track matching is to be done with ITS tracks standing alone.
Definition: AliEMCALRecoUtils.h:527
AliEMCALRecoUtils::SetITSTrackSA
void SetITSTrackSA(Bool_t isITS)
Definition: AliEMCALRecoUtils.h:350
AliEMCALRecoUtils::IsExoticCluster
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
Definition: AliEMCALRecoUtils.cxx:703
AliEMCALRecoUtils::fMisalRotShift
Float_t fMisalRotShift[15]
Cluster position rotation shift parameters.
Definition: AliEMCALRecoUtils.h:447
AliEMCALRecoUtils::SetMinNClustersTPC
void SetMinNClustersTPC(Int_t min=-1)
Definition: AliEMCALRecoUtils.h:369
AliEMCALRecoUtils::fRecalDistToBadChannels
Bool_t fRecalDistToBadChannels
Calculate distance from highest energy tower of cluster to closes bad channel.
Definition: AliEMCALRecoUtils.h:487
AliEMCALRecoUtils::fRecalibration
Bool_t fRecalibration
Switch on or off the recalibration.
Definition: AliEMCALRecoUtils.h:469
AliEMCALRecoUtils::fMatchedTrackIndex
TArrayI * fMatchedTrackIndex
Array that stores indexes of matched tracks.
Definition: AliEMCALRecoUtils.h:514
AliEMCALRecoUtils::SetMaxDCAToVertexZ
void SetMaxDCAToVertexZ(Float_t dist=1e10)
Definition: AliEMCALRecoUtils.h:377
AliEMCALRecoUtils::IsBadChannelsRemovalSwitchedOn
Bool_t IsBadChannelsRemovalSwitchedOn() const
Definition: AliEMCALRecoUtils.h:253
AliEMCALRecoUtils::IsExoticCell
Bool_t IsExoticCell(Int_t absId, AliVCaloCells *cells, Int_t bc=-1)
Definition: AliEMCALRecoUtils.cxx:669
AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal
void RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Definition: AliEMCALRecoUtils.cxx:1931
AliEMCALRecoUtils::fExoticCellDiffTime
Float_t fExoticCellDiffTime
If time of candidate to exotic and close cell is too different (in ns), it must be noisy...
Definition: AliEMCALRecoUtils.h:503
AliEMCALRecoUtils::fResidualPhi
TArrayF * fResidualPhi
Array that stores the residual phi.
Definition: AliEMCALRecoUtils.h:517
AliEMCALRecoUtils::RecalibrateCells
void RecalibrateCells(AliVCaloCells *cells, Int_t bc)
Definition: AliEMCALRecoUtils.cxx:1700
AliEMCALRecoUtils::fCutMinTrackPt
Double_t fCutMinTrackPt
Cut on track pT.
Definition: AliEMCALRecoUtils.h:532
AliEMCALRecoUtils::IsRecalibrationOn
Bool_t IsRecalibrationOn() const
Definition: AliEMCALRecoUtils.h:165
AliEMCALRecoUtils
Some utilities for cluster and cell treatment.
Definition: AliEMCALRecoUtils.h:60
AliEMCALRecoUtils::fCutMaxChi2PerClusterITS
Float_t fCutMaxChi2PerClusterITS
Max its fit chi2 per its cluster.
Definition: AliEMCALRecoUtils.h:536
AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts
void RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster, Float_t cellEcut, Float_t cellTimeCut, Int_t bc, Float_t &enAfterCuts)
Definition: AliEMCALRecoUtils.cxx:2578
AliEMCALRecoUtils::fRemoveBadChannels
Bool_t fRemoveBadChannels
Check the channel status provided and remove clusters with bad channels.
Definition: AliEMCALRecoUtils.h:486
AliEMCALRecoUtils::IsAccepted
Bool_t IsAccepted(AliESDtrack *track)
Definition: AliEMCALRecoUtils.cxx:3221
AliEMCALRecoUtils::fMisalTransShift
Float_t fMisalTransShift[15]
Cluster position translation shift parameters.
Definition: AliEMCALRecoUtils.h:446
AliEMCALRecoUtils::SetAcceptKinkDaughters
void SetAcceptKinkDaughters(Bool_t b=kTRUE)
Definition: AliEMCALRecoUtils.h:375
AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
Definition: AliEMCALRecoUtils.cxx:2101
AliEMCALRecoUtils::fClusterWindow
Double_t fClusterWindow
Select clusters in the window to be matched.
Definition: AliEMCALRecoUtils.h:523
AliEMCALRecoUtils::fCutRequireITSRefit
Bool_t fCutRequireITSRefit
Require ITS refit.
Definition: AliEMCALRecoUtils.h:538
AliEMCALRecoUtils::fNonLinearityParams
Float_t fNonLinearityParams[10]
Parameters for the non linearity function.
Definition: AliEMCALRecoUtils.h:459
AliEMCALRecoUtils::GetCellWeight
Float_t GetCellWeight(Float_t eCell, Float_t eCluster) const
Definition: AliEMCALRecoUtils.cxx:1384
AliEMCALRecoUtils::fUseL1PhaseInTimeRecalibration
Bool_t fUseL1PhaseInTimeRecalibration
Switch on or off the L1 phase in time recalibration.
Definition: AliEMCALRecoUtils.h:479
AliEMCALRecoUtils::GetEMCALChannelRecalibrationFactor
Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow) const
Definition: AliEMCALRecoUtils.h:174
AliEMCALRecoUtils::SetClusterMatchedToTrack
void SetClusterMatchedToTrack(const AliVEvent *event)
Definition: AliEMCALRecoUtils.cxx:3444
AliEMCALRecoUtils::fCutDCAToVertex2D
Bool_t fCutDCAToVertex2D
If true a 2D DCA cut is made.
Definition: AliEMCALRecoUtils.h:542
AliEMCALRecoUtils::GetMatchedTrackIndex
Int_t GetMatchedTrackIndex(Int_t clsIndex)
Definition: AliEMCALRecoUtils.cxx:3066
AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex
void RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Definition: AliEMCALRecoUtils.cxx:2018
Int_t
int Int_t
Definition: External.C:63
AliEMCALRecoUtils::Print
void Print(const Option_t *) const
Definition: AliEMCALRecoUtils.cxx:3692
AliEMCALRecoUtils::FindMatchedPosForCluster
UInt_t FindMatchedPosForCluster(Int_t clsIndex) const
Definition: AliEMCALRecoUtils.cxx:3129
AliEMCALRecoUtils::AcceptCalibrateCell
Bool_t AcceptCalibrateCell(Int_t absId, Int_t bc, Float_t &amp, Double_t &time, AliVCaloCells *cells)
Definition: AliEMCALRecoUtils.cxx:397
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliEMCALRecoUtils::GetEMCALChannelStatus
Bool_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow, Int_t &status) const
Definition: AliEMCALRecoUtils.cxx:1254
AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters
void RecalculateClusterShowerShapeParameters(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
Definition: AliEMCALRecoUtils.cxx:2548
Float_t
float Float_t
Definition: External.C:68
AliEMCALRecoUtils::fCutR
Float_t fCutR
sqrt(dEta^2+dPhi^2) cut on matching
Definition: AliEMCALRecoUtils.h:520
AliEMCALRecoUtils::fNCellsFromEMCALBorder
Int_t fNCellsFromEMCALBorder
Number of cells from EMCAL border the cell with maximum amplitude has to be.
Definition: AliEMCALRecoUtils.h:496
AliEMCALRecoUtils::fCellsRecalibrated
Bool_t fCellsRecalibrated
Internal bool to check if cells (time/energy) where recalibrated and not recalibrate them when recalc...
Definition: AliEMCALRecoUtils.h:468
AliEMCALRecoUtils::fRejectExoticCluster
Bool_t fRejectExoticCluster
Switch on or off exotic cluster rejection.
Definition: AliEMCALRecoUtils.h:500
AliEMCALRecoUtils::IsTimeRecalibrationOn
Bool_t IsTimeRecalibrationOn() const
Definition: AliEMCALRecoUtils.h:192
AliEMCALRecoUtils::fParticleType
Int_t fParticleType
Particle type for depth calculation, see enum ParticleType.
Definition: AliEMCALRecoUtils.h:448
AliEMCALRecoUtils::SetEMCALBadChannelStatusSelection
void SetEMCALBadChannelStatusSelection(Bool_t all, Bool_t dead, Bool_t hot, Bool_t warm)
Definition: AliEMCALRecoUtils.cxx:1235
TH2D
Definition: External.C:228
AliEMCALRecoUtils::SetMaxChi2PerClusterTPC
void SetMaxChi2PerClusterTPC(Float_t max=1e10)
Definition: AliEMCALRecoUtils.h:371
AliEMCALRecoUtils::fConstantTimeShift
Float_t fConstantTimeShift
Apply a 600 ns (+15.8) time shift in case of simulation, shift in ns.
Definition: AliEMCALRecoUtils.h:473
AliEMCALRecoUtils::GetMatchedClusterIndex
Int_t GetMatchedClusterIndex(Int_t trkIndex)
Definition: AliEMCALRecoUtils.cxx:3081
AliEMCALRecoUtils::fNonLinearityFunction
Int_t fNonLinearityFunction
Non linearity function choice, see enum NonlinearityFunctions.
Definition: AliEMCALRecoUtils.h:458
AliEMCALRecoUtils::fCutRequireITSStandAlone
Bool_t fCutRequireITSStandAlone
Require ITSStandAlone.
Definition: AliEMCALRecoUtils.h:543
AliEMCALRecoUtils::fEMCALL1PhaseInTimeRecalibration
TObjArray * fEMCALL1PhaseInTimeRecalibration
Histogram with map of L1 phase per SM, EMCAL.
Definition: AliEMCALRecoUtils.h:480
AliEMCALRecoUtils::fEMCALRecalibrationFactors
TObjArray * fEMCALRecalibrationFactors
Array of histograms with map of recalibration factors, EMCAL.
Definition: AliEMCALRecoUtils.h:470
AliEMCALRecoUtils::fUseRunCorrectionFactors
Bool_t fUseRunCorrectionFactors
Use Run Dependent Correction.
Definition: AliEMCALRecoUtils.h:483
AliEMCALRecoUtils::SetTracksMatchedToCluster
void SetTracksMatchedToCluster(const AliVEvent *event)
Definition: AliEMCALRecoUtils.cxx:3487
AliEMCALRecoUtils::fCutMaxChi2PerClusterTPC
Float_t fCutMaxChi2PerClusterTPC
Max tpc fit chi2 per tpc cluster.
Definition: AliEMCALRecoUtils.h:535
AliEMCALRecoUtils::fExoticCellMinAmplitude
Float_t fExoticCellMinAmplitude
Check for exotic only if amplitud is larger than this value.
Definition: AliEMCALRecoUtils.h:504
AliEMCALRecoUtils::RecalibrateClusterEnergy
void RecalibrateClusterEnergy(const AliEMCALGeometry *geom, AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=-1)
Definition: AliEMCALRecoUtils.cxx:1615
AliEMCALRecoUtils::fW0
Float_t fW0
Energy weight used in cluster position and shower shape calculations.
Definition: AliEMCALRecoUtils.h:451
AliEMCALRecoUtils::fAODHybridTracks
Bool_t fAODHybridTracks
Match with hybrid.
Definition: AliEMCALRecoUtils.h:511
AliEMCALRecoUtils::GetMatchedResiduals
void GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
Definition: AliEMCALRecoUtils.cxx:3020
AliEMCALRecoUtils::operator=
AliEMCALRecoUtils & operator=(const AliEMCALRecoUtils &)
Definition: AliEMCALRecoUtils.cxx:175
AliEMCALRecoUtils::SetRequireTPCRefit
void SetRequireTPCRefit(Bool_t b=kFALSE)
Definition: AliEMCALRecoUtils.h:373
AliEMCALRecoUtils::fResidualEta
TArrayF * fResidualEta
Array that stores the residual eta.
Definition: AliEMCALRecoUtils.h:516
Short_t
short Short_t
Definition: External.C:23
AliEMCALRecoUtils::fMCGenerToAccept
TString fMCGenerToAccept[5]
List with name of generators that should not be included.
Definition: AliEMCALRecoUtils.h:548
AliEMCALRecoUtils::fSmearClusterEnergy
Bool_t fSmearClusterEnergy
Smear cluster energy, to be done only for simulated data to match real data.
Definition: AliEMCALRecoUtils.h:463
AliEMCALRecoUtils::IsGoodCluster
Bool_t IsGoodCluster(AliVCluster *cluster, const AliEMCALGeometry *geom, AliVCaloCells *cells, Int_t bc=-1)
Definition: AliEMCALRecoUtils.cxx:3196
AliEMCALRecoUtils::GetMaxEnergyCell
void GetMaxEnergyCell(const AliEMCALGeometry *geom, AliVCaloCells *cells, const AliVCluster *clu, Int_t &absId, Int_t &iSupMod, Int_t &ieta, Int_t &iphi, Bool_t &shared)
Definition: AliEMCALRecoUtils.cxx:1308
AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors
void InitEMCALTimeRecalibrationFactors()
Definition: AliEMCALRecoUtils.cxx:1513
AliEMCALRecoUtils::fStepSurface
Double_t fStepSurface
Length of step to extrapolate tracks to EMCal surface.
Definition: AliEMCALRecoUtils.h:525
AliEMCALRecoUtils::fAODFilterMask
UInt_t fAODFilterMask
Filter mask to select AOD tracks. Refer to $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C.
Definition: AliEMCALRecoUtils.h:510
AliEMCALRecoUtils::ExtrapolateTrackToCluster
static Bool_t ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, const AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
Definition: AliEMCALRecoUtils.h:312
AliEMCALRecoUtils::GetMatchedClusterResiduals
void GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
Definition: AliEMCALRecoUtils.cxx:3045
AliEMCALRecoUtils::fSmearClusterParam
Float_t fSmearClusterParam[3]
Energy smearing parameters.
Definition: AliEMCALRecoUtils.h:464
AliEMCALRecoUtils::fCutAcceptKinkDaughters
Bool_t fCutAcceptKinkDaughters
Accepting kink daughters?
Definition: AliEMCALRecoUtils.h:539
AliEMCALRecoUtils::RecalculateClusterPosition
void RecalculateClusterPosition(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Definition: AliEMCALRecoUtils.cxx:1907
AliEMCALRecoUtils::fExoticCellFraction
Float_t fExoticCellFraction
Good cell if fraction < 1-ecross/ecell.
Definition: AliEMCALRecoUtils.h:502
AliEMCALRecoUtils::RecalculateClusterPID
void RecalculateClusterPID(AliVCluster *cluster)
Definition: AliEMCALRecoUtils.cxx:2189
AliEMCALRecoUtils::fPosAlgo
Int_t fPosAlgo
Position recalculation algorithm, see enum PositionAlgorithms.
Definition: AliEMCALRecoUtils.h:449
AliEMCALRecoUtils::SetEMCALL1PhaseInTimeRecalibrationForAllSM
void SetEMCALL1PhaseInTimeRecalibrationForAllSM(const TObjArray *map)
Definition: AliEMCALRecoUtils.cxx:3660
AliEMCALRecoUtils::fPIDUtils
AliEMCALPIDUtils * fPIDUtils
Recalculate PID parameters.
Definition: AliEMCALRecoUtils.h:507
AliEMCALRecoUtils::fRejectExoticCells
Bool_t fRejectExoticCells
Remove exotic cells.
Definition: AliEMCALRecoUtils.h:501
AliEMCALRecoUtils::FindMatchedClusterInEvent
Int_t FindMatchedClusterInEvent(const AliESDtrack *track, const AliVEvent *event, const AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi)
Definition: AliEMCALRecoUtils.cxx:2851
AliEMCALRecoUtils::IsTrackMatched
Bool_t IsTrackMatched(Int_t trkIndex) const
Definition: AliEMCALRecoUtils.cxx:3113
AliEMCALRecoUtils::RecalibrateCellTimeL1Phase
void RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t &time) const
Definition: AliEMCALRecoUtils.cxx:1766
AliEMCALRecoUtils::fCutMaxDCAToVertexXY
Float_t fCutMaxDCAToVertexXY
Track-to-vertex cut in max absolute distance in xy-plane.
Definition: AliEMCALRecoUtils.h:540
AliEMCALRecoUtils::FindMatches
void FindMatches(AliVEvent *event, TObjArray *clusterArr=0x0, const AliEMCALGeometry *geom=0x0, AliMCEvent *mc=0x0)
Definition: AliEMCALRecoUtils.cxx:2612
UShort_t
unsigned short UShort_t
Definition: External.C:28
AliEMCALRecoUtils::fEMCALBadChannelMap
TObjArray * fEMCALBadChannelMap
Array of histograms with map of bad channels, EMCAL.
Definition: AliEMCALRecoUtils.h:488
AliEMCALRecoUtils::SetEMCALChannelTimeRecalibrationFactors
void SetEMCALChannelTimeRecalibrationFactors(const TObjArray *map)
Definition: AliEMCALRecoUtils.cxx:3631
AliEMCALRecoUtils::fStepCluster
Double_t fStepCluster
Length of step to extrapolate tracks to clusters.
Definition: AliEMCALRecoUtils.h:526
AliEMCALRecoUtils::fNoEMCALBorderAtEta0
Bool_t fNoEMCALBorderAtEta0
Do fiducial cut in EMCAL region eta = 0?
Definition: AliEMCALRecoUtils.h:497
Option_t
const char Option_t
Definition: External.C:48
AliEMCALRecoUtils::fNMCGenerToAccept
Int_t fNMCGenerToAccept
Number of MC generators that should not be included in analysis.
Definition: AliEMCALRecoUtils.h:547
AliEMCALRecoUtils::fBadStatusSelection
Bool_t fBadStatusSelection[4]
Definition: AliEMCALRecoUtils.h:489
AliEMCALRecoUtils::SetEMCALChannelTimeRecalibrationFactor
void SetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Double_t c=0, Bool_t isLGon=kFALSE)
Definition: AliEMCALRecoUtils.h:203
AliEMCALRecoUtils::fCutRequireITSpureSA
Bool_t fCutRequireITSpureSA
ITS pure standalone tracks.
Definition: AliEMCALRecoUtils.h:544
AliEMCALRecoUtils::GetEMCALL1PhaseInTimeRecalibrationForSM
Int_t GetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM) const
Definition: AliEMCALRecoUtils.h:225
Bool_t
bool Bool_t
Definition: External.C:53
AliEMCALRecoUtils::fCutMinNClusterITS
Int_t fCutMinNClusterITS
Min number of its clusters.
Definition: AliEMCALRecoUtils.h:534
AliEMCALRecoUtils::CorrectClusterEnergyLinearity
Float_t CorrectClusterEnergyLinearity(AliVCluster *clu)
Definition: AliEMCALRecoUtils.cxx:771
AliEMCALRecoUtils::FindMatchedPosForTrack
UInt_t FindMatchedPosForTrack(Int_t trkIndex) const
Definition: AliEMCALRecoUtils.cxx:3159
AliEMCALRecoUtils::ClusterContainsBadChannel
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)
Definition: AliEMCALRecoUtils.cxx:551
AliEMCALRecoUtils::fEMCalSurfaceDistance
Double_t fEMCalSurfaceDistance
EMCal surface distance (= 430 by default, the last 10 cm are propagated on a cluster-track pair basis...
Definition: AliEMCALRecoUtils.h:528
AliEMCALRecoUtils::SetEMCALChannelRecalibrationFactors
void SetEMCALChannelRecalibrationFactors(const TObjArray *map)
Definition: AliEMCALRecoUtils.cxx:3573
AliEMCALRecoUtils::GetECross
Float_t GetECross(Int_t absID, Double_t tcell, AliVCaloCells *cells, Int_t bc)
Definition: AliEMCALRecoUtils.cxx:594
AliEMCALRecoUtils::fLowGain
Bool_t fLowGain
Switch on or off calibration with low gain channels.
Definition: AliEMCALRecoUtils.h:476
AliEMCALRecoUtils::SetEMCALChannelStatusMap
void SetEMCALChannelStatusMap(const TObjArray *map)
Definition: AliEMCALRecoUtils.cxx:3602
AliEMCALRecoUtils::fCutPhi
Float_t fCutPhi
dPhi cut on matching
Definition: AliEMCALRecoUtils.h:522
AliEMCALRecoUtils::fTrackCutsType
Int_t fTrackCutsType
ESD track cuts type for matching, see enum TrackCutsType.
Definition: AliEMCALRecoUtils.h:531
AliEMCALRecoUtils::RecalibrateCellTime
void RecalibrateCellTime(Int_t absId, Int_t bc, Double_t &time, Bool_t isLGon=kFALSE) const
Definition: AliEMCALRecoUtils.cxx:1748
AliEMCALRecoUtils::fCutEtaPhiSeparate
Bool_t fCutEtaPhiSeparate
Cut on dEta and dPhi separately.
Definition: AliEMCALRecoUtils.h:519
AliEMCALRecoUtils::fCutEta
Float_t fCutEta
dEta cut on matching
Definition: AliEMCALRecoUtils.h:521
AliEMCALRecoUtils::SmearClusterEnergy
Float_t SmearClusterEnergy(const AliVCluster *clu)
Definition: AliEMCALRecoUtils.cxx:740