Line data Source code
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 <TGeoManager.h>
18 : #include <TGeoMatrix.h>
19 : #include <TGeoBBox.h>
20 : #include <TH2F.h>
21 : #include <TArrayI.h>
22 : #include <TArrayF.h>
23 : #include <TObjArray.h>
24 :
25 : // STEER includes
26 : #include "AliVCluster.h"
27 : #include "AliVCaloCells.h"
28 : #include "AliLog.h"
29 : #include "AliPID.h"
30 : #include "AliESDEvent.h"
31 : #include "AliAODEvent.h"
32 : #include "AliESDtrack.h"
33 : #include "AliAODTrack.h"
34 : #include "AliExternalTrackParam.h"
35 : #include "AliESDfriendTrack.h"
36 : #include "AliTrackerBase.h"
37 :
38 : // EMCAL includes
39 : #include "AliEMCALRecoUtils.h"
40 : #include "AliEMCALGeometry.h"
41 : #include "AliTrackerBase.h"
42 : #include "AliEMCALPIDUtils.h"
43 :
44 : /// \cond CLASSIMP
45 72 : ClassImp(AliEMCALRecoUtils) ;
46 : /// \endcond
47 :
48 : ///
49 : /// Constructor.
50 : /// Initialize all constant values which have to be used
51 : /// during Reco algorithm execution
52 : ///
53 : //_____________________________________
54 0 : AliEMCALRecoUtils::AliEMCALRecoUtils():
55 0 : fParticleType(0), fPosAlgo(0), fW0(0),
56 0 : fNonLinearityFunction(0), fNonLinearThreshold(0),
57 0 : fSmearClusterEnergy(kFALSE), fRandom(),
58 0 : fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
59 0 : fConstantTimeShift(0), fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fLowGain(kFALSE),
60 0 : fUseL1PhaseInTimeRecalibration(kFALSE), fEMCALL1PhaseInTimeRecalibration(),
61 0 : fUseRunCorrectionFactors(kFALSE),
62 0 : fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
63 0 : fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
64 0 : fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
65 0 : fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
66 0 : fPIDUtils(), fAODFilterMask(0),
67 0 : fAODHybridTracks(0), fAODTPCOnlyTracks(0),
68 0 : fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
69 0 : fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
70 0 : fCutR(0), fCutEta(0), fCutPhi(0),
71 0 : fClusterWindow(0), fMass(0),
72 0 : fStepSurface(0), fStepCluster(0),
73 0 : fITSTrackSA(kFALSE), fEMCalSurfaceDistance(440.),
74 0 : fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
75 0 : fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
76 0 : fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
77 0 : fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
78 0 : fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE)
79 0 : {
80 : // Init parameters
81 0 : InitParameters();
82 :
83 : // Track matching arrays init
84 0 : fMatchedTrackIndex = new TArrayI();
85 0 : fMatchedClusterIndex = new TArrayI();
86 0 : fResidualPhi = new TArrayF();
87 0 : fResidualEta = new TArrayF();
88 0 : fPIDUtils = new AliEMCALPIDUtils();
89 0 : }
90 :
91 : //
92 : // Copy constructor.
93 : //
94 : //______________________________________________________________________
95 : AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
96 0 : : TNamed(reco),
97 0 : fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
98 0 : fNonLinearityFunction(reco.fNonLinearityFunction), fNonLinearThreshold(reco.fNonLinearThreshold),
99 0 : fSmearClusterEnergy(reco.fSmearClusterEnergy), fRandom(),
100 0 : fCellsRecalibrated(reco.fCellsRecalibrated),
101 0 : fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
102 0 : fConstantTimeShift(reco.fConstantTimeShift),
103 0 : fTimeRecalibration(reco.fTimeRecalibration), fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
104 0 : fLowGain(reco.fLowGain),
105 0 : fUseL1PhaseInTimeRecalibration(reco.fUseL1PhaseInTimeRecalibration),
106 0 : fEMCALL1PhaseInTimeRecalibration(reco.fEMCALL1PhaseInTimeRecalibration),
107 0 : fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),
108 0 : fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
109 0 : fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
110 0 : fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
111 0 : fRejectExoticCluster(reco.fRejectExoticCluster), fRejectExoticCells(reco.fRejectExoticCells),
112 0 : fExoticCellFraction(reco.fExoticCellFraction), fExoticCellDiffTime(reco.fExoticCellDiffTime),
113 0 : fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
114 0 : fPIDUtils(reco.fPIDUtils), fAODFilterMask(reco.fAODFilterMask),
115 0 : fAODHybridTracks(reco.fAODHybridTracks), fAODTPCOnlyTracks(reco.fAODTPCOnlyTracks),
116 0 : fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
117 0 : fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
118 0 : fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
119 0 : fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
120 0 : fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
121 0 : fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
122 0 : fClusterWindow(reco.fClusterWindow),
123 0 : fMass(reco.fMass), fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
124 0 : fITSTrackSA(reco.fITSTrackSA), fEMCalSurfaceDistance(440.),
125 0 : fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
126 0 : fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
127 0 : fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
128 0 : fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
129 0 : fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
130 0 : fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
131 0 : fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone), fCutRequireITSpureSA(reco.fCutRequireITSpureSA)
132 0 : {
133 0 : for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
134 0 : fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
135 0 : for (Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
136 0 : for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
137 0 : }
138 :
139 : ///
140 : /// Assignment operator.
141 : ///
142 : //______________________________________________________________________
143 : AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
144 : {
145 0 : if (this == &reco)return *this;
146 0 : ((TNamed *)this)->operator=(reco);
147 :
148 0 : for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
149 0 : fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
150 0 : for (Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
151 0 : for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
152 :
153 0 : fParticleType = reco.fParticleType;
154 0 : fPosAlgo = reco.fPosAlgo;
155 0 : fW0 = reco.fW0;
156 :
157 0 : fNonLinearityFunction = reco.fNonLinearityFunction;
158 0 : fNonLinearThreshold = reco.fNonLinearThreshold;
159 0 : fSmearClusterEnergy = reco.fSmearClusterEnergy;
160 :
161 0 : fCellsRecalibrated = reco.fCellsRecalibrated;
162 0 : fRecalibration = reco.fRecalibration;
163 0 : fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
164 :
165 0 : fConstantTimeShift = reco.fConstantTimeShift;
166 0 : fTimeRecalibration = reco.fTimeRecalibration;
167 0 : fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
168 0 : fLowGain = reco.fLowGain;
169 :
170 0 : fUseL1PhaseInTimeRecalibration = reco.fUseL1PhaseInTimeRecalibration;
171 0 : fEMCALL1PhaseInTimeRecalibration = reco.fEMCALL1PhaseInTimeRecalibration;
172 :
173 0 : fUseRunCorrectionFactors = reco.fUseRunCorrectionFactors;
174 :
175 0 : fRemoveBadChannels = reco.fRemoveBadChannels;
176 0 : fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
177 0 : fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
178 :
179 0 : fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
180 0 : fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
181 :
182 0 : fRejectExoticCluster = reco.fRejectExoticCluster;
183 0 : fRejectExoticCells = reco.fRejectExoticCells;
184 0 : fExoticCellFraction = reco.fExoticCellFraction;
185 0 : fExoticCellDiffTime = reco.fExoticCellDiffTime;
186 0 : fExoticCellMinAmplitude = reco.fExoticCellMinAmplitude;
187 :
188 0 : fPIDUtils = reco.fPIDUtils;
189 :
190 0 : fAODFilterMask = reco.fAODFilterMask;
191 0 : fAODHybridTracks = reco.fAODHybridTracks;
192 0 : fAODTPCOnlyTracks = reco.fAODTPCOnlyTracks;
193 :
194 0 : fCutEtaPhiSum = reco.fCutEtaPhiSum;
195 0 : fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
196 0 : fCutR = reco.fCutR;
197 0 : fCutEta = reco.fCutEta;
198 0 : fCutPhi = reco.fCutPhi;
199 0 : fClusterWindow = reco.fClusterWindow;
200 0 : fMass = reco.fMass;
201 0 : fStepSurface = reco.fStepSurface;
202 0 : fStepCluster = reco.fStepCluster;
203 0 : fITSTrackSA = reco.fITSTrackSA;
204 0 : fEMCalSurfaceDistance = reco.fEMCalSurfaceDistance;
205 :
206 0 : fTrackCutsType = reco.fTrackCutsType;
207 0 : fCutMinTrackPt = reco.fCutMinTrackPt;
208 0 : fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
209 0 : fCutMinNClusterITS = reco.fCutMinNClusterITS;
210 0 : fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
211 0 : fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
212 0 : fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
213 0 : fCutRequireITSRefit = reco.fCutRequireITSRefit;
214 0 : fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
215 0 : fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
216 0 : fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
217 0 : fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
218 0 : fCutRequireITSStandAlone = reco.fCutRequireITSStandAlone;
219 0 : fCutRequireITSpureSA = reco.fCutRequireITSpureSA;
220 :
221 : //
222 : // Assign or copy construct the different TArrays
223 : //
224 0 : if (reco.fResidualEta)
225 : {
226 0 : if (fResidualEta)
227 0 : *fResidualEta = *reco.fResidualEta;
228 : else
229 0 : fResidualEta = new TArrayF(*reco.fResidualEta);
230 : }
231 : else
232 : {
233 0 : if (fResidualEta) delete fResidualEta;
234 0 : fResidualEta = 0;
235 : }
236 :
237 0 : if (reco.fResidualPhi)
238 : {
239 0 : if (fResidualPhi)
240 0 : *fResidualPhi = *reco.fResidualPhi;
241 : else
242 0 : fResidualPhi = new TArrayF(*reco.fResidualPhi);
243 : }
244 : else
245 : {
246 0 : if (fResidualPhi) delete fResidualPhi;
247 0 : fResidualPhi = 0;
248 : }
249 :
250 0 : if (reco.fMatchedTrackIndex)
251 : {
252 0 : if (fMatchedTrackIndex)
253 0 : *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
254 : else
255 0 : fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
256 : }
257 : else
258 : {
259 0 : if (fMatchedTrackIndex) delete fMatchedTrackIndex;
260 0 : fMatchedTrackIndex = 0;
261 : }
262 :
263 0 : if (reco.fMatchedClusterIndex)
264 : {
265 0 : if (fMatchedClusterIndex)
266 0 : *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
267 : else
268 0 : fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
269 : }
270 : else
271 : {
272 0 : if (fMatchedClusterIndex) delete fMatchedClusterIndex;
273 0 : fMatchedClusterIndex = 0;
274 : }
275 :
276 0 : return *this;
277 0 : }
278 :
279 : ///
280 : /// Destructor.
281 : ///
282 : //_____________________________________
283 : AliEMCALRecoUtils::~AliEMCALRecoUtils()
284 0 : {
285 0 : if (fEMCALRecalibrationFactors)
286 : {
287 0 : fEMCALRecalibrationFactors->Clear();
288 0 : delete fEMCALRecalibrationFactors;
289 : }
290 :
291 0 : if (fEMCALTimeRecalibrationFactors)
292 : {
293 0 : fEMCALTimeRecalibrationFactors->Clear();
294 0 : delete fEMCALTimeRecalibrationFactors;
295 : }
296 :
297 0 : if(fEMCALL1PhaseInTimeRecalibration)
298 : {
299 0 : fEMCALL1PhaseInTimeRecalibration->Clear();
300 0 : delete fEMCALL1PhaseInTimeRecalibration;
301 : }
302 :
303 0 : if (fEMCALBadChannelMap)
304 : {
305 0 : fEMCALBadChannelMap->Clear();
306 0 : delete fEMCALBadChannelMap;
307 : }
308 :
309 0 : delete fMatchedTrackIndex ;
310 0 : delete fMatchedClusterIndex ;
311 0 : delete fResidualEta ;
312 0 : delete fResidualPhi ;
313 0 : delete fPIDUtils ;
314 :
315 0 : InitTrackCuts();
316 0 : }
317 :
318 : ///
319 : /// Reject cell if acceptance criteria not passed (correct cell number, is it bad channel)
320 : /// and calibrate it in energy and time.
321 : ///
322 : /// \param absID: absolute cell ID number
323 : /// \param bc: bunch crossing number
324 : /// \param amp: input cell energy amplitude, output calibrated amplitude
325 : /// \param time: input cell time, output calibrated time
326 : /// \param cells: list of cells
327 : ///
328 : /// \return bool quality of cell, exists or not
329 : ///
330 : //_______________________________________________________________________________
331 : Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(Int_t absID, Int_t bc,
332 : Float_t & amp, Double_t & time,
333 : AliVCaloCells* cells)
334 : {
335 0 : AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
336 :
337 0 : if(!geom)
338 : {
339 0 : AliError("No instance of the geometry is available");
340 0 : return kFALSE;
341 : }
342 :
343 0 : if ( absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules() )
344 0 : return kFALSE;
345 :
346 0 : Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
347 :
348 0 : if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
349 : {
350 : // cell absID does not exist
351 0 : amp=0; time = 1.e9;
352 0 : return kFALSE;
353 : }
354 :
355 0 : geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
356 :
357 : // Do not include bad channels found in analysis,
358 0 : if (IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi))
359 0 : return kFALSE;
360 :
361 : //Recalibrate energy
362 0 : amp = cells->GetCellAmplitude(absID);
363 0 : if (!fCellsRecalibrated && IsRecalibrationOn())
364 0 : amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
365 :
366 : // Recalibrate time
367 0 : time = cells->GetCellTime(absID);
368 0 : time-=fConstantTimeShift*1e-9; // only in case of old Run1 simulation
369 0 : Bool_t isLowGain = !(cells->GetCellHighGain(absID));//HG = false -> LG = true
370 :
371 0 : RecalibrateCellTime(absID,bc,time,isLowGain);
372 :
373 : //Recalibrate time with L1 phase
374 0 : RecalibrateCellTimeL1Phase(imod, bc, time);
375 :
376 : return kTRUE;
377 0 : }
378 :
379 : ///
380 : /// Given the list of AbsId cells of the cluster, get the maximum cell and
381 : /// check if there are fNCellsFromBorder from the calorimeter border.
382 : ///
383 : /// \param geom: AliEMCALGeometry pointer
384 : /// \param cluster: AliVCluster pointer
385 : /// \param cells: list of cells
386 : ///
387 : /// \return bool, true if cluster in defined fiducial region
388 : //_____________________________________________________________________________
389 : Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom,
390 : const AliVCluster* cluster,
391 : AliVCaloCells* cells)
392 : {
393 0 : if (!cluster)
394 : {
395 0 : AliInfo("Cluster pointer null!");
396 0 : return kFALSE;
397 : }
398 :
399 : //If the distance to the border is 0 or negative just exit accept all clusters
400 0 : if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
401 0 : return kTRUE;
402 :
403 0 : Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
404 0 : Bool_t shared = kFALSE;
405 0 : GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
406 :
407 0 : AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
408 : absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
409 :
410 0 : if (absIdMax==-1) return kFALSE;
411 :
412 : //Check if the cell is close to the borders:
413 : Bool_t okrow = kFALSE;
414 : Bool_t okcol = kFALSE;
415 :
416 0 : if (iSM < 0 || iphi < 0 || ieta < 0 ) {
417 0 : AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
418 : iSM,ieta,iphi));
419 0 : return kFALSE; // trick coverity
420 : }
421 :
422 : //Check rows/phi
423 : Int_t iPhiLast = 24;
424 0 : if( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
425 0 : else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
426 :
427 0 : if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
428 :
429 : //Check columns/eta
430 : Int_t iEtaLast = 48;
431 0 : if(!fNoEMCALBorderAtEta0 || geom->IsDCALSM(iSM)) {// conside inner border
432 0 : if( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard ) iEtaLast = iEtaLast*2/3;
433 0 : if(ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
434 : } else {
435 0 : if (iSM%2==0) {
436 0 : if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
437 : } else {
438 0 : if(ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
439 : }
440 : }//eta 0 not checked
441 :
442 0 : AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
443 : fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
444 :
445 0 : if (okcol && okrow) {
446 : //printf("Accept\n");
447 0 : return kTRUE;
448 : } else {
449 : //printf("Reject\n");
450 0 : AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
451 0 : return kFALSE;
452 : }
453 0 : }
454 :
455 : ///
456 : /// Check that in the cluster cells, there is no bad channel of those stored
457 : /// in fEMCALBadChannelMap
458 : ///
459 : /// \param geom: AliEMCALGeometry pointer
460 : /// \param cellsList: list of cells absolute ID in cluster
461 : /// \param nCells: number of cells in cluster
462 : ///
463 : /// \return bool, true if cluster contains a bad channel
464 : ///
465 : //_______________________________________________________________________________
466 : Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom,
467 : const UShort_t* cellList,
468 : Int_t nCells)
469 : {
470 0 : if (!fRemoveBadChannels) return kFALSE;
471 0 : if (!fEMCALBadChannelMap) return kFALSE;
472 :
473 0 : Int_t icol = -1;
474 0 : Int_t irow = -1;
475 0 : Int_t imod = -1;
476 0 : for (Int_t iCell = 0; iCell<nCells; iCell++)
477 : {
478 : //Get the column and row
479 0 : Int_t iTower = -1, iIphi = -1, iIeta = -1;
480 0 : geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
481 :
482 0 : if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
483 :
484 0 : geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
485 :
486 0 : if (GetEMCALChannelStatus(imod, icol, irow))
487 : {
488 0 : AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
489 0 : return kTRUE;
490 : }
491 0 : }// cell cluster loop
492 :
493 0 : return kFALSE;
494 0 : }
495 :
496 : ///
497 : /// Calculate the energy in the cross around the energy of a given cell.
498 : /// Used in exotic clusters/cells rejection.
499 : ///
500 : /// \param absID: controlled cell absolute ID number
501 : /// \param tcell: time of cell under control
502 : /// \param cells: full list of cells
503 : /// \param bc: bunch crossing number
504 : ///
505 : /// \return float E_cross
506 : ///
507 : //___________________________________________________________________________
508 : Float_t AliEMCALRecoUtils::GetECross(Int_t absID, Double_t tcell,
509 : AliVCaloCells* cells, Int_t bc)
510 : {
511 0 : AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
512 :
513 0 : if(!geom)
514 : {
515 0 : AliError("No instance of the geometry is available");
516 0 : return -1;
517 : }
518 :
519 0 : Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
520 0 : geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
521 0 : geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
522 :
523 : // Get close cells index, energy and time, not in corners
524 :
525 : Int_t absID1 = -1;
526 : Int_t absID2 = -1;
527 :
528 0 : if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
529 0 : if ( iphi > 0 ) absID2 = geom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
530 :
531 : // In case of cell in eta = 0 border, depending on SM shift the cross cell index
532 :
533 : Int_t absID3 = -1;
534 : Int_t absID4 = -1;
535 :
536 0 : if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
537 : {
538 0 : absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
539 0 : absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
540 0 : }
541 0 : else if ( ieta == 0 && imod%2 )
542 : {
543 0 : absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
544 0 : absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
545 0 : }
546 : else
547 : {
548 0 : if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
549 0 : absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
550 0 : if ( ieta > 0 )
551 0 : absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
552 : }
553 :
554 : //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
555 :
556 0 : Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
557 0 : Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
558 :
559 0 : AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
560 0 : AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
561 0 : AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
562 0 : AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
563 :
564 0 : if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
565 0 : if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
566 0 : if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
567 0 : if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
568 :
569 0 : return ecell1+ecell2+ecell3+ecell4;
570 0 : }
571 :
572 : ///
573 : /// Look to cell neighbourhood and reject if it seems exotic
574 : /// Do before recalibrating the cells.
575 : ///
576 : /// \param absID: controlled cell absolute ID number
577 : /// \param tcell: time of cell under control
578 : /// \param cells: full list of cells
579 : /// \param bc: bunch crossing number
580 : ///
581 : /// \return bool true if cell is found exotic
582 : ///
583 : //_____________________________________________________________________________________________
584 : Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
585 : {
586 0 : if (!fRejectExoticCells) return kFALSE;
587 :
588 0 : Float_t ecell = 0;
589 0 : Double_t tcell = 0;
590 0 : Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
591 :
592 0 : if (!accept) return kTRUE; // reject this cell
593 :
594 0 : if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
595 :
596 0 : Float_t eCross = GetECross(absID,tcell,cells,bc);
597 :
598 0 : if (1-eCross/ecell > fExoticCellFraction)
599 : {
600 0 : AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
601 : absID,ecell,eCross,1-eCross/ecell));
602 0 : return kTRUE;
603 : }
604 :
605 0 : return kFALSE;
606 0 : }
607 :
608 : ///
609 : /// Check if the cluster highest energy tower is exotic.
610 : ///
611 : /// \param cluster: pointer to AliVCluster
612 : /// \param cells: full list of cells
613 : /// \param bc: bunch crossing number
614 : ///
615 : /// \param cluster: pointer to AliVCluster
616 : ///
617 : //___________________________________________________________________
618 : Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
619 : AliVCaloCells *cells,
620 : Int_t bc)
621 : {
622 0 : if (!cluster)
623 : {
624 0 : AliInfo("Cluster pointer null!");
625 0 : return kFALSE;
626 : }
627 :
628 0 : if (!fRejectExoticCluster) return kFALSE;
629 :
630 : // Get highest energy tower
631 0 : AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
632 :
633 0 : if(!geom)
634 : {
635 0 : AliError("No instance of the geometry is available");
636 0 : return kFALSE;
637 : }
638 :
639 0 : Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
640 0 : Bool_t shared = kFALSE;
641 0 : GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
642 :
643 0 : return IsExoticCell(absId,cells,bc);
644 0 : }
645 :
646 : ///
647 : /// In case of MC analysis, smear energy to match resolution/calibration in real data
648 : /// (old, in principle not needed anymore).
649 : ///
650 : /// \param cluster: pointer to AliVCluster
651 : ///
652 : /// \return float with smeared energy
653 : ///
654 : //_______________________________________________________________________
655 : Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
656 : {
657 0 : if (!cluster)
658 : {
659 0 : AliInfo("Cluster pointer null!");
660 0 : return 0;
661 : }
662 :
663 0 : Float_t energy = cluster->E() ;
664 : Float_t rdmEnergy = energy ;
665 :
666 0 : if (fSmearClusterEnergy)
667 : {
668 0 : rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
669 0 : fSmearClusterParam[1] * energy +
670 0 : fSmearClusterParam[2] );
671 0 : AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
672 : }
673 :
674 : return rdmEnergy;
675 0 : }
676 :
677 : ///
678 : /// Correct cluster energy from non linearity functions, defined in enum NonlinearityFunctions
679 : /// Recomended for data kBeamTestCorrectedv3 and for simulation kPi0MCv3
680 : ///
681 : /// \param cluster: pointer to AliVCluster
682 : ///
683 : /// \return float with corrected cluster energy
684 : ///
685 : //____________________________________________________________________________
686 : Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
687 : {
688 0 : if (!cluster)
689 : {
690 0 : AliInfo("Cluster pointer null!");
691 0 : return 0;
692 : }
693 :
694 0 : Float_t energy = cluster->E();
695 0 : if (energy < 0.100)
696 : {
697 : // Clusters with less than 100 MeV or negative are not possible
698 0 : AliInfo(Form("Too low cluster energy!, E = %f < 0.100 GeV",energy));
699 0 : return 0;
700 : }
701 0 : else if(energy > 300.)
702 : {
703 0 : AliInfo(Form("Too high cluster energy!, E = %f GeV, do not apply non linearity",energy));
704 0 : return energy;
705 : }
706 :
707 0 : switch (fNonLinearityFunction)
708 : {
709 : case kPi0MC:
710 : {
711 : //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
712 : //fNonLinearityParams[0] = 1.014;
713 : //fNonLinearityParams[1] =-0.03329;
714 : //fNonLinearityParams[2] =-0.3853;
715 : //fNonLinearityParams[3] = 0.5423;
716 : //fNonLinearityParams[4] =-0.4335;
717 0 : energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
718 0 : ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
719 0 : exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
720 0 : break;
721 : }
722 :
723 : case kPi0MCv2:
724 : {
725 : //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
726 : //fNonLinearityParams[0] = 3.11111e-02;
727 : //fNonLinearityParams[1] =-5.71666e-02;
728 : //fNonLinearityParams[2] = 5.67995e-01;
729 :
730 0 : energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
731 0 : break;
732 : }
733 :
734 : case kPi0MCv3:
735 : {
736 : //Same as beam test corrected, change parameters
737 : //fNonLinearityParams[0] = 9.81039e-01
738 : //fNonLinearityParams[1] = 1.13508e-01;
739 : //fNonLinearityParams[2] = 1.00173e+00;
740 : //fNonLinearityParams[3] = 9.67998e-02;
741 : //fNonLinearityParams[4] = 2.19381e+02;
742 : //fNonLinearityParams[5] = 6.31604e+01;
743 : //fNonLinearityParams[6] = 1;
744 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
745 :
746 0 : break;
747 : }
748 :
749 :
750 : case kPi0GammaGamma:
751 : {
752 : //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
753 : //fNonLinearityParams[0] = 1.04;
754 : //fNonLinearityParams[1] = -0.1445;
755 : //fNonLinearityParams[2] = 1.046;
756 0 : energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
757 0 : break;
758 : }
759 :
760 : case kPi0GammaConversion:
761 : {
762 : //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
763 : //fNonLinearityParams[0] = 0.139393/0.1349766;
764 : //fNonLinearityParams[1] = 0.0566186;
765 : //fNonLinearityParams[2] = 0.982133;
766 0 : energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
767 :
768 0 : break;
769 : }
770 :
771 : case kBeamTest:
772 : {
773 : //From beam test, Alexei's results, for different ZS thresholds
774 : // th=30 MeV; th = 45 MeV; th = 75 MeV
775 : //fNonLinearityParams[0] = 1.007; 1.003; 1.002
776 : //fNonLinearityParams[1] = 0.894; 0.719; 0.797
777 : //fNonLinearityParams[2] = 0.246; 0.334; 0.358
778 : //Rescale the param[0] with 1.03
779 0 : energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
780 :
781 0 : break;
782 : }
783 :
784 : case kBeamTestCorrected:
785 : {
786 : // From beam test, corrected for material between beam and EMCAL
787 : //fNonLinearityParams[0] = 0.99078
788 : //fNonLinearityParams[1] = 0.161499;
789 : //fNonLinearityParams[2] = 0.655166;
790 : //fNonLinearityParams[3] = 0.134101;
791 : //fNonLinearityParams[4] = 163.282;
792 : //fNonLinearityParams[5] = 23.6904;
793 : //fNonLinearityParams[6] = 0.978;
794 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
795 :
796 0 : break;
797 : }
798 :
799 : case kBeamTestCorrectedv2:
800 : {
801 : // From beam test, corrected for material between beam and EMCAL
802 : // Different function to kBeamTestCorrected
803 : //fNonLinearityParams[0] = 0.983504;
804 : //fNonLinearityParams[1] = 0.210106;
805 : //fNonLinearityParams[2] = 0.897274;
806 : //fNonLinearityParams[3] = 0.0829064;
807 : //fNonLinearityParams[4] = 152.299;
808 : //fNonLinearityParams[5] = 31.5028;
809 : //fNonLinearityParams[6] = 0.968;
810 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
811 :
812 0 : break;
813 : }
814 :
815 : case kBeamTestCorrectedv3:
816 : {
817 : // Same function as kBeamTestCorrectedv2, different default parametrization.
818 : //fNonLinearityParams[0] = 0.976941;
819 : //fNonLinearityParams[1] = 0.162310;
820 : //fNonLinearityParams[2] = 1.08689;
821 : //fNonLinearityParams[3] = 0.0819592;
822 : //fNonLinearityParams[4] = 152.338;
823 : //fNonLinearityParams[5] = 30.9594;
824 : //fNonLinearityParams[6] = 0.9615;
825 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
826 :
827 0 : break;
828 : }
829 :
830 : case kSDMv5:
831 : {
832 : //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
833 : //fNonLinearityParams[0] = 1.0;
834 : //fNonLinearityParams[1] = 6.64778e-02;
835 : //fNonLinearityParams[2] = 1.570;
836 : //fNonLinearityParams[3] = 9.67998e-02;
837 : //fNonLinearityParams[4] = 2.19381e+02;
838 : //fNonLinearityParams[5] = 6.31604e+01;
839 : //fNonLinearityParams[6] = 1.01286;
840 0 : 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));
841 :
842 0 : break;
843 : }
844 :
845 : case kPi0MCv5:
846 : {
847 : //Based on comparing MC truth information to the reconstructed energy of clusters.
848 : //fNonLinearityParams[0] = 1.0;
849 : //fNonLinearityParams[1] = 6.64778e-02;
850 : //fNonLinearityParams[2] = 1.570;
851 : //fNonLinearityParams[3] = 9.67998e-02;
852 : //fNonLinearityParams[4] = 2.19381e+02;
853 : //fNonLinearityParams[5] = 6.31604e+01;
854 : //fNonLinearityParams[6] = 1.01286;
855 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
856 :
857 0 : break;
858 : }
859 :
860 : case kSDMv6:
861 : {
862 : //Based on fit to the MC/data using kNoCorrection on the data
863 : // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
864 : // - parameters constrained by the test beam data as well
865 : // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
866 : //fNonLinearityParams[0] = 1.0;
867 : //fNonLinearityParams[1] = 0.237767;
868 : //fNonLinearityParams[2] = 0.651203;
869 : //fNonLinearityParams[3] = 0.183741;
870 : //fNonLinearityParams[4] = 155.427;
871 : //fNonLinearityParams[5] = 17.0335;
872 : //fNonLinearityParams[6] = 0.987054;
873 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
874 :
875 0 : break;
876 : }
877 :
878 : case kPi0MCv6:
879 : {
880 : //Based on comparing MC truth information to the reconstructed energy of clusters.
881 : // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
882 : //fNonLinearityParams[0] = 1.0;
883 : //fNonLinearityParams[1] = 0.0797873;
884 : //fNonLinearityParams[2] = 1.68322;
885 : //fNonLinearityParams[3] = 0.0806098;
886 : //fNonLinearityParams[4] = 244.586;
887 : //fNonLinearityParams[5] = 116.938;
888 : //fNonLinearityParams[6] = 1.00437;
889 0 : energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
890 :
891 0 : break;
892 : }
893 :
894 : case kNoCorrection:
895 0 : AliDebug(2,"No correction on the energy\n");
896 : break;
897 :
898 : }
899 :
900 0 : return energy;
901 0 : }
902 :
903 : ///
904 : /// Initialising non Linearity Parameters for the different
905 : /// parametrizations available, defined in enum NonlinearityFunctions
906 : ///
907 : //__________________________________________________
908 : void AliEMCALRecoUtils::InitNonLinearityParam()
909 : {
910 0 : if (fNonLinearityFunction == kPi0MC) {
911 0 : fNonLinearityParams[0] = 1.014;
912 0 : fNonLinearityParams[1] = -0.03329;
913 0 : fNonLinearityParams[2] = -0.3853;
914 0 : fNonLinearityParams[3] = 0.5423;
915 0 : fNonLinearityParams[4] = -0.4335;
916 0 : }
917 :
918 0 : if (fNonLinearityFunction == kPi0MCv2) {
919 0 : fNonLinearityParams[0] = 3.11111e-02;
920 0 : fNonLinearityParams[1] =-5.71666e-02;
921 0 : fNonLinearityParams[2] = 5.67995e-01;
922 0 : }
923 :
924 0 : if (fNonLinearityFunction == kPi0MCv3) {
925 0 : fNonLinearityParams[0] = 9.81039e-01;
926 0 : fNonLinearityParams[1] = 1.13508e-01;
927 0 : fNonLinearityParams[2] = 1.00173e+00;
928 0 : fNonLinearityParams[3] = 9.67998e-02;
929 0 : fNonLinearityParams[4] = 2.19381e+02;
930 0 : fNonLinearityParams[5] = 6.31604e+01;
931 0 : fNonLinearityParams[6] = 1;
932 0 : }
933 :
934 0 : if (fNonLinearityFunction == kPi0GammaGamma) {
935 0 : fNonLinearityParams[0] = 1.04;
936 0 : fNonLinearityParams[1] = -0.1445;
937 0 : fNonLinearityParams[2] = 1.046;
938 0 : }
939 :
940 0 : if (fNonLinearityFunction == kPi0GammaConversion) {
941 0 : fNonLinearityParams[0] = 0.139393;
942 0 : fNonLinearityParams[1] = 0.0566186;
943 0 : fNonLinearityParams[2] = 0.982133;
944 0 : }
945 :
946 0 : if (fNonLinearityFunction == kBeamTest) {
947 0 : if (fNonLinearThreshold == 30) {
948 0 : fNonLinearityParams[0] = 1.007;
949 0 : fNonLinearityParams[1] = 0.894;
950 0 : fNonLinearityParams[2] = 0.246;
951 0 : }
952 0 : if (fNonLinearThreshold == 45) {
953 0 : fNonLinearityParams[0] = 1.003;
954 0 : fNonLinearityParams[1] = 0.719;
955 0 : fNonLinearityParams[2] = 0.334;
956 0 : }
957 0 : if (fNonLinearThreshold == 75) {
958 0 : fNonLinearityParams[0] = 1.002;
959 0 : fNonLinearityParams[1] = 0.797;
960 0 : fNonLinearityParams[2] = 0.358;
961 0 : }
962 : }
963 :
964 0 : if (fNonLinearityFunction == kBeamTestCorrected) {
965 0 : fNonLinearityParams[0] = 0.99078;
966 0 : fNonLinearityParams[1] = 0.161499;
967 0 : fNonLinearityParams[2] = 0.655166;
968 0 : fNonLinearityParams[3] = 0.134101;
969 0 : fNonLinearityParams[4] = 163.282;
970 0 : fNonLinearityParams[5] = 23.6904;
971 0 : fNonLinearityParams[6] = 0.978;
972 0 : }
973 :
974 0 : if (fNonLinearityFunction == kBeamTestCorrectedv2) {
975 : // Parameters until November 2015, use now kBeamTestCorrectedv3
976 0 : fNonLinearityParams[0] = 0.983504;
977 0 : fNonLinearityParams[1] = 0.210106;
978 0 : fNonLinearityParams[2] = 0.897274;
979 0 : fNonLinearityParams[3] = 0.0829064;
980 0 : fNonLinearityParams[4] = 152.299;
981 0 : fNonLinearityParams[5] = 31.5028;
982 0 : fNonLinearityParams[6] = 0.968;
983 0 : }
984 :
985 0 : if (fNonLinearityFunction == kBeamTestCorrectedv3) {
986 :
987 : // New parametrization of kBeamTestCorrectedv2
988 : // excluding point at 0.5 GeV from Beam Test Data
989 : // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
990 :
991 0 : fNonLinearityParams[0] = 0.976941;
992 0 : fNonLinearityParams[1] = 0.162310;
993 0 : fNonLinearityParams[2] = 1.08689;
994 0 : fNonLinearityParams[3] = 0.0819592;
995 0 : fNonLinearityParams[4] = 152.338;
996 0 : fNonLinearityParams[5] = 30.9594;
997 0 : fNonLinearityParams[6] = 0.9615;
998 0 : }
999 :
1000 0 : if (fNonLinearityFunction == kSDMv5) {
1001 0 : fNonLinearityParams[0] = 1.0;
1002 0 : fNonLinearityParams[1] = 6.64778e-02;
1003 0 : fNonLinearityParams[2] = 1.570;
1004 0 : fNonLinearityParams[3] = 9.67998e-02;
1005 0 : fNonLinearityParams[4] = 2.19381e+02;
1006 0 : fNonLinearityParams[5] = 6.31604e+01;
1007 0 : fNonLinearityParams[6] = 1.01286;
1008 0 : }
1009 :
1010 0 : if (fNonLinearityFunction == kPi0MCv5) {
1011 0 : fNonLinearityParams[0] = 1.0;
1012 0 : fNonLinearityParams[1] = 6.64778e-02;
1013 0 : fNonLinearityParams[2] = 1.570;
1014 0 : fNonLinearityParams[3] = 9.67998e-02;
1015 0 : fNonLinearityParams[4] = 2.19381e+02;
1016 0 : fNonLinearityParams[5] = 6.31604e+01;
1017 0 : fNonLinearityParams[6] = 1.01286;
1018 0 : }
1019 :
1020 0 : if (fNonLinearityFunction == kSDMv6) {
1021 0 : fNonLinearityParams[0] = 1.0;
1022 0 : fNonLinearityParams[1] = 0.237767;
1023 0 : fNonLinearityParams[2] = 0.651203;
1024 0 : fNonLinearityParams[3] = 0.183741;
1025 0 : fNonLinearityParams[4] = 155.427;
1026 0 : fNonLinearityParams[5] = 17.0335;
1027 0 : fNonLinearityParams[6] = 0.987054;
1028 0 : }
1029 :
1030 0 : if (fNonLinearityFunction == kPi0MCv6) {
1031 0 : fNonLinearityParams[0] = 1.0;
1032 0 : fNonLinearityParams[1] = 0.0797873;
1033 0 : fNonLinearityParams[2] = 1.68322;
1034 0 : fNonLinearityParams[3] = 0.0806098;
1035 0 : fNonLinearityParams[4] = 244.586;
1036 0 : fNonLinearityParams[5] = 116.938;
1037 0 : fNonLinearityParams[6] = 1.00437;
1038 0 : }
1039 0 : }
1040 :
1041 : ///
1042 : /// Calculate shower depth for a given cluster energy and particle type.
1043 : /// Needed to calculate cluster position.
1044 : ///
1045 : /// \param energy: cluster energy
1046 : /// \param iParticle: particle assumption defined in enum ParticleType
1047 : /// \param iSM: supermodule number
1048 : ///
1049 : //_________________________________________________________
1050 : Float_t AliEMCALRecoUtils::GetDepth(Float_t energy,
1051 : Int_t iParticle,
1052 : Int_t iSM) const
1053 : {
1054 : // parameters
1055 : Float_t x0 = 1.31;
1056 : Float_t ecr = 8;
1057 : Float_t depth = 0;
1058 0 : Float_t arg = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV
1059 :
1060 0 : switch ( iParticle )
1061 : {
1062 : case kPhoton:
1063 0 : if (arg < 1)
1064 0 : depth = 0;
1065 : else
1066 0 : depth = x0 * (TMath::Log(arg) + 0.5);
1067 : break;
1068 :
1069 : case kElectron:
1070 0 : if (arg < 1)
1071 0 : depth = 0;
1072 : else
1073 0 : depth = x0 * (TMath::Log(arg) - 0.5);
1074 : break;
1075 :
1076 : case kHadron:
1077 : // hadron
1078 : // boxes anc. here
1079 0 : if (gGeoManager)
1080 : {
1081 0 : gGeoManager->cd("ALIC_1/XEN1_1");
1082 0 : TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
1083 0 : TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
1084 0 : if (geoSM)
1085 : {
1086 0 : TGeoVolume *geoSMVol = geoSM->GetVolume();
1087 0 : TGeoShape *geoSMShape = geoSMVol->GetShape();
1088 0 : TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
1089 0 : if (geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
1090 0 : else AliFatal("Null GEANT box");
1091 0 : }
1092 0 : else AliFatal("NULL GEANT node matrix");
1093 0 : }
1094 : else
1095 : {//electron
1096 0 : if (arg < 1)
1097 0 : depth = 0;
1098 : else
1099 0 : depth = x0 * (TMath::Log(arg) - 0.5);
1100 : }
1101 :
1102 : break;
1103 :
1104 : default://photon
1105 0 : if (arg < 1)
1106 0 : depth = 0;
1107 : else
1108 0 : depth = x0 * (TMath::Log(arg) + 0.5);
1109 : }
1110 :
1111 0 : return depth;
1112 : }
1113 :
1114 : ///
1115 : /// For a given CaloCluster gets the absId of the cell
1116 : /// with maximum energy deposit.
1117 : ///
1118 : /// \param geom: AliEMCALGeometry pointer
1119 : /// \param cells: full list of cells
1120 : /// \param clu: pointer to AliVCluster
1121 : /// \param absId: absolute id number of cell with highest energy in cluster
1122 : /// \param iSupmod: supermodule number of cell with highest energy in cluster
1123 : /// \param ieta: column number of cell with highest energy in cluster
1124 : /// \param iphi: row number of cell with highest energy in cluster
1125 : /// \param shared: cluster is shared between 2 supermodules
1126 : ///
1127 : //____________________________________________________________________
1128 : void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1129 : AliVCaloCells* cells,
1130 : const AliVCluster* clu,
1131 : Int_t & absId,
1132 : Int_t & iSupMod,
1133 : Int_t & ieta,
1134 : Int_t & iphi,
1135 : Bool_t & shared)
1136 : {
1137 : Double_t eMax = -1.;
1138 : Double_t eCell = -1.;
1139 : Float_t fraction = 1.;
1140 : Float_t recalFactor = 1.;
1141 : Int_t cellAbsId = -1 ;
1142 :
1143 0 : Int_t iTower = -1;
1144 0 : Int_t iIphi = -1;
1145 0 : Int_t iIeta = -1;
1146 : Int_t iSupMod0= -1;
1147 :
1148 0 : if (!clu)
1149 : {
1150 0 : AliInfo("Cluster pointer null!");
1151 0 : absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1152 0 : return;
1153 : }
1154 :
1155 0 : for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1156 : {
1157 0 : cellAbsId = clu->GetCellAbsId(iDig);
1158 0 : fraction = clu->GetCellAmplitudeFraction(iDig);
1159 : //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1160 0 : if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1161 :
1162 0 : geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1163 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1164 :
1165 0 : if (iDig==0)
1166 : {
1167 0 : iSupMod0=iSupMod;
1168 0 : } else if (iSupMod0!=iSupMod)
1169 : {
1170 0 : shared = kTRUE;
1171 : //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1172 0 : }
1173 :
1174 0 : if (!fCellsRecalibrated && IsRecalibrationOn())
1175 0 : recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1176 :
1177 0 : eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1178 : //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1179 0 : if (eCell > eMax)
1180 : {
1181 : eMax = eCell;
1182 0 : absId = cellAbsId;
1183 : //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1184 0 : }
1185 : }// cell loop
1186 :
1187 : //Get from the absid the supermodule, tower and eta/phi numbers
1188 0 : geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1189 : //Gives SuperModule and Tower numbers
1190 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1191 : //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1192 : //printf("Max end---\n");
1193 0 : }
1194 :
1195 : ///
1196 : /// Initialize data members with default values
1197 : ///
1198 : //______________________________________
1199 : void AliEMCALRecoUtils::InitParameters()
1200 : {
1201 0 : fParticleType = kPhoton;
1202 0 : fPosAlgo = kUnchanged;
1203 0 : fW0 = 4.5;
1204 :
1205 0 : fNonLinearityFunction = kNoCorrection;
1206 0 : fNonLinearThreshold = 30;
1207 :
1208 0 : fExoticCellFraction = 0.97;
1209 0 : fExoticCellDiffTime = 1e6;
1210 0 : fExoticCellMinAmplitude = 4.0;
1211 :
1212 0 : fAODFilterMask = 128;
1213 0 : fAODHybridTracks = kFALSE;
1214 0 : fAODTPCOnlyTracks = kTRUE;
1215 :
1216 0 : fCutEtaPhiSum = kTRUE;
1217 0 : fCutEtaPhiSeparate = kFALSE;
1218 :
1219 0 : fCutR = 0.05;
1220 0 : fCutEta = 0.025;
1221 0 : fCutPhi = 0.05;
1222 :
1223 0 : fClusterWindow = 100;
1224 0 : fMass = 0.139;
1225 :
1226 0 : fStepSurface = 20.;
1227 0 : fStepCluster = 5.;
1228 0 : fTrackCutsType = kLooseCut;
1229 :
1230 0 : fCutMinTrackPt = 0;
1231 0 : fCutMinNClusterTPC = -1;
1232 0 : fCutMinNClusterITS = -1;
1233 :
1234 0 : fCutMaxChi2PerClusterTPC = 1e10;
1235 0 : fCutMaxChi2PerClusterITS = 1e10;
1236 :
1237 0 : fCutRequireTPCRefit = kFALSE;
1238 0 : fCutRequireITSRefit = kFALSE;
1239 0 : fCutAcceptKinkDaughters = kFALSE;
1240 :
1241 0 : fCutMaxDCAToVertexXY = 1e10;
1242 0 : fCutMaxDCAToVertexZ = 1e10;
1243 0 : fCutDCAToVertex2D = kFALSE;
1244 :
1245 0 : fCutRequireITSStandAlone = kFALSE; //MARCEL
1246 0 : fCutRequireITSpureSA = kFALSE; //Marcel
1247 :
1248 : // Misalignment matrices
1249 0 : for (Int_t i = 0; i < 15 ; i++)
1250 : {
1251 0 : fMisalTransShift[i] = 0.;
1252 0 : fMisalRotShift[i] = 0.;
1253 : }
1254 :
1255 : // Non linearity
1256 0 : for (Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
1257 :
1258 : // For kBeamTestCorrectedv2 case, but default is no correction
1259 0 : fNonLinearityParams[0] = 0.983504;
1260 0 : fNonLinearityParams[1] = 0.210106;
1261 0 : fNonLinearityParams[2] = 0.897274;
1262 0 : fNonLinearityParams[3] = 0.0829064;
1263 0 : fNonLinearityParams[4] = 152.299;
1264 0 : fNonLinearityParams[5] = 31.5028;
1265 0 : fNonLinearityParams[6] = 0.968;
1266 :
1267 : // Cluster energy smearing
1268 0 : fSmearClusterEnergy = kFALSE;
1269 0 : fSmearClusterParam[0] = 0.07; // * sqrt E term
1270 0 : fSmearClusterParam[1] = 0.00; // * E term
1271 0 : fSmearClusterParam[2] = 0.00; // constant
1272 0 : }
1273 :
1274 : ///
1275 : /// Init EMCAL energy calibration factors container
1276 : ///
1277 : //_____________________________________________________
1278 : void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
1279 : {
1280 0 : AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1281 :
1282 : // In order to avoid rewriting the same histograms
1283 0 : Bool_t oldStatus = TH1::AddDirectoryStatus();
1284 0 : TH1::AddDirectory(kFALSE);
1285 :
1286 0 : fEMCALRecalibrationFactors = new TObjArray(22);
1287 0 : for (int i = 0; i < 22; i++)
1288 0 : fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1289 0 : Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1290 : //Init the histograms with 1
1291 0 : for (Int_t sm = 0; sm < 22; sm++)
1292 : {
1293 0 : for (Int_t i = 0; i < 48; i++)
1294 : {
1295 0 : for (Int_t j = 0; j < 24; j++)
1296 : {
1297 0 : SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1298 : }
1299 : }
1300 : }
1301 :
1302 0 : fEMCALRecalibrationFactors->SetOwner(kTRUE);
1303 0 : fEMCALRecalibrationFactors->Compress();
1304 :
1305 : // In order to avoid rewriting the same histograms
1306 0 : TH1::AddDirectory(oldStatus);
1307 0 : }
1308 :
1309 : ///
1310 : /// Init EMCAL time calibration shifts container
1311 : ///
1312 : //_________________________________________________________
1313 : void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1314 : {
1315 0 : AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1316 :
1317 : // In order to avoid rewriting the same histograms
1318 0 : Bool_t oldStatus = TH1::AddDirectoryStatus();
1319 0 : TH1::AddDirectory(kFALSE);
1320 :
1321 0 : if(fLowGain) fEMCALTimeRecalibrationFactors = new TObjArray(8);
1322 0 : else fEMCALTimeRecalibrationFactors = new TObjArray(4);
1323 :
1324 0 : for (int i = 0; i < 4; i++)
1325 0 : fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1326 0 : Form("hAllTimeAvBC%d",i),
1327 : 48*24*22,0.,48*24*22) );
1328 : // Init the histograms with 0
1329 0 : for (Int_t iBC = 0; iBC < 4; iBC++)
1330 : {
1331 0 : for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1332 0 : SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1333 : }
1334 :
1335 0 : if(fLowGain) {
1336 0 : for (int iBC = 0; iBC < 4; iBC++) {
1337 0 : fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1338 0 : Form("hAllTimeAvLGBC%d",iBC),
1339 : 48*24*22,0.,48*24*22) );
1340 0 : for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1341 0 : SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1342 : }
1343 0 : }
1344 :
1345 :
1346 :
1347 0 : fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1348 0 : fEMCALTimeRecalibrationFactors->Compress();
1349 :
1350 : // In order to avoid rewriting the same histograms
1351 0 : TH1::AddDirectory(oldStatus);
1352 0 : }
1353 :
1354 : ///
1355 : /// Init EMCAL bad channels map container
1356 : ///
1357 : //____________________________________________________
1358 : void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1359 : {
1360 0 : AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1361 :
1362 : // In order to avoid rewriting the same histograms
1363 0 : Bool_t oldStatus = TH1::AddDirectoryStatus();
1364 0 : TH1::AddDirectory(kFALSE);
1365 :
1366 0 : fEMCALBadChannelMap = new TObjArray(22);
1367 : //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1368 :
1369 0 : for (int i = 0; i < 22; i++)
1370 0 : fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1371 :
1372 0 : fEMCALBadChannelMap->SetOwner(kTRUE);
1373 0 : fEMCALBadChannelMap->Compress();
1374 :
1375 : // In order to avoid rewriting the same histograms
1376 0 : TH1::AddDirectory(oldStatus);
1377 0 : }
1378 :
1379 : ///
1380 : /// Init EMCAL L1 phase shifts container
1381 : ///
1382 : //___________________________________________________________
1383 : void AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibration()
1384 : {
1385 0 : AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
1386 :
1387 : // In order to avoid rewriting the same histograms
1388 0 : Bool_t oldStatus = TH1::AddDirectoryStatus();
1389 0 : TH1::AddDirectory(kFALSE);
1390 :
1391 0 : fEMCALL1PhaseInTimeRecalibration = new TObjArray(1);
1392 :
1393 0 : fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
1394 :
1395 0 : for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
1396 0 : SetEMCALL1PhaseInTimeRecalibrationForSM(i,0);
1397 :
1398 0 : fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
1399 0 : fEMCALL1PhaseInTimeRecalibration->Compress();
1400 :
1401 : // In order to avoid rewriting the same histograms
1402 0 : TH1::AddDirectory(oldStatus);
1403 0 : }
1404 :
1405 : ///
1406 : /// Recalibrate the cluster energy and time, considering the recalibration map
1407 : /// and the time and energy of the cells that compose the cluster.
1408 : ///
1409 : /// \param geom: pointer to geometry
1410 : /// \param cluster: pointer to cluster
1411 : /// \param cells: list of cells
1412 : /// \param bc: bunch crossing number returned by esdevent->GetBunchCrossNumber()
1413 : ///
1414 : //____________________________________________________________________________
1415 : void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1416 : AliVCluster * cluster,
1417 : AliVCaloCells * cells,
1418 : Int_t bc)
1419 : {
1420 0 : if (!cluster)
1421 : {
1422 0 : AliInfo("Cluster pointer null!");
1423 0 : return;
1424 : }
1425 :
1426 : // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1427 0 : UShort_t * index = cluster->GetCellsAbsId() ;
1428 0 : Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1429 0 : Int_t ncells = cluster->GetNCells();
1430 :
1431 : // Initialize some used variables
1432 : Float_t energy = 0;
1433 : Int_t absId =-1;
1434 0 : Int_t icol =-1, irow =-1, imod=1;
1435 : Float_t factor = 1, frac = 0;
1436 : Int_t absIdMax = -1;
1437 : Float_t emax = 0;
1438 :
1439 : // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1440 0 : for (Int_t icell = 0; icell < ncells; icell++)
1441 : {
1442 0 : absId = index[icell];
1443 0 : frac = fraction[icell];
1444 0 : if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1445 :
1446 0 : if (!fCellsRecalibrated && IsRecalibrationOn())
1447 : {
1448 : // Energy
1449 0 : Int_t iTower = -1, iIphi = -1, iIeta = -1;
1450 0 : geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1451 0 : if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1452 0 : continue;
1453 0 : geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1454 0 : factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1455 :
1456 0 : AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1457 : imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1458 :
1459 0 : }
1460 :
1461 0 : energy += cells->GetCellAmplitude(absId)*factor*frac;
1462 :
1463 0 : if (emax < cells->GetCellAmplitude(absId)*factor*frac)
1464 : {
1465 0 : emax = cells->GetCellAmplitude(absId)*factor*frac;
1466 : absIdMax = absId;
1467 0 : }
1468 : }
1469 :
1470 0 : AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1471 :
1472 0 : cluster->SetE(energy);
1473 :
1474 : // Recalculate time of cluster
1475 0 : Double_t timeorg = cluster->GetTOF();
1476 0 : Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
1477 :
1478 0 : Double_t time = cells->GetCellTime(absIdMax);
1479 0 : if (!fCellsRecalibrated && IsTimeRecalibrationOn())
1480 0 : RecalibrateCellTime(absIdMax,bc,time,isLowGain);
1481 0 : time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
1482 :
1483 : // Recalibrate time with L1 phase
1484 0 : if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn())
1485 0 : RecalibrateCellTimeL1Phase(imod, bc, time);
1486 :
1487 0 : cluster->SetTOF(time);
1488 :
1489 0 : AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1490 0 : }
1491 :
1492 : ///
1493 : /// Recalibrate all the cells time and energy, considering the recalibration map and
1494 : /// the energy and time of each cell.
1495 : ///
1496 : /// \param cells: list of cells
1497 : /// \param bc: bunch crossing number returned by esdevent->GetBunchCrossNumber()
1498 : ///
1499 : //_______________________________________________________________________
1500 : void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
1501 : {
1502 0 : if (!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn())
1503 : return;
1504 :
1505 0 : if (!cells)
1506 : {
1507 0 : AliInfo("Cells pointer null!");
1508 0 : return;
1509 : }
1510 :
1511 0 : Short_t absId =-1;
1512 : Bool_t accept = kFALSE;
1513 0 : Float_t ecell = 0;
1514 0 : Double_t tcell = 0;
1515 0 : Double_t ecellin = 0;
1516 0 : Double_t tcellin = 0;
1517 0 : Int_t mclabel = -1;
1518 0 : Double_t efrac = 0;
1519 :
1520 0 : Int_t nEMcell = cells->GetNumberOfCells() ;
1521 0 : for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1522 : {
1523 0 : cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1524 :
1525 0 : accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1526 0 : if (!accept)
1527 : {
1528 0 : ecell = 0;
1529 0 : tcell = -1;
1530 0 : }
1531 :
1532 : // Set new values
1533 0 : cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1534 : }
1535 :
1536 0 : fCellsRecalibrated = kTRUE;
1537 0 : }
1538 :
1539 : ///
1540 : /// Recalibrate time of cell from AbsID number considering cell calibration map
1541 : ///
1542 : /// \param absID: cell absolute ID number
1543 : /// \param bc: bunch crossing number returned by esdevent->GetBunchCrossNumber()
1544 : /// \param celltime: cell time to be returned calibrated
1545 : ///
1546 : //_______________________________________________________________________________________________________
1547 : void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
1548 : {
1549 0 : if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
1550 0 : if(fLowGain)
1551 0 : celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
1552 : else
1553 0 : celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
1554 : }
1555 0 : }
1556 :
1557 : ///
1558 : /// Recalibrate time of cell from SM number considering the L1 phase shift
1559 : ///
1560 : /// \param iSM: supermodule number
1561 : /// \param bc: bunch crossing number returned by esdevent->GetBunchCrossNumber()
1562 : /// \param celltime: cell time to be returned calibrated
1563 : ///
1564 : //_______________________________________________________________________________________________________
1565 : void AliEMCALRecoUtils::RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t & celltime) const
1566 : {
1567 0 : if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
1568 : {
1569 0 : bc=bc%4;
1570 :
1571 : Float_t offsetPerSM=0.;
1572 0 : Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM);
1573 0 : Int_t l1Phase=l1PhaseShift & 3; //bit operation
1574 :
1575 0 : if(bc >= l1Phase)
1576 0 : offsetPerSM = (bc - l1Phase)*25;
1577 : else
1578 0 : offsetPerSM = (bc - l1Phase + 4)*25;
1579 :
1580 0 : Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
1581 0 : l1shiftOffset*=25;
1582 :
1583 0 : celltime -= offsetPerSM*1.e-9;
1584 0 : celltime -= l1shiftOffset*1.e-9;
1585 0 : }
1586 0 : }
1587 :
1588 : ///
1589 : /// For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1590 : ///
1591 : /// \param geom: EMCal geometry pointer
1592 : /// \param cells: list of EMCal cells with signal
1593 : /// \param cluster: EMCal cluster subject to position recalculation
1594 : ///
1595 : //______________________________________________________________________________
1596 : void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1597 : AliVCaloCells* cells,
1598 : AliVCluster* clu)
1599 : {
1600 0 : if (!clu)
1601 : {
1602 0 : AliInfo("Cluster pointer null!");
1603 0 : return;
1604 : }
1605 :
1606 0 : if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1607 0 : else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1608 0 : else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1609 0 : }
1610 :
1611 : ///
1612 : /// For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1613 : /// The algorithm is a copy of what is done in AliEMCALRecPoint.
1614 : ///
1615 : /// \param geom: EMCal geometry pointer
1616 : /// \param cells: list of EMCal cells with signal
1617 : /// \param cluster: EMCal cluster subject to position recalculation
1618 : ///
1619 : //_____________________________________________________________________________________________
1620 : void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom,
1621 : AliVCaloCells* cells,
1622 : AliVCluster* clu)
1623 : {
1624 : Double_t eCell = 0.;
1625 : Float_t fraction = 1.;
1626 : Float_t recalFactor = 1.;
1627 :
1628 0 : Int_t absId = -1;
1629 0 : Int_t iTower = -1, iIphi = -1, iIeta = -1;
1630 0 : Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1631 : Float_t weight = 0., totalWeight=0.;
1632 0 : Float_t newPos[3] = {-1.,-1.,-1.};
1633 0 : Double_t pLocal[3], pGlobal[3];
1634 0 : Bool_t shared = kFALSE;
1635 :
1636 0 : Float_t clEnergy = clu->E(); //Energy already recalibrated previously
1637 0 : if (clEnergy <= 0) return;
1638 :
1639 0 : GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1640 0 : Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1641 :
1642 : //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1643 :
1644 0 : for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1645 : {
1646 0 : absId = clu->GetCellAbsId(iDig);
1647 0 : fraction = clu->GetCellAmplitudeFraction(iDig);
1648 0 : if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1649 :
1650 0 : if (!fCellsRecalibrated)
1651 : {
1652 0 : geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1653 0 : geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1654 0 : if (IsRecalibrationOn()) {
1655 0 : recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1656 0 : }
1657 : }
1658 :
1659 0 : eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1660 :
1661 0 : weight = GetCellWeight(eCell,clEnergy);
1662 0 : totalWeight += weight;
1663 :
1664 0 : geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1665 : //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1666 0 : geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1667 : //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1668 :
1669 0 : for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1670 : }// cell loop
1671 :
1672 0 : if (totalWeight>0)
1673 : {
1674 0 : for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1675 0 : }
1676 :
1677 : //Float_t pos[]={0,0,0};
1678 : //clu->GetPosition(pos);
1679 : //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1680 : //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1681 :
1682 0 : if (iSupModMax > 1) { //sector 1
1683 0 : newPos[0] +=fMisalTransShift[3];//-=3.093;
1684 0 : newPos[1] +=fMisalTransShift[4];//+=6.82;
1685 0 : newPos[2] +=fMisalTransShift[5];//+=1.635;
1686 : //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1687 0 : } else { //sector 0
1688 0 : newPos[0] +=fMisalTransShift[0];//+=1.134;
1689 0 : newPos[1] +=fMisalTransShift[1];//+=8.2;
1690 0 : newPos[2] +=fMisalTransShift[2];//+=1.197;
1691 : //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1692 : }
1693 : //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1694 :
1695 0 : clu->SetPosition(newPos);
1696 0 : }
1697 :
1698 : ///
1699 : /// For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1700 : /// The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position.
1701 : ///
1702 : /// \param geom: EMCal geometry pointer
1703 : /// \param cells: list of EMCal cells with signal
1704 : /// \param cluster: EMCal cluster subject to position recalculation
1705 : ///
1706 : //____________________________________________________________________________________________
1707 : void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom,
1708 : AliVCaloCells* cells,
1709 : AliVCluster* clu)
1710 : {
1711 : Double_t eCell = 1.;
1712 : Float_t fraction = 1.;
1713 : Float_t recalFactor = 1.;
1714 :
1715 0 : Int_t absId = -1;
1716 0 : Int_t iTower = -1;
1717 0 : Int_t iIphi = -1, iIeta = -1;
1718 0 : Int_t iSupMod = -1, iSupModMax = -1;
1719 0 : Int_t iphi = -1, ieta =-1;
1720 0 : Bool_t shared = kFALSE;
1721 :
1722 0 : Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
1723 :
1724 0 : if (clEnergy <= 0)
1725 0 : return;
1726 0 : GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1727 0 : Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1728 :
1729 : Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
1730 : Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1731 : Int_t startingSM = -1;
1732 :
1733 0 : for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1734 : {
1735 0 : absId = clu->GetCellAbsId(iDig);
1736 0 : fraction = clu->GetCellAmplitudeFraction(iDig);
1737 0 : if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1738 :
1739 0 : if (iDig==0) startingSM = iSupMod;
1740 0 : else if (iSupMod != startingSM) areInSameSM = kFALSE;
1741 :
1742 0 : eCell = cells->GetCellAmplitude(absId);
1743 :
1744 0 : geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1745 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1746 :
1747 0 : if (!fCellsRecalibrated)
1748 : {
1749 0 : if (IsRecalibrationOn()) {
1750 0 : recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1751 0 : }
1752 : }
1753 :
1754 0 : eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1755 :
1756 0 : weight = GetCellWeight(eCell,clEnergy);
1757 0 : if (weight < 0) weight = 0;
1758 0 : totalWeight += weight;
1759 0 : weightedCol += ieta*weight;
1760 0 : weightedRow += iphi*weight;
1761 :
1762 : //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1763 : }// cell loop
1764 :
1765 0 : Float_t xyzNew[]={-1.,-1.,-1.};
1766 0 : if (areInSameSM == kTRUE)
1767 : {
1768 : //printf("In Same SM\n");
1769 0 : weightedCol = weightedCol/totalWeight;
1770 0 : weightedRow = weightedRow/totalWeight;
1771 0 : geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1772 0 : }
1773 : else
1774 : {
1775 : //printf("In Different SM\n");
1776 0 : geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1777 : }
1778 :
1779 0 : clu->SetPosition(xyzNew);
1780 0 : }
1781 :
1782 : ///
1783 : /// Re-evaluate distance to bad channel with updated bad map
1784 : ///
1785 : /// \param geom: EMCal geometry pointer
1786 : /// \param cells: list of EMCal cells with signal
1787 : /// \param cluster: EMCal cluster subject to distance to bad channel recalculation
1788 : ///
1789 : //___________________________________________________________________________________________
1790 : void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
1791 : AliVCaloCells* cells,
1792 : AliVCluster * cluster)
1793 : {
1794 0 : if (!fRecalDistToBadChannels) return;
1795 :
1796 0 : if (!cluster)
1797 : {
1798 0 : AliInfo("Cluster pointer null!");
1799 0 : return;
1800 : }
1801 :
1802 : // Get channels map of the supermodule where the cluster is.
1803 0 : Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
1804 0 : Bool_t shared = kFALSE;
1805 0 : GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1806 0 : TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1807 :
1808 : Int_t dRrow, dRcol;
1809 : Float_t minDist = 10000.;
1810 : Float_t dist = 0.;
1811 :
1812 : // Loop on tower status map
1813 0 : for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1814 : {
1815 0 : for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1816 : {
1817 : // Check if tower is bad.
1818 0 : if (hMap->GetBinContent(icol,irow)==0) continue;
1819 : //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
1820 : // iSupMod,icol, irow, icolM,irowM);
1821 :
1822 0 : dRrow=TMath::Abs(irowM-irow);
1823 0 : dRcol=TMath::Abs(icolM-icol);
1824 0 : dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1825 0 : if (dist < minDist)
1826 : {
1827 : //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1828 : minDist = dist;
1829 0 : }
1830 : }
1831 : }
1832 :
1833 : // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
1834 0 : if (shared)
1835 : {
1836 : TH2D* hMap2 = 0;
1837 : Int_t iSupMod2 = -1;
1838 :
1839 : // The only possible combinations are (0,1), (2,3) ... (8,9)
1840 0 : if (iSupMod%2) iSupMod2 = iSupMod-1;
1841 0 : else iSupMod2 = iSupMod+1;
1842 0 : hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
1843 :
1844 : // Loop on tower status map of second super module
1845 0 : for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1846 : {
1847 0 : for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1848 : {
1849 : // Check if tower is bad.
1850 0 : if (hMap2->GetBinContent(icol,irow)==0) continue;
1851 :
1852 : //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",
1853 : // iSupMod2,icol, irow,iSupMod,icolM,irowM);
1854 0 : dRrow=TMath::Abs(irow-irowM);
1855 :
1856 0 : if (iSupMod%2)
1857 0 : dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1858 : else
1859 0 : dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1860 :
1861 0 : dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1862 0 : if (dist < minDist) minDist = dist;
1863 : }
1864 : }
1865 0 : }// shared cluster in 2 SuperModules
1866 :
1867 0 : AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1868 0 : cluster->SetDistanceToBadChannel(minDist);
1869 0 : }
1870 :
1871 : ///
1872 : /// Re-evaluate identification parameters with bayesian
1873 : ///
1874 : /// \param cluster: EMCal cluster subject to PID recalculation
1875 : ///
1876 : ///
1877 : //__________________________________________________________________
1878 : void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1879 : {
1880 0 : if (!cluster)
1881 : {
1882 0 : AliInfo("Cluster pointer null!");
1883 0 : return;
1884 : }
1885 :
1886 0 : if (cluster->GetM02() != 0)
1887 0 : fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1888 :
1889 0 : Float_t pidlist[AliPID::kSPECIESCN+1];
1890 0 : for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1891 :
1892 0 : cluster->SetPID(pidlist);
1893 0 : }
1894 :
1895 : ///
1896 : /// Calculates different types of shower shape parameters, dispersion, shower shape eigenvalues and other.
1897 : ///
1898 : /// \param geom: EMCal geometry pointer
1899 : /// \param cells: list of EMCal cells with signal
1900 : /// \param cluster: EMCal cluster subject to shower shape recalculation
1901 : /// \param cellEcut: minimum cell energy to be considered in the shower shape recalculation
1902 : /// \param cellTimeCut: time window of cells to be considered in shower recalculation
1903 : /// \param bc: event bunch crossing number
1904 : /// \param enAfterCuts: cluster energy when applying the cell cuts cellEcut and cellTime cut
1905 : /// \param l0: main shower shape eigen value
1906 : /// \param l1: second eigenvalue of shower shape
1907 : /// \param disp: dispersion
1908 : /// \param dEta: dispersion in eta (cols) direction
1909 : /// \param dPhi: disperion in phi (rows) direction
1910 : /// \param sEta: shower shape in eta (cols) direction
1911 : /// \param sPhi: shower shape in phi (rows) direction
1912 : /// \param sEtaPhi: shower shape on phi / eta directions term
1913 : ///
1914 : ///
1915 : //___________________________________________________________________________________________________________________
1916 : void AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry * geom,
1917 : AliVCaloCells* cells, AliVCluster * cluster,
1918 : Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
1919 : Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
1920 : Float_t & disp, Float_t & dEta, Float_t & dPhi,
1921 : Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
1922 : {
1923 0 : if (!cluster)
1924 : {
1925 0 : AliInfo("Cluster pointer null!");
1926 0 : return;
1927 : }
1928 :
1929 : Double_t eCell = 0.;
1930 0 : Double_t tCell = 0.;
1931 : Float_t fraction = 1.;
1932 : Float_t recalFactor = 1.;
1933 : Bool_t isLowGain = kFALSE;
1934 :
1935 0 : Int_t iSupMod = -1;
1936 0 : Int_t iTower = -1;
1937 0 : Int_t iIphi = -1;
1938 0 : Int_t iIeta = -1;
1939 0 : Int_t iphi = -1;
1940 0 : Int_t ieta = -1;
1941 : Double_t etai = -1.;
1942 : Double_t phii = -1.;
1943 :
1944 : Int_t nstat = 0 ;
1945 : Float_t wtot = 0.;
1946 : Double_t w = 0.;
1947 : Double_t etaMean = 0.;
1948 : Double_t phiMean = 0.;
1949 :
1950 : // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
1951 : // or the non linearity correction was applied
1952 : // and to check if the cluster is between 2 SM in eta
1953 : Int_t iSM0 = -1;
1954 : Bool_t shared = kFALSE;
1955 : Float_t energy = 0;
1956 0 : enAfterCuts = 0;
1957 :
1958 0 : for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1959 : {
1960 : // Get from the absid the supermodule, tower and eta/phi numbers
1961 :
1962 0 : Int_t absId = cluster->GetCellAbsId(iDigit);
1963 :
1964 0 : geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1965 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1966 :
1967 : // Check if there are cells of different SM
1968 0 : if (iDigit == 0 ) iSM0 = iSupMod;
1969 0 : else if (iSupMod!= iSM0) shared = kTRUE;
1970 :
1971 : // Get the cell energy, if recalibration is on, apply factors
1972 0 : fraction = cluster->GetCellAmplitudeFraction(iDigit);
1973 0 : if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1974 :
1975 0 : if (IsRecalibrationOn())
1976 0 : recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1977 :
1978 :
1979 0 : eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1980 0 : tCell = cells->GetCellTime (absId);
1981 0 : isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
1982 :
1983 0 : RecalibrateCellTime(absId, bc, tCell,isLowGain);
1984 0 : tCell*=1e9;
1985 0 : tCell-=fConstantTimeShift;
1986 :
1987 0 : if(eCell > 0.05) // at least a minimum 50 MeV cut
1988 0 : energy += eCell;
1989 :
1990 0 : if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
1991 0 : enAfterCuts += eCell;
1992 : } // cell loop
1993 :
1994 : // Loop on cells to calculate weights and shower shape terms parameters
1995 0 : for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1996 : {
1997 : // Get from the absid the supermodule, tower and eta/phi numbers
1998 0 : Int_t absId = cluster->GetCellAbsId(iDigit);
1999 :
2000 0 : geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2001 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2002 :
2003 : //Get the cell energy, if recalibration is on, apply factors
2004 0 : fraction = cluster->GetCellAmplitudeFraction(iDigit);
2005 0 : if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2006 :
2007 0 : if (!fCellsRecalibrated && IsRecalibrationOn())
2008 0 : recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2009 :
2010 0 : eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2011 0 : tCell = cells->GetCellTime (absId);
2012 0 : isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2013 :
2014 0 : RecalibrateCellTime(absId, bc, tCell,isLowGain);
2015 0 : tCell*=1e9;
2016 0 : tCell-=fConstantTimeShift;
2017 :
2018 : // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2019 : // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2020 0 : if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2021 :
2022 0 : if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2023 : {
2024 0 : w = GetCellWeight(eCell, energy);
2025 :
2026 0 : etai=(Double_t)ieta;
2027 0 : phii=(Double_t)iphi;
2028 :
2029 0 : if (w > 0.0)
2030 : {
2031 0 : wtot += w ;
2032 0 : nstat++;
2033 :
2034 : // Shower shape
2035 0 : sEta += w * etai * etai ;
2036 0 : etaMean += w * etai ;
2037 0 : sPhi += w * phii * phii ;
2038 0 : phiMean += w * phii ;
2039 0 : sEtaPhi += w * etai * phii ;
2040 0 : }
2041 : }
2042 0 : else if(eCell > 0.05)
2043 0 : AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2044 : } // cell loop
2045 :
2046 : // Normalize to the weight
2047 0 : if (wtot > 0)
2048 : {
2049 0 : etaMean /= wtot ;
2050 0 : phiMean /= wtot ;
2051 0 : }
2052 : else
2053 0 : AliDebug(2,Form("Wrong weight %f\n", wtot));
2054 :
2055 : // Loop on cells to calculate dispersion
2056 0 : for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2057 : {
2058 : // Get from the absid the supermodule, tower and eta/phi numbers
2059 0 : Int_t absId = cluster->GetCellAbsId(iDigit);
2060 :
2061 0 : geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2062 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2063 :
2064 : //Get the cell energy, if recalibration is on, apply factors
2065 0 : fraction = cluster->GetCellAmplitudeFraction(iDigit);
2066 0 : if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2067 0 : if (IsRecalibrationOn())
2068 0 : recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2069 :
2070 0 : eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2071 0 : tCell = cells->GetCellTime (absId);
2072 0 : isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2073 :
2074 0 : RecalibrateCellTime(absId, bc, tCell,isLowGain);
2075 0 : tCell*=1e9;
2076 0 : tCell-=fConstantTimeShift;
2077 :
2078 : // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2079 : // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2080 0 : if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2081 :
2082 0 : if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2083 : {
2084 0 : w = GetCellWeight(eCell,cluster->E());
2085 :
2086 0 : etai=(Double_t)ieta;
2087 0 : phii=(Double_t)iphi;
2088 0 : if (w > 0.0)
2089 : {
2090 0 : disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2091 0 : dEta += w * (etai-etaMean)*(etai-etaMean) ;
2092 0 : dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2093 0 : }
2094 : }
2095 : else
2096 0 : AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2097 : }// cell loop
2098 :
2099 : // Normalize to the weigth and set shower shape parameters
2100 0 : if (wtot > 0 && nstat > 1)
2101 : {
2102 0 : disp /= wtot ;
2103 0 : dEta /= wtot ;
2104 0 : dPhi /= wtot ;
2105 0 : sEta /= wtot ;
2106 0 : sPhi /= wtot ;
2107 0 : sEtaPhi /= wtot ;
2108 :
2109 0 : sEta -= etaMean * etaMean ;
2110 0 : sPhi -= phiMean * phiMean ;
2111 0 : sEtaPhi -= etaMean * phiMean ;
2112 :
2113 0 : l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2114 0 : l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2115 0 : }
2116 : else
2117 : {
2118 0 : l0 = 0. ;
2119 0 : l1 = 0. ;
2120 0 : dEta = 0. ; dPhi = 0. ; disp = 0. ;
2121 0 : sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2122 : }
2123 0 : }
2124 :
2125 : ///
2126 : /// Calculates different types of shower shape parameters, dispersion, shower shape eigenvalues and other.
2127 : /// Call to AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts
2128 : /// with default cell cuts (50 MeV minimum cell energy and not cut on time)
2129 : ///
2130 : /// \param geom: EMCal geometry pointer
2131 : /// \param cells: list of EMCal cells with signal
2132 : /// \param cluster: EMCal cluster subject to shower shape recalculation
2133 : /// \param l0: main shower shape eigen value
2134 : /// \param l1: second eigenvalue of shower shape
2135 : /// \param disp: dispersion
2136 : /// \param dEta: dispersion in eta (cols) direction
2137 : /// \param dPhi: disperion in phi (rows) direction
2138 : /// \param sEta: shower shape in eta (cols) direction
2139 : /// \param sPhi: shower shape in phi (rows) direction
2140 : /// \param sEtaPhi: shower shape on phi / eta directions term
2141 : ///
2142 : ///
2143 : //___________________________________________________________________________________________________________________
2144 : void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
2145 : AliVCaloCells* cells, AliVCluster * cluster,
2146 : Float_t & l0, Float_t & l1,
2147 : Float_t & disp, Float_t & dEta, Float_t & dPhi,
2148 : Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2149 : {
2150 0 : Float_t newEnergy = 0;
2151 : Float_t cellEmin = 0.05; // 50 MeV
2152 : Float_t cellTimeWindow = 1000; // open cut
2153 : Int_t bc = 0;
2154 0 : AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(geom, cells, cluster,
2155 : cellEmin, cellTimeWindow, bc,
2156 : newEnergy, l0, l1, disp,
2157 : dEta, dPhi, sEta, sPhi, sEtaPhi);
2158 0 : }
2159 :
2160 : ///
2161 : /// Calculates Dispersion and main axis and puts them into the cluster
2162 : /// Call to method RecalculateClusterShowerShapeParameters
2163 : ///
2164 : /// \param geom: EMCal geometry pointer
2165 : /// \param cells: list of EMCal cells with signal
2166 : /// \param cluster: EMCal cluster subject to shower shape recalculation
2167 : ///
2168 : //____________________________________________________________________________________________
2169 : void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
2170 : AliVCaloCells* cells,
2171 : AliVCluster * cluster)
2172 : {
2173 0 : Float_t l0 = 0., l1 = 0.;
2174 0 : Float_t disp = 0., dEta = 0., dPhi = 0.;
2175 0 : Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2176 :
2177 0 : AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
2178 : dEta, dPhi, sEta, sPhi, sEtaPhi);
2179 :
2180 0 : cluster->SetM02(l0);
2181 0 : cluster->SetM20(l1);
2182 0 : if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2183 :
2184 0 : }
2185 :
2186 : ///
2187 : /// Calculates Dispersion and main axis and puts them into the cluster.
2188 : /// Possibility to restrict the cell Energy and time window in the calculations
2189 : ///
2190 : /// \param geom: EMCal geometry pointer
2191 : /// \param cells: list of EMCal cells with signal
2192 : /// \param cluster: EMCal cluster subject to shower shape recalculation
2193 : /// \param cellEcut: minimum cell energy to be considered in the shower shape recalculation
2194 : /// \param cellTimeCut: time window of cells to be considered in shower recalculation
2195 : /// \param bc: event bunch crossing number
2196 : /// \param enAfterCuts: cluster energy when applying the cell cuts cellEcut and cellTime cut
2197 : ///
2198 : //____________________________________________________________________________________________
2199 : void AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry * geom,
2200 : AliVCaloCells* cells, AliVCluster * cluster,
2201 : Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2202 : Float_t & enAfterCuts)
2203 : {
2204 0 : Float_t l0 = 0., l1 = 0.;
2205 0 : Float_t disp = 0., dEta = 0., dPhi = 0.;
2206 0 : Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2207 :
2208 0 : AliEMCALRecoUtils::RecalculateClusterShowerShapeParametersWithCellCuts(geom, cells, cluster,
2209 : cellEcut, cellTimeCut, bc,
2210 : enAfterCuts, l0, l1, disp,
2211 : dEta, dPhi, sEta, sPhi, sEtaPhi);
2212 :
2213 0 : cluster->SetM02(l0);
2214 0 : cluster->SetM20(l1);
2215 0 : if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2216 :
2217 0 : }
2218 :
2219 : ///
2220 : /// Find the candidate cluster-track matchs.
2221 : ///
2222 : /// This function should be called before the cluster loop.
2223 : /// Before call this function, please recalculate the cluster positions.
2224 : /// Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR.
2225 : /// Store matched cluster indexes and residuals.
2226 : ///
2227 : /// \param event: event pointer
2228 : /// \param clusterArr: list of clusters
2229 : /// \param geom: AliEMCALGeometry pointer
2230 : ///
2231 : //____________________________________________________________________________
2232 : void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
2233 : TObjArray * clusterArr,
2234 : const AliEMCALGeometry *geom)
2235 : {
2236 0 : fMatchedTrackIndex ->Reset();
2237 0 : fMatchedClusterIndex->Reset();
2238 0 : fResidualPhi->Reset();
2239 0 : fResidualEta->Reset();
2240 :
2241 0 : fMatchedTrackIndex ->Set(1000);
2242 0 : fMatchedClusterIndex->Set(1000);
2243 0 : fResidualPhi->Set(1000);
2244 0 : fResidualEta->Set(1000);
2245 :
2246 0 : AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
2247 0 : AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
2248 :
2249 : // Init the magnetic field if not already on
2250 0 : if (!TGeoGlobalMagField::Instance()->GetField())
2251 : {
2252 0 : if (!event->InitMagneticField())
2253 : {
2254 0 : AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
2255 0 : }
2256 : } // Init mag field
2257 :
2258 0 : if (esdevent)
2259 : {
2260 0 : UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
2261 0 : UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
2262 0 : Bool_t desc1 = (mask1 >> 3) & 0x1;
2263 0 : Bool_t desc2 = (mask2 >> 3) & 0x1;
2264 0 : if (desc1==0 || desc2==0) {
2265 : // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
2266 : // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
2267 : // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
2268 0 : fITSTrackSA=kTRUE;
2269 0 : }
2270 0 : }
2271 :
2272 : TObjArray *clusterArray = 0x0;
2273 0 : if (!clusterArr)
2274 : {
2275 0 : clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
2276 0 : for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2277 : {
2278 0 : AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2279 0 : if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
2280 0 : continue;
2281 0 : clusterArray->AddAt(cluster,icl);
2282 0 : }
2283 0 : }
2284 :
2285 : Int_t matched=0;
2286 0 : Double_t cv[21];
2287 0 : for (Int_t i=0; i<21;i++) cv[i]=0;
2288 0 : for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
2289 : {
2290 : AliExternalTrackParam *trackParam = 0;
2291 :
2292 : // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
2293 : AliESDtrack *esdTrack = 0;
2294 : AliAODTrack *aodTrack = 0;
2295 0 : if (esdevent)
2296 : {
2297 0 : esdTrack = esdevent->GetTrack(itr);
2298 0 : if (!esdTrack) continue;
2299 0 : if (!IsAccepted(esdTrack)) continue;
2300 0 : if (esdTrack->Pt()<fCutMinTrackPt) continue;
2301 :
2302 0 : if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
2303 :
2304 : // Save some time and memory in case of no DCal present
2305 0 : if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2306 : {
2307 0 : Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
2308 0 : if ( phi <= 10 || phi >= 250 ) continue;
2309 0 : }
2310 :
2311 0 : if (!fITSTrackSA)
2312 0 : trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
2313 : else
2314 0 : trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
2315 : }
2316 :
2317 : // If the input event is AOD, the starting point for extrapolation is at vertex
2318 : // AOD tracks are selected according to its filterbit.
2319 0 : else if (aodevent)
2320 : {
2321 0 : aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
2322 :
2323 0 : if (!aodTrack) AliFatal("Not a standard AOD");
2324 :
2325 0 : if (!aodTrack) continue;
2326 :
2327 0 : if (fAODTPCOnlyTracks)
2328 : { // Match with TPC only tracks, default from May 2013, before filter bit 32
2329 : //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
2330 0 : if (!aodTrack->IsTPCConstrained()) continue ;
2331 : }
2332 0 : else if (fAODHybridTracks)
2333 : { // Match with hybrid tracks
2334 : //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
2335 0 : if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
2336 : } else
2337 : { // Match with tracks on a mask
2338 : //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
2339 0 : if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
2340 : }
2341 :
2342 0 : if (aodTrack->Pt() < fCutMinTrackPt) continue;
2343 :
2344 0 : if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
2345 :
2346 : // Save some time and memory in case of no DCal present
2347 0 : if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2348 : {
2349 0 : Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
2350 0 : if ( phi <= 10 || phi >= 250 ) continue;
2351 0 : }
2352 :
2353 0 : Double_t pos[3],mom[3];
2354 0 : aodTrack->GetXYZ(pos);
2355 0 : aodTrack->GetPxPyPz(mom);
2356 0 : AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
2357 : itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
2358 :
2359 0 : trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
2360 0 : }
2361 :
2362 : //Return if the input data is not "AOD" or "ESD"
2363 : else
2364 : {
2365 0 : printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
2366 0 : if (clusterArray)
2367 : {
2368 0 : clusterArray->Clear();
2369 0 : delete clusterArray;
2370 : }
2371 0 : return;
2372 : }
2373 :
2374 0 : if (!trackParam) continue;
2375 :
2376 : // Extrapolate the track to EMCal surface
2377 0 : AliExternalTrackParam emcalParam(*trackParam);
2378 0 : Float_t eta, phi, pt;
2379 0 : if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2380 : {
2381 0 : if (aodevent && trackParam) delete trackParam;
2382 0 : if (fITSTrackSA && trackParam) delete trackParam;
2383 0 : continue;
2384 : }
2385 :
2386 0 : if ( TMath::Abs(eta) > 0.75 )
2387 : {
2388 0 : if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2389 0 : continue;
2390 : }
2391 :
2392 : // Save some time and memory in case of no DCal present
2393 0 : if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2394 0 : ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2395 : {
2396 0 : if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2397 0 : continue;
2398 : }
2399 :
2400 : //Find matched clusters
2401 : Int_t index = -1;
2402 0 : Float_t dEta = -999, dPhi = -999;
2403 0 : if (!clusterArr)
2404 0 : index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
2405 : else
2406 0 : index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
2407 :
2408 :
2409 0 : if (index>-1)
2410 : {
2411 0 : fMatchedTrackIndex ->AddAt(itr,matched);
2412 0 : fMatchedClusterIndex ->AddAt(index,matched);
2413 0 : fResidualEta ->AddAt(dEta,matched);
2414 0 : fResidualPhi ->AddAt(dPhi,matched);
2415 0 : matched++;
2416 0 : }
2417 :
2418 0 : if (aodevent && trackParam) delete trackParam;
2419 0 : if (fITSTrackSA && trackParam) delete trackParam;
2420 0 : }//track loop
2421 :
2422 0 : if (clusterArray)
2423 : {
2424 0 : clusterArray->Clear();
2425 0 : delete clusterArray;
2426 : }
2427 :
2428 0 : AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2429 :
2430 0 : fMatchedTrackIndex ->Set(matched);
2431 0 : fMatchedClusterIndex ->Set(matched);
2432 0 : fResidualPhi ->Set(matched);
2433 0 : fResidualEta ->Set(matched);
2434 0 : }
2435 :
2436 : ///
2437 : /// Find matched cluster in event. See Find MatchedClusterInClusterArr().
2438 : ///
2439 : /// \param track: track pointer
2440 : /// \param event: event pointer
2441 : /// \param geom: AliEMCALGeometry pointer
2442 : /// \param dEta: found track-cluster match residual in eta direction
2443 : /// \param dPhi: found track-cluster match residual in phi direction
2444 : ///
2445 : /// \return the index of matched cluster to input track, returns -1 if no match is found
2446 : ///
2447 : //________________________________________________________________________________
2448 : Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
2449 : const AliVEvent *event,
2450 : const AliEMCALGeometry *geom,
2451 : Float_t &dEta, Float_t &dPhi)
2452 : {
2453 : Int_t index = -1;
2454 :
2455 0 : if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
2456 :
2457 : // Save some time and memory in case of no DCal present
2458 0 : if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2459 : {
2460 0 : Double_t phiV = track->Phi()*TMath::RadToDeg();
2461 0 : if ( phiV <= 10 || phiV >= 250 ) return index;
2462 0 : }
2463 :
2464 : AliExternalTrackParam *trackParam = 0;
2465 0 : if (!fITSTrackSA)
2466 0 : trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
2467 : else
2468 0 : trackParam = new AliExternalTrackParam(*track);
2469 :
2470 0 : if (!trackParam) return index;
2471 0 : AliExternalTrackParam emcalParam(*trackParam);
2472 0 : Float_t eta, phi, pt;
2473 :
2474 0 : if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2475 : {
2476 0 : if (fITSTrackSA) delete trackParam;
2477 0 : return index;
2478 : }
2479 :
2480 0 : if ( TMath::Abs(eta) > 0.75 )
2481 : {
2482 0 : if (fITSTrackSA) delete trackParam;
2483 0 : return index;
2484 : }
2485 :
2486 : // Save some time and memory in case of no DCal present
2487 0 : if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2488 0 : ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2489 : {
2490 0 : if (fITSTrackSA) delete trackParam;
2491 0 : return index;
2492 : }
2493 :
2494 0 : TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2495 :
2496 0 : for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2497 : {
2498 0 : AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2499 0 : if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2500 0 : clusterArr->AddAt(cluster,icl);
2501 0 : }
2502 :
2503 0 : index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2504 0 : clusterArr->Clear();
2505 0 : delete clusterArr;
2506 0 : if (fITSTrackSA) delete trackParam;
2507 :
2508 : return index;
2509 0 : }
2510 :
2511 : ///
2512 : /// Find matched cluster in input array of clusters.
2513 : ///
2514 : /// \param emcalParam: emcal track parameters container?
2515 : /// \param trkParam: track parameters container?
2516 : /// \param clusterArr: input array of clusters
2517 : /// \param dEta: found track-cluster match residual in eta direction
2518 : /// \param dPhi: found track-cluster match residual in phi direction
2519 : ///
2520 : /// \return the index of matched cluster to input track.
2521 : //_______________________________________________________________________________________________
2522 : Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2523 : AliExternalTrackParam *trkParam,
2524 : const TObjArray * clusterArr,
2525 : Float_t &dEta, Float_t &dPhi)
2526 : {
2527 0 : dEta=-999, dPhi=-999;
2528 0 : Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2529 : Int_t index = -1;
2530 0 : Float_t tmpEta=-999, tmpPhi=-999;
2531 :
2532 0 : Double_t exPos[3] = {0.,0.,0.};
2533 0 : if (!emcalParam->GetXYZ(exPos)) return index;
2534 :
2535 0 : Float_t clsPos[3] = {0.,0.,0.};
2536 0 : for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2537 : {
2538 0 : AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2539 :
2540 0 : if (!cluster || !cluster->IsEMCAL()) continue;
2541 :
2542 0 : cluster->GetPosition(clsPos);
2543 :
2544 0 : 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));
2545 0 : if (dR > fClusterWindow) continue;
2546 :
2547 0 : AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
2548 :
2549 0 : if (!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2550 :
2551 0 : if (fCutEtaPhiSum)
2552 : {
2553 0 : Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2554 0 : if (tmpR<dRMax)
2555 : {
2556 : dRMax=tmpR;
2557 0 : dEtaMax=tmpEta;
2558 0 : dPhiMax=tmpPhi;
2559 : index=icl;
2560 0 : }
2561 0 : }
2562 0 : else if (fCutEtaPhiSeparate)
2563 : {
2564 0 : if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
2565 : {
2566 0 : dEtaMax = tmpEta;
2567 0 : dPhiMax = tmpPhi;
2568 : index=icl;
2569 0 : }
2570 : }
2571 : else
2572 : {
2573 0 : printf("Error: please specify your cut criteria\n");
2574 0 : printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2575 0 : printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2576 0 : return index;
2577 : }
2578 0 : }
2579 :
2580 0 : dEta=dEtaMax;
2581 0 : dPhi=dPhiMax;
2582 :
2583 0 : return index;
2584 0 : }
2585 :
2586 : ///
2587 : /// Extrapolate track to EMCAL surface.
2588 : ///
2589 : /// \param track: pointer to track
2590 : /// \param emcalR: distance
2591 : /// \param mass: mass hypothesis
2592 : /// \param step: ...
2593 : /// \param minpt: minimum track pT for matching
2594 : /// \param useMassForTracking: switch to use the mass hypothesis
2595 : ///
2596 : /// \return bool true if track could be extrapolated
2597 : ///
2598 : //------------------------------------------------------------------------------------
2599 : Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliVTrack *track,
2600 : Double_t emcalR, Double_t mass,
2601 : Double_t step, Double_t minpt,
2602 : Bool_t useMassForTracking)
2603 : {
2604 142 : track->SetTrackPhiEtaPtOnEMCal(-999, -999, -999);
2605 :
2606 142 : if ( track->Pt() < minpt )
2607 24 : return kFALSE;
2608 :
2609 118 : if ( TMath::Abs(track->Eta()) > 0.9 ) return kFALSE;
2610 :
2611 : // Save some time and memory in case of no DCal present
2612 118 : AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
2613 :
2614 : Int_t nSupMod = AliEMCALGeoParams::fgkEMCALModules;
2615 236 : if ( geom ) nSupMod = geom->GetNumberOfSuperModules();
2616 :
2617 118 : if ( nSupMod < 13 ) // Run1 10 (12, 2 not active but present)
2618 : {
2619 118 : Double_t phi = track->Phi()*TMath::RadToDeg();
2620 138 : if ( phi <= 10 || phi >= 250 ) return kFALSE;
2621 98 : }
2622 :
2623 294 : AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(track);
2624 : AliAODTrack *aodt = 0;
2625 98 : if (!esdt)
2626 : {
2627 0 : aodt = dynamic_cast<AliAODTrack*>(track);
2628 :
2629 0 : if (!aodt)
2630 0 : return kFALSE;
2631 : }
2632 :
2633 : // Select the mass hypothesis
2634 98 : if ( mass < 0 )
2635 : {
2636 : Bool_t onlyTPC = kFALSE;
2637 0 : if ( mass == -99 ) onlyTPC=kTRUE;
2638 :
2639 0 : if (esdt)
2640 : {
2641 0 : if ( useMassForTracking ) mass = esdt->GetMassForTracking();
2642 0 : else mass = esdt->GetMass(onlyTPC);
2643 : }
2644 : else
2645 : {
2646 0 : if ( useMassForTracking ) mass = aodt->GetMassForTracking();
2647 0 : else mass = aodt->M();
2648 : }
2649 0 : }
2650 :
2651 : AliExternalTrackParam *trackParam = 0;
2652 98 : if (esdt)
2653 : {
2654 98 : const AliExternalTrackParam *in = esdt->GetInnerParam();
2655 :
2656 98 : if (!in)
2657 14 : return kFALSE;
2658 :
2659 84 : trackParam = new AliExternalTrackParam(*in);
2660 84 : }
2661 : else
2662 : {
2663 0 : Double_t xyz[3] = {0}, pxpypz[3] = {0}, cv[21] = {0};
2664 0 : aodt->PxPyPz(pxpypz);
2665 0 : aodt->XvYvZv(xyz);
2666 0 : aodt->GetCovarianceXYZPxPyPz(cv);
2667 0 : trackParam = new AliExternalTrackParam(xyz,pxpypz,cv,aodt->Charge());
2668 0 : }
2669 :
2670 84 : if (!trackParam)
2671 0 : return kFALSE;
2672 :
2673 84 : Float_t etaout=-999, phiout=-999, ptout=-999;
2674 84 : Bool_t ret = ExtrapolateTrackToEMCalSurface(trackParam,
2675 : emcalR,
2676 : mass,
2677 : step,
2678 : etaout,
2679 : phiout,
2680 : ptout);
2681 :
2682 168 : delete trackParam;
2683 :
2684 86 : if (!ret) return kFALSE;
2685 :
2686 82 : if ( TMath::Abs(etaout) > 0.75 ) return kFALSE;
2687 :
2688 : // Save some time and memory in case of no DCal present
2689 82 : if ( nSupMod < 13 ) // Run1 10 (12, 2 not active but present)
2690 : {
2691 170 : if ( (phiout < 70*TMath::DegToRad()) || (phiout > 190*TMath::DegToRad()) ) return kFALSE;
2692 : }
2693 :
2694 26 : track->SetTrackPhiEtaPtOnEMCal(phiout, etaout, ptout);
2695 :
2696 26 : return kTRUE;
2697 226 : }
2698 :
2699 :
2700 : ///
2701 : /// Extrapolate track to EMCAL surface.
2702 : ///
2703 : /// \param trkParam: pointer to external track param
2704 : /// \param emcalR: distance
2705 : /// \param mass: mass hypothesis
2706 : /// \param step: ...
2707 : /// \param eta: track extrapolated eta
2708 : /// \param phi: track extrapolated phi
2709 : /// \param pt: track extrapolated pt
2710 : ///
2711 : /// \return bool true if track could be extrapolated
2712 : ///
2713 : //---------------------------------------------------------------------------------------
2714 : Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
2715 : Double_t emcalR,
2716 : Double_t mass,
2717 : Double_t step,
2718 : Float_t &eta,
2719 : Float_t &phi,
2720 : Float_t &pt)
2721 : {
2722 380 : eta = -999, phi = -999, pt = -999;
2723 :
2724 190 : if (!trkParam) return kFALSE;
2725 :
2726 217 : if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
2727 :
2728 163 : Double_t trkPos[3] = {0.,0.,0.};
2729 :
2730 163 : if (!trkParam->GetXYZ(trkPos)) return kFALSE;
2731 :
2732 163 : TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2733 :
2734 326 : eta = trkPosVec.Eta();
2735 326 : phi = trkPosVec.Phi();
2736 326 : pt = trkParam->Pt();
2737 :
2738 163 : if ( phi < 0 )
2739 19 : phi += TMath::TwoPi();
2740 :
2741 : return kTRUE;
2742 516 : }
2743 :
2744 : ///
2745 : /// Return the residual by extrapolating a track param to a global position.
2746 : ///
2747 : /// \param trkParam: pointer to external track param
2748 : /// \param clsPos: cluster position array xyz
2749 : /// \param mass: mass hypothesis
2750 : /// \param step: ...
2751 : /// \param tmpEta: residual eta
2752 : /// \param tmpPhi: residual phi
2753 : ///
2754 : /// \return bool true if track could be extrapolated
2755 : //-----------------------------------------------------------------------------------
2756 : Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2757 : const Float_t *clsPos,
2758 : Double_t mass,
2759 : Double_t step,
2760 : Float_t &tmpEta,
2761 : Float_t &tmpPhi)
2762 : {
2763 112 : tmpEta = -999;
2764 56 : tmpPhi = -999;
2765 :
2766 56 : if (!trkParam) return kFALSE;
2767 :
2768 56 : Double_t trkPos[3] = {0.,0.,0.};
2769 56 : TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2770 :
2771 112 : Float_t phi = vec.Phi();
2772 :
2773 56 : if ( phi < 0 )
2774 0 : phi += TMath::TwoPi();
2775 :
2776 : // Rotate the cluster to the local extrapolation coordinate system
2777 56 : Double_t alpha = ((int)(phi*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2778 :
2779 56 : vec.RotateZ(-alpha);
2780 :
2781 112 : if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1))
2782 2 : return kFALSE;
2783 :
2784 108 : if (!trkParam->GetXYZ(trkPos))
2785 0 : return kFALSE; // Get the extrapolated global position
2786 :
2787 54 : TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2788 54 : TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2789 :
2790 : // Track cluster matching
2791 108 : tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
2792 162 : tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2793 :
2794 : return kTRUE;
2795 166 : }
2796 :
2797 : ///
2798 : /// Return the residual by extrapolating a track param to a cluster.
2799 : ///
2800 : /// \param trkParam: pointer to external track param
2801 : /// \param cluster: pointer to AliVCluster
2802 : /// \param mass: mass hypothesis
2803 : /// \param step: ...
2804 : /// \param tmpEta: residual eta
2805 : /// \param tmpPhi: residual phi
2806 : ///
2807 : /// \return bool true if track could be extrapolated
2808 : ///
2809 : //----------------------------------------------------------------------------------
2810 : Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
2811 : const AliVCluster *cluster,
2812 : Double_t mass,
2813 : Double_t step,
2814 : Float_t &tmpEta,
2815 : Float_t &tmpPhi)
2816 : {
2817 36 : tmpEta = -999;
2818 18 : tmpPhi = -999;
2819 :
2820 18 : if (!cluster || !trkParam)
2821 0 : return kFALSE;
2822 :
2823 18 : Float_t clsPos[3] = {0.,0.,0.};
2824 18 : cluster->GetPosition(clsPos);
2825 :
2826 18 : return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2827 36 : }
2828 :
2829 : ///
2830 : ///
2831 : /// Return the residual by extrapolating a track param to a cluster.
2832 : /// Mass and step hypothesis are set via data members fSteCluster and fMass and
2833 : /// not passed like in method avobe.
2834 : ///
2835 : /// \param trkParam: pointer to external track param
2836 : /// \param cluster: pointer to AliVCluster
2837 : /// \param tmpEta: residual eta
2838 : /// \param tmpPhi: residual phi
2839 : ///
2840 : /// \return bool true if track could be extrapolated
2841 : ///
2842 : //---------------------------------------------------------------------------------
2843 : Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
2844 : const AliVCluster *cluster,
2845 : Float_t &tmpEta,
2846 : Float_t &tmpPhi)
2847 : {
2848 0 : return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
2849 : }
2850 :
2851 : ///
2852 : /// Given a cluster index, get the residuals dEta and dPhi for this cluster
2853 : /// to the closest track.
2854 : /// It works with ESDs and AODs.
2855 : ///
2856 : /// \param clsIndex: cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2857 : /// \param dEta: residual eta
2858 : /// \param dPhi: residual phi
2859 : ///
2860 : //_______________________________________________________________________
2861 : void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex,
2862 : Float_t &dEta, Float_t &dPhi)
2863 : {
2864 0 : if (FindMatchedPosForCluster(clsIndex) >= 999)
2865 : {
2866 0 : AliDebug(2,"No matched tracks found!\n");
2867 0 : dEta=999.;
2868 0 : dPhi=999.;
2869 0 : return;
2870 : }
2871 :
2872 0 : dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2873 0 : dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
2874 0 : }
2875 :
2876 : ///
2877 : /// Given a track index, get the residuals dEta and dPhi for this track
2878 : /// to the closest cluster.
2879 : /// It works with ESDs and AODs.
2880 : ///
2881 : /// \param trkIndex: cluster index as in AliESDEvent::GetTrack(trkIndex)
2882 : /// \param dEta: residual eta
2883 : /// \param dPhi: residual phi
2884 : ///
2885 : //______________________________________________________________________________________________
2886 : void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
2887 : {
2888 0 : if (FindMatchedPosForTrack(trkIndex) >= 999)
2889 : {
2890 0 : AliDebug(2,"No matched cluster found!\n");
2891 0 : dEta=999.;
2892 0 : dPhi=999.;
2893 0 : return;
2894 : }
2895 :
2896 0 : dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2897 0 : dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
2898 0 : }
2899 :
2900 : ///
2901 : /// Given a cluster index , get the index of matched track to this cluster.
2902 : /// It works with ESDs and AODs.
2903 : ///
2904 : /// \param clsIndex: cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2905 : ///
2906 : //__________________________________________________________
2907 : Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2908 : {
2909 0 : if (IsClusterMatched(clsIndex))
2910 0 : return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2911 : else
2912 0 : return -1;
2913 0 : }
2914 :
2915 : ///
2916 : /// Given a track index, get the index of matched cluster to this track.
2917 : /// It works with ESDs and AODs.
2918 : ///
2919 : /// \param trkIndex: cluster index as in AliESDEvent::GetTrack(trkIndex)
2920 : ///
2921 : //__________________________________________________________
2922 : Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
2923 : {
2924 0 : if (IsTrackMatched(trkIndex))
2925 0 : return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
2926 : else
2927 0 : return -1;
2928 0 : }
2929 :
2930 : ///
2931 : /// Given a cluster index, it returns if the cluster has a match.
2932 : ///
2933 : /// \param trkIndex: cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2934 : ///
2935 : /// \return bool true if cluster is matched
2936 : ///
2937 : //______________________________________________________________
2938 : Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
2939 : {
2940 0 : if (FindMatchedPosForCluster(clsIndex) < 999)
2941 0 : return kTRUE;
2942 : else
2943 0 : return kFALSE;
2944 0 : }
2945 :
2946 : ///
2947 : /// Given a track index, it returns if the track has a match.
2948 : ///
2949 : /// \param clsIndex: cluster index as in AliESDEvent::GetTrack(trkIndex)
2950 : ///
2951 : /// \return bool true if cluster is matched
2952 : ///
2953 : //____________________________________________________________
2954 : Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
2955 : {
2956 0 : if (FindMatchedPosForTrack(trkIndex) < 999)
2957 0 : return kTRUE;
2958 : else
2959 0 : return kFALSE;
2960 0 : }
2961 :
2962 : ///
2963 : /// Given a cluster index, it returns the position of the match in the fMatchedClusterIndex array
2964 : ///
2965 : /// \param clsIndex: cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2966 : ///
2967 : /// \return cluster position index
2968 : ///
2969 : //______________________________________________________________________
2970 : UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
2971 : {
2972 0 : Float_t tmpR = fCutR;
2973 : UInt_t pos = 999;
2974 :
2975 0 : for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
2976 : {
2977 0 : if (fMatchedClusterIndex->At(i)==clsIndex)
2978 : {
2979 0 : Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
2980 0 : if (r<tmpR)
2981 : {
2982 : pos=i;
2983 : tmpR=r;
2984 0 : AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2985 : fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
2986 : }
2987 0 : }
2988 : }
2989 0 : return pos;
2990 : }
2991 :
2992 : ///
2993 : /// Given a cluster index, it returns the position of the match in the fMatchedTrackIndex array
2994 : ///
2995 : /// \param trkIndex: cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2996 : ///
2997 : /// \return track position index
2998 : ///
2999 : //____________________________________________________________________
3000 : UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
3001 : {
3002 0 : Float_t tmpR = fCutR;
3003 : UInt_t pos = 999;
3004 :
3005 0 : for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3006 : {
3007 0 : if (fMatchedTrackIndex->At(i)==trkIndex)
3008 : {
3009 0 : Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3010 0 : if (r<tmpR)
3011 : {
3012 : pos=i;
3013 : tmpR=r;
3014 0 : AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3015 : fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3016 : }
3017 0 : }
3018 : }
3019 0 : return pos;
3020 : }
3021 :
3022 : ///
3023 : /// Check if the cluster survives the following quality cuts:
3024 : /// * Cluster pointer is non null and is EMCal (or DCal)
3025 : /// * There is no bad channel cell inside
3026 : /// * The cluster is not too close to a fiducial border
3027 : /// * The cluster is not exotic
3028 : ///
3029 : /// \param cluster: pointer to AliVCluster
3030 : /// \param geom: pointer to AliEMCALGeometry
3031 : /// \param cells: full list of cells
3032 : /// \param bc: event bunch crossing number
3033 : ///
3034 : /// \return true if cluster passes all selection criteria
3035 : ///
3036 : //__________________________________________________________________________
3037 : Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
3038 : const AliEMCALGeometry *geom,
3039 : AliVCaloCells* cells, Int_t bc)
3040 : {
3041 : Bool_t isGood=kTRUE;
3042 :
3043 0 : if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3044 0 : if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3045 0 : if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3046 0 : if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3047 :
3048 0 : return isGood;
3049 0 : }
3050 :
3051 : ///
3052 : /// Given a esd track, return whether the track survive all the cuts.
3053 : ///
3054 : /// The different quality parameter are first.
3055 : /// retrieved from the track. then it is found out what cuts the
3056 : /// track did not survive and finally the cuts are imposed.
3057 : ///
3058 : /// \param esdTrack: pointer to ESD track
3059 : ///
3060 : /// \return true if track passes all selection criteria
3061 : //__________________________________________________________
3062 : Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
3063 : {
3064 0 : UInt_t status = esdTrack->GetStatus();
3065 :
3066 0 : Int_t nClustersITS = esdTrack->GetITSclusters(0);
3067 0 : Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3068 :
3069 : Float_t chi2PerClusterITS = -1;
3070 : Float_t chi2PerClusterTPC = -1;
3071 0 : if (nClustersITS!=0)
3072 0 : chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3073 0 : if (nClustersTPC!=0)
3074 0 : chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3075 :
3076 : //
3077 : // DCA cuts
3078 : // Only to be used for primary
3079 : //
3080 0 : if ( fTrackCutsType == kGlobalCut )
3081 : {
3082 0 : Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3083 : // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3084 :
3085 : //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3086 :
3087 0 : SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3088 0 : }
3089 0 : else if( fTrackCutsType == kGlobalCut2011 )
3090 : {
3091 0 : Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3092 : // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3093 :
3094 : //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3095 :
3096 0 : SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3097 0 : }
3098 :
3099 0 : Float_t b[2];
3100 0 : Float_t bCov[3];
3101 0 : esdTrack->GetImpactParameters(b,bCov);
3102 0 : if (bCov[0]<=0 || bCov[2]<=0)
3103 : {
3104 0 : AliDebug(1, "Estimated b resolution lower or equal zero!");
3105 0 : bCov[0]=0; bCov[2]=0;
3106 0 : }
3107 :
3108 0 : Float_t dcaToVertexXY = b[0];
3109 0 : Float_t dcaToVertexZ = b[1];
3110 : Float_t dcaToVertex = -1;
3111 :
3112 0 : if (fCutDCAToVertex2D)
3113 0 : dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3114 0 : dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3115 : else
3116 0 : dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3117 :
3118 : // cut the track?
3119 :
3120 0 : Bool_t cuts[kNCuts];
3121 0 : for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3122 :
3123 : // track quality cuts
3124 0 : if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3125 0 : cuts[0]=kTRUE;
3126 0 : if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3127 0 : cuts[1]=kTRUE;
3128 0 : if (nClustersTPC<fCutMinNClusterTPC)
3129 0 : cuts[2]=kTRUE;
3130 0 : if (nClustersITS<fCutMinNClusterITS)
3131 0 : cuts[3]=kTRUE;
3132 0 : if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3133 0 : cuts[4]=kTRUE;
3134 0 : if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3135 0 : cuts[5]=kTRUE;
3136 0 : if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3137 0 : cuts[6]=kTRUE;
3138 0 : if (fCutDCAToVertex2D && dcaToVertex > 1)
3139 0 : cuts[7] = kTRUE;
3140 0 : if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3141 0 : cuts[8] = kTRUE;
3142 0 : if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3143 0 : cuts[9] = kTRUE;
3144 :
3145 0 : if (fTrackCutsType == kGlobalCut || fTrackCutsType == kGlobalCut2011)
3146 : {
3147 : //Require at least one SPD point + anything else in ITS
3148 0 : if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3149 0 : cuts[10] = kTRUE;
3150 : }
3151 :
3152 : // ITS
3153 0 : if (fCutRequireITSStandAlone || fCutRequireITSpureSA)
3154 : {
3155 0 : if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3156 : {
3157 : // TPC tracks
3158 0 : cuts[11] = kTRUE;
3159 0 : }
3160 : else
3161 : {
3162 : // ITS standalone tracks
3163 0 : if (fCutRequireITSStandAlone && !fCutRequireITSpureSA)
3164 : {
3165 0 : if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3166 : }
3167 0 : else if (fCutRequireITSpureSA)
3168 : {
3169 0 : if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3170 : }
3171 : }
3172 : }
3173 :
3174 : Bool_t cut=kFALSE;
3175 0 : for (Int_t i=0; i<kNCuts; i++)
3176 0 : if (cuts[i]) { cut = kTRUE ; }
3177 :
3178 : // cut the track
3179 0 : if (cut)
3180 0 : return kFALSE;
3181 : else
3182 0 : return kTRUE;
3183 0 : }
3184 :
3185 : ///
3186 : /// Initialize the track cut criteria.
3187 : /// By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts().
3188 : /// Also, you can customize the cuts using the setters.
3189 : ///
3190 : //_____________________________________
3191 : void AliEMCALRecoUtils::InitTrackCuts()
3192 : {
3193 0 : switch (fTrackCutsType)
3194 : {
3195 : case kTPCOnlyCut:
3196 : {
3197 0 : AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3198 : //TPC
3199 0 : SetMinNClustersTPC(70);
3200 0 : SetMaxChi2PerClusterTPC(4);
3201 0 : SetAcceptKinkDaughters(kFALSE);
3202 0 : SetRequireTPCRefit(kFALSE);
3203 :
3204 : //ITS
3205 0 : SetRequireITSRefit(kFALSE);
3206 0 : SetMaxDCAToVertexZ(3.2);
3207 0 : SetMaxDCAToVertexXY(2.4);
3208 0 : SetDCAToVertex2D(kTRUE);
3209 :
3210 0 : break;
3211 : }
3212 :
3213 : case kGlobalCut:
3214 : {
3215 0 : AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3216 : //TPC
3217 0 : SetMinNClustersTPC(70);
3218 0 : SetMaxChi2PerClusterTPC(4);
3219 0 : SetAcceptKinkDaughters(kFALSE);
3220 0 : SetRequireTPCRefit(kTRUE);
3221 :
3222 : //ITS
3223 0 : SetRequireITSRefit(kTRUE);
3224 0 : SetMaxDCAToVertexZ(2);
3225 0 : SetMaxDCAToVertexXY();
3226 0 : SetDCAToVertex2D(kFALSE);
3227 :
3228 0 : break;
3229 : }
3230 :
3231 : case kLooseCut:
3232 : {
3233 0 : AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3234 0 : SetMinNClustersTPC(50);
3235 0 : SetAcceptKinkDaughters(kTRUE);
3236 :
3237 0 : break;
3238 : }
3239 :
3240 : case kITSStandAlone:
3241 : {
3242 0 : AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3243 0 : SetRequireITSRefit(kTRUE);
3244 0 : SetRequireITSStandAlone(kTRUE);
3245 0 : SetITSTrackSA(kTRUE);
3246 0 : break;
3247 : }
3248 :
3249 : case kGlobalCut2011:
3250 : {
3251 0 : AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3252 : //TPC
3253 0 : SetMinNClustersTPC(50);
3254 0 : SetMaxChi2PerClusterTPC(4);
3255 0 : SetAcceptKinkDaughters(kFALSE);
3256 0 : SetRequireTPCRefit(kTRUE);
3257 :
3258 : //ITS
3259 0 : SetRequireITSRefit(kTRUE);
3260 0 : SetMaxDCAToVertexZ(2);
3261 0 : SetMaxDCAToVertexXY();
3262 0 : SetDCAToVertex2D(kFALSE);
3263 :
3264 0 : break;
3265 : }
3266 :
3267 : case kLooseCutWithITSrefit:
3268 : {
3269 0 : AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3270 0 : SetMinNClustersTPC(50);
3271 0 : SetAcceptKinkDaughters(kTRUE);
3272 0 : SetRequireITSRefit(kTRUE);
3273 :
3274 0 : break;
3275 : }
3276 : }
3277 0 : }
3278 :
3279 : ///
3280 : /// Check if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
3281 : ///
3282 : /// \param event: pointer to event
3283 : ///
3284 : //________________________________________________________________________
3285 : void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
3286 : {
3287 0 : Int_t nTracks = event->GetNumberOfTracks();
3288 0 : for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3289 : {
3290 0 : AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3291 0 : if (!track)
3292 : {
3293 0 : AliWarning(Form("Could not receive track %d", iTrack));
3294 0 : continue;
3295 : }
3296 :
3297 0 : Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3298 0 : track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3299 :
3300 : /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3301 0 : AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3302 0 : if (esdtrack)
3303 : {
3304 0 : if (matchClusIndex != -1)
3305 0 : esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3306 : else
3307 0 : esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3308 : }
3309 : else
3310 : {
3311 0 : AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3312 0 : if (matchClusIndex != -1)
3313 0 : aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3314 : else
3315 0 : aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3316 : }
3317 0 : }
3318 0 : AliDebug(2,"Track matched to closest cluster");
3319 0 : }
3320 :
3321 : ///
3322 : /// Checks if EMC clusters are matched to ESD track.
3323 : /// Adds track indexes of all the tracks matched to a cluster within residuals in AliVCluster.
3324 : ///
3325 : /// \param event: pointer to event
3326 : ///
3327 : //_________________________________________________________________________
3328 : void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
3329 : {
3330 0 : for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
3331 : {
3332 0 : AliVCluster *cluster = event->GetCaloCluster(iClus);
3333 0 : if (!cluster->IsEMCAL())
3334 0 : continue;
3335 :
3336 : //
3337 : // Remove old matches in cluster
3338 : //
3339 0 : if(cluster->GetNTracksMatched() > 0)
3340 : {
3341 0 : if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
3342 : {
3343 0 : TArrayI arrayTrackMatched(0);
3344 0 : ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
3345 0 : }
3346 : else
3347 : {
3348 0 : for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
3349 : {
3350 0 : ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
3351 : }
3352 : }
3353 : }
3354 :
3355 : //
3356 : // Find new matches and put them in the cluster
3357 : //
3358 0 : Int_t nTracks = event->GetNumberOfTracks();
3359 0 : TArrayI arrayTrackMatched(nTracks);
3360 :
3361 : // Get the closest track matched to the cluster
3362 : Int_t nMatched = 0;
3363 0 : Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
3364 0 : if (matchTrackIndex != -1)
3365 : {
3366 0 : arrayTrackMatched[nMatched] = matchTrackIndex;
3367 : nMatched++;
3368 0 : }
3369 :
3370 : // Get all other tracks matched to the cluster
3371 0 : for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
3372 : {
3373 0 : AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
3374 :
3375 0 : if( !track ) continue;
3376 :
3377 0 : if ( iTrk == matchTrackIndex ) continue;
3378 :
3379 0 : if ( track->GetEMCALcluster() == iClus )
3380 : {
3381 0 : arrayTrackMatched[nMatched] = iTrk;
3382 0 : ++nMatched;
3383 0 : }
3384 0 : }
3385 :
3386 0 : AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
3387 :
3388 0 : arrayTrackMatched.Set(nMatched);
3389 0 : AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
3390 0 : if (esdcluster)
3391 0 : esdcluster->AddTracksMatched(arrayTrackMatched);
3392 0 : else if ( nMatched > 0 )
3393 : {
3394 0 : AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
3395 0 : if (aodcluster)
3396 : {
3397 0 : aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
3398 : //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
3399 : //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());
3400 : //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
3401 : }
3402 0 : }
3403 :
3404 0 : Float_t eta= -999, phi = -999;
3405 0 : if (matchTrackIndex != -1)
3406 0 : GetMatchedResiduals(iClus, eta, phi);
3407 :
3408 0 : cluster->SetTrackDistance(phi, eta);
3409 0 : }
3410 :
3411 0 : AliDebug(2,"Cluster matched to tracks");
3412 0 : }
3413 :
3414 : ///
3415 : /// Print Parameters.
3416 : ///
3417 : //___________________________________________________
3418 : void AliEMCALRecoUtils::Print(const Option_t *) const
3419 : {
3420 0 : printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3421 0 : printf("AliEMCALRecoUtils Settings: \n");
3422 0 : printf("\tMisalignment shifts\n");
3423 0 : 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,
3424 0 : fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
3425 0 : fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
3426 0 : printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
3427 0 : if (fNonLinearityFunction != 3) // print only if not kNoCorrection
3428 0 : for (Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
3429 :
3430 0 : printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
3431 :
3432 0 : printf("\tMatching criteria: ");
3433 0 : if (fCutEtaPhiSum) {
3434 0 : printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
3435 0 : } else if (fCutEtaPhiSeparate) {
3436 0 : printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
3437 0 : } else {
3438 0 : printf("\t\tError\n");
3439 0 : printf("\t\tplease specify your cut criteria\n");
3440 0 : printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
3441 0 : printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
3442 : }
3443 :
3444 0 : 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);
3445 0 : printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
3446 :
3447 0 : printf("\tTrack cuts: \n");
3448 0 : printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
3449 0 : printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
3450 0 : printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
3451 0 : printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
3452 0 : printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
3453 0 : printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
3454 0 : printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
3455 0 : printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3456 0 : }
|
