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