]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/EMCALUtils/AliEMCALRecoUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / EMCAL / EMCALUtils / AliEMCALRecoUtils.cxx
CommitLineData
d9b3567c 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
40891fa2 16/* $Id: AliEMCALRecoUtils.cxx | Sun Dec 8 06:56:48 2013 +0100 | Constantin Loizides $ */
d9b3567c 17
18///////////////////////////////////////////////////////////////////////////////
19//
20// Class AliEMCALRecoUtils
21// Some utilities to recalculate the cluster position or energy linearity
22//
23//
24// Author: Gustavo Conesa (LPSC- Grenoble)
b540d03f 25// Track matching part: Rongrong Ma (Yale)
26
d9b3567c 27///////////////////////////////////////////////////////////////////////////////
d9b3567c 28// --- standard c ---
29
30// standard C++ includes
31//#include <Riostream.h>
32
33// ROOT includes
094786cc 34#include <TGeoManager.h>
35#include <TGeoMatrix.h>
36#include <TGeoBBox.h>
7cdec71f 37#include <TH2F.h>
38#include <TArrayI.h>
39#include <TArrayF.h>
01d44f1f 40#include <TObjArray.h>
d9b3567c 41
42// STEER includes
d9b3567c 43#include "AliVCluster.h"
44#include "AliVCaloCells.h"
45#include "AliLog.h"
83bfd77a 46#include "AliPID.h"
a520bcd0 47#include "AliESDEvent.h"
bb6f5f0b 48#include "AliAODEvent.h"
bd8c7aef 49#include "AliESDtrack.h"
bb6f5f0b 50#include "AliAODTrack.h"
51#include "AliExternalTrackParam.h"
52#include "AliESDfriendTrack.h"
53#include "AliTrackerBase.h"
b540d03f 54
55// EMCAL includes
56#include "AliEMCALRecoUtils.h"
57#include "AliEMCALGeometry.h"
ee602376 58#include "AliTrackerBase.h"
b540d03f 59#include "AliEMCALPIDUtils.h"
ee602376 60
d9b3567c 61ClassImp(AliEMCALRecoUtils)
62
88b96ad8 63//_____________________________________
d9b3567c 64AliEMCALRecoUtils::AliEMCALRecoUtils():
88b96ad8 65 fParticleType(0), fPosAlgo(0), fW0(0),
66 fNonLinearityFunction(0), fNonLinearThreshold(0),
01d44f1f 67 fSmearClusterEnergy(kFALSE), fRandom(),
3bfc4732 68 fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
7bf608c9 69 fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fUseRunCorrectionFactors(kFALSE),
01d44f1f 70 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
71 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
a7e5a381 72 fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
88b96ad8 73 fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
a6a1e3ab 74 fPIDUtils(), fAODFilterMask(0),
75 fAODHybridTracks(0), fAODTPCOnlyTracks(0),
01d44f1f 76 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
88b96ad8 77 fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
78 fCutR(0), fCutEta(0), fCutPhi(0),
79 fClusterWindow(0), fMass(0),
80 fStepSurface(0), fStepCluster(0),
8517fbd8 81 fITSTrackSA(kFALSE), fEMCalSurfaceDistance(440.),
88b96ad8 82 fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
83 fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
01d44f1f 84 fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
42ceff04 85 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
18a566cc 86 fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE)
d9b3567c 87{
88//
89 // Constructor.
90 // Initialize all constant values which have to be used
91 // during Reco algorithm execution
92 //
93
88b96ad8 94 // Init parameters
95 InitParameters();
01d44f1f 96
b540d03f 97 //Track matching
7cdec71f 98 fMatchedTrackIndex = new TArrayI();
99 fMatchedClusterIndex = new TArrayI();
bd36717e 100 fResidualPhi = new TArrayF();
101 fResidualEta = new TArrayF();
7cdec71f 102 fPIDUtils = new AliEMCALPIDUtils();
01d44f1f 103
d9b3567c 104}
105
106//______________________________________________________________________
107AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
01d44f1f 108: TNamed(reco),
109 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
110 fNonLinearityFunction(reco.fNonLinearityFunction), fNonLinearThreshold(reco.fNonLinearThreshold),
111 fSmearClusterEnergy(reco.fSmearClusterEnergy), fRandom(),
3bfc4732 112 fCellsRecalibrated(reco.fCellsRecalibrated),
01d44f1f 113 fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
3bfc4732 114 fTimeRecalibration(reco.fTimeRecalibration), fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
7bf608c9 115 fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),
01d44f1f 116 fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
78467229 117 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
01d44f1f 118 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
a7e5a381 119 fRejectExoticCluster(reco.fRejectExoticCluster), fRejectExoticCells(reco.fRejectExoticCells),
120 fExoticCellFraction(reco.fExoticCellFraction), fExoticCellDiffTime(reco.fExoticCellDiffTime),
121 fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
122 fPIDUtils(reco.fPIDUtils), fAODFilterMask(reco.fAODFilterMask),
a6a1e3ab 123 fAODHybridTracks(reco.fAODHybridTracks), fAODTPCOnlyTracks(reco.fAODTPCOnlyTracks),
01d44f1f 124 fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
bd8c7aef 125 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
01d44f1f 126 fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
127 fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
128 fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
129 fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
8fc351e3 130 fClusterWindow(reco.fClusterWindow),
131 fMass(reco.fMass), fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
8517fbd8 132 fITSTrackSA(reco.fITSTrackSA), fEMCalSurfaceDistance(440.),
01d44f1f 133 fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
134 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
135 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
136 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
137 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
42ceff04 138 fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
18a566cc 139 fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone), fCutRequireITSpureSA(reco.fCutRequireITSpureSA)
d9b3567c 140{
141 //Copy ctor
142
8517fbd8 143 for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
a31af82c 144 fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
8517fbd8 145 for (Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
146 for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
bd8c7aef 147
d9b3567c 148}
149
150
151//______________________________________________________________________
152AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
153{
154 //Assignment operator
155
8517fbd8 156 if (this == &reco)return *this;
d9b3567c 157 ((TNamed *)this)->operator=(reco);
158
8517fbd8 159 for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
01d44f1f 160 fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
8517fbd8 161 for (Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
162 for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
01d44f1f 163
96957075 164 fParticleType = reco.fParticleType;
165 fPosAlgo = reco.fPosAlgo;
166 fW0 = reco.fW0;
01d44f1f 167
168 fNonLinearityFunction = reco.fNonLinearityFunction;
7e0ecb89 169 fNonLinearThreshold = reco.fNonLinearThreshold;
01d44f1f 170 fSmearClusterEnergy = reco.fSmearClusterEnergy;
171
3bfc4732 172 fCellsRecalibrated = reco.fCellsRecalibrated;
96957075 173 fRecalibration = reco.fRecalibration;
094786cc 174 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
3bfc4732 175
176 fTimeRecalibration = reco.fTimeRecalibration;
177 fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
178
179 fUseRunCorrectionFactors = reco.fUseRunCorrectionFactors;
01d44f1f 180
96957075 181 fRemoveBadChannels = reco.fRemoveBadChannels;
182 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
183 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
01d44f1f 184
96957075 185 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
186 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
a7e5a381 187
01d44f1f 188 fRejectExoticCluster = reco.fRejectExoticCluster;
a7e5a381 189 fRejectExoticCells = reco.fRejectExoticCells;
190 fExoticCellFraction = reco.fExoticCellFraction;
191 fExoticCellDiffTime = reco.fExoticCellDiffTime;
192 fExoticCellMinAmplitude = reco.fExoticCellMinAmplitude;
193
01d44f1f 194 fPIDUtils = reco.fPIDUtils;
83bfd77a 195
01d44f1f 196 fAODFilterMask = reco.fAODFilterMask;
1af378e6 197 fAODHybridTracks = reco.fAODHybridTracks;
a6a1e3ab 198 fAODTPCOnlyTracks = reco.fAODTPCOnlyTracks;
d9b3567c 199
fa4287a2 200 fCutEtaPhiSum = reco.fCutEtaPhiSum;
201 fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
96957075 202 fCutR = reco.fCutR;
fa4287a2 203 fCutEta = reco.fCutEta;
204 fCutPhi = reco.fCutPhi;
8fc351e3 205 fClusterWindow = reco.fClusterWindow;
bb6f5f0b 206 fMass = reco.fMass;
8fc351e3 207 fStepSurface = reco.fStepSurface;
208 fStepCluster = reco.fStepCluster;
a29b2a8a 209 fITSTrackSA = reco.fITSTrackSA;
210 fEMCalSurfaceDistance = reco.fEMCalSurfaceDistance;
42ceff04 211
5f7714ad 212 fTrackCutsType = reco.fTrackCutsType;
fa4287a2 213 fCutMinTrackPt = reco.fCutMinTrackPt;
96957075 214 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
215 fCutMinNClusterITS = reco.fCutMinNClusterITS;
216 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
217 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
218 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
219 fCutRequireITSRefit = reco.fCutRequireITSRefit;
220 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
221 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
222 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
223 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
42ceff04 224 fCutRequireITSStandAlone = reco.fCutRequireITSStandAlone;
18a566cc 225 fCutRequireITSpureSA = reco.fCutRequireITSpureSA;
226
8517fbd8 227 if (reco.fResidualEta) {
bd8c7aef 228 // assign or copy construct
8517fbd8 229 if (fResidualEta) {
fa4287a2 230 *fResidualEta = *reco.fResidualEta;
8517fbd8 231 } else {
7f5392da 232 fResidualEta = new TArrayF(*reco.fResidualEta);
233 }
8517fbd8 234 } else {
235 if (fResidualEta) delete fResidualEta;
fa4287a2 236 fResidualEta = 0;
bd8c7aef 237 }
238
8517fbd8 239 if (reco.fResidualPhi) {
bd8c7aef 240 // assign or copy construct
8517fbd8 241 if (fResidualPhi) {
fa4287a2 242 *fResidualPhi = *reco.fResidualPhi;
8517fbd8 243 } else {
7f5392da 244 fResidualPhi = new TArrayF(*reco.fResidualPhi);
245 }
8517fbd8 246 } else {
247 if (fResidualPhi) delete fResidualPhi;
fa4287a2 248 fResidualPhi = 0;
bd8c7aef 249 }
250
8517fbd8 251 if (reco.fMatchedTrackIndex) {
b540d03f 252 // assign or copy construct
8517fbd8 253 if (fMatchedTrackIndex) {
b540d03f 254 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
8517fbd8 255 } else {
7f5392da 256 fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
257 }
8517fbd8 258 } else {
259 if (fMatchedTrackIndex) delete fMatchedTrackIndex;
b540d03f 260 fMatchedTrackIndex = 0;
261 }
bd8c7aef 262
8517fbd8 263 if (reco.fMatchedClusterIndex) {
bd8c7aef 264 // assign or copy construct
8517fbd8 265 if (fMatchedClusterIndex) {
bd8c7aef 266 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
8517fbd8 267 } else {
7f5392da 268 fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
269 }
8517fbd8 270 } else {
271 if (fMatchedClusterIndex) delete fMatchedClusterIndex;
bd8c7aef 272 fMatchedClusterIndex = 0;
273 }
18a566cc 274
d9b3567c 275 return *this;
276}
277
a7e5a381 278//_____________________________________
094786cc 279AliEMCALRecoUtils::~AliEMCALRecoUtils()
280{
281 //Destructor.
fbddd006 282
8517fbd8 283 if (fEMCALRecalibrationFactors) {
b6557fd1 284 fEMCALRecalibrationFactors->Clear();
841dbf60 285 delete fEMCALRecalibrationFactors;
fbddd006 286 }
fd6df01c 287
8517fbd8 288 if (fEMCALTimeRecalibrationFactors) {
841dbf60 289 fEMCALTimeRecalibrationFactors->Clear();
290 delete fEMCALTimeRecalibrationFactors;
fbddd006 291 }
3bfc4732 292
8517fbd8 293 if (fEMCALBadChannelMap) {
b6557fd1 294 fEMCALBadChannelMap->Clear();
841dbf60 295 delete fEMCALBadChannelMap;
b6557fd1 296 }
bd8c7aef 297
7cdec71f 298 delete fMatchedTrackIndex ;
299 delete fMatchedClusterIndex ;
300 delete fResidualEta ;
301 delete fResidualPhi ;
b6557fd1 302 delete fPIDUtils ;
bd8c7aef 303
b6557fd1 304 InitTrackCuts();
094786cc 305}
306
a7e5a381 307//_______________________________________________________________________________
fb8eab96 308Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(Int_t absID, Int_t bc,
a7e5a381 309 Float_t & amp, Double_t & time,
310 AliVCaloCells* cells)
311{
312 // Reject cell if criteria not passed and calibrate it
313
314 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
315
8517fbd8 316 if (absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules())
317 return kFALSE;
a7e5a381 318
319 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
c3d9f926 320
8517fbd8 321 if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta)) {
c3d9f926 322 // cell absID does not exist
323 amp=0; time = 1.e9;
324 return kFALSE;
325 }
326
fbddd006 327 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
84cb7b3f 328
a7e5a381 329 // Do not include bad channels found in analysis,
8517fbd8 330 if (IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) {
a7e5a381 331 return kFALSE;
332 }
333
334 //Recalibrate energy
335 amp = cells->GetCellAmplitude(absID);
8517fbd8 336 if (!fCellsRecalibrated && IsRecalibrationOn())
a7e5a381 337 amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
338
a7e5a381 339 // Recalibrate time
340 time = cells->GetCellTime(absID);
341
342 RecalibrateCellTime(absID,bc,time);
343
344 return kTRUE;
345}
346
a520bcd0 347//_____________________________________________________________________________
348Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom,
349 const AliVCluster* cluster,
350 AliVCaloCells* cells)
fd6df01c 351{
fbddd006 352 // Given the list of AbsId of the cluster, get the maximum cell and
353 // check if there are fNCellsFromBorder from the calorimeter border
354
8517fbd8 355 if (!cluster)
7f5392da 356 {
2aeb4226 357 AliInfo("Cluster pointer null!");
358 return kFALSE;
359 }
360
fd6df01c 361 //If the distance to the border is 0 or negative just exit accept all clusters
8517fbd8 362 if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
363 return kTRUE;
fd6df01c 364
fbddd006 365 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
cb231979 366 Bool_t shared = kFALSE;
367 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
7f5392da 368
83bfd77a 369 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
7f5392da 370 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
fbddd006 371
8517fbd8 372 if (absIdMax==-1) return kFALSE;
fbddd006 373
374 //Check if the cell is close to the borders:
375 Bool_t okrow = kFALSE;
376 Bool_t okcol = kFALSE;
fd6df01c 377
8517fbd8 378 if (iSM < 0 || iphi < 0 || ieta < 0 ) {
fd6df01c 379 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
380 iSM,ieta,iphi));
c9a7e94b 381 return kFALSE; // trick coverity
fd6df01c 382 }
383
384 //Check rows/phi
8cc543cb 385 Int_t iPhiLast = 24;
386 if( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
387 else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
fd6df01c 388
8cc543cb 389 if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
390
fd6df01c 391 //Check columns/eta
8cc543cb 392 Int_t iEtaLast = 48;
393 if(!fNoEMCALBorderAtEta0 || geom->IsDCALSM(iSM)) {// conside inner border
394 if( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard ) iEtaLast = iEtaLast*2/3;
395 if(ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
8517fbd8 396 } else {
397 if (iSM%2==0) {
8cc543cb 398 if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
8517fbd8 399 } else {
8cc543cb 400 if(ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
fd6df01c 401 }
402 }//eta 0 not checked
7f5392da 403
83bfd77a 404 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
7f5392da 405 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
fbddd006 406
8517fbd8 407 if (okcol && okrow) {
83bfd77a 408 //printf("Accept\n");
409 return kTRUE;
8517fbd8 410 } else {
83bfd77a 411 //printf("Reject\n");
412 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
413 return kFALSE;
414 }
fbddd006 415}
fd6df01c 416
a520bcd0 417//_______________________________________________________________________________
418Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom,
419 const UShort_t* cellList,
fb8eab96 420 Int_t nCells)
841dbf60 421{
422 // Check that in the cluster cells, there is no bad channel of those stored
423 // in fEMCALBadChannelMap or fPHOSBadChannelMap
fbddd006 424
8517fbd8 425 if (!fRemoveBadChannels) return kFALSE;
426 if (!fEMCALBadChannelMap) return kFALSE;
fbddd006 427
841dbf60 428 Int_t icol = -1;
429 Int_t irow = -1;
430 Int_t imod = -1;
8517fbd8 431 for (Int_t iCell = 0; iCell<nCells; iCell++) {
841dbf60 432 //Get the column and row
fd6df01c 433 Int_t iTower = -1, iIphi = -1, iIeta = -1;
434 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
8517fbd8 435 if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
fbddd006 436 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
8517fbd8 437 if (GetEMCALChannelStatus(imod, icol, irow)) {
83bfd77a 438 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
439 return kTRUE;
440 }
841dbf60 441 }// cell cluster loop
8517fbd8 442
841dbf60 443 return kFALSE;
fd6df01c 444}
094786cc 445
841dbf60 446
ba19aaf1 447//___________________________________________________________________________
fb8eab96 448Float_t AliEMCALRecoUtils::GetECross(Int_t absID, Double_t tcell,
449 AliVCaloCells* cells, Int_t bc)
ba19aaf1 450{
451 //Calculate the energy in the cross around the energy given cell
7f5392da 452
a7e5a381 453 AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
454
455 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
456 geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
fbddd006 457 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
a7e5a381 458
459 //Get close cells index, energy and time, not in corners
ba19aaf1 460
84cb7b3f 461 Int_t absID1 = -1;
462 Int_t absID2 = -1;
7f5392da 463
8517fbd8 464 if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
465 if ( iphi > 0 ) absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
a7e5a381 466
7f5392da 467 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
ba19aaf1 468
7f5392da 469 Int_t absID3 = -1;
470 Int_t absID4 = -1;
471
8517fbd8 472 if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
84cb7b3f 473 absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
474 absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
8517fbd8 475 } else if ( ieta == 0 && imod%2 ) {
84cb7b3f 476 absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
477 absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
8517fbd8 478 } else {
479 if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
84cb7b3f 480 absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
8517fbd8 481 if ( ieta > 0 )
84cb7b3f 482 absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
7f5392da 483 }
a7e5a381 484
ba19aaf1 485 //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
a7e5a381 486
ba19aaf1 487 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
488 Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
8517fbd8 489
490 AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
491 AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
492 AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
493 AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
ba19aaf1 494
8517fbd8 495 if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
496 if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
497 if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
498 if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
ba19aaf1 499
8517fbd8 500 return ecell1+ecell2+ecell3+ecell4;
ba19aaf1 501}
a7e5a381 502
ba19aaf1 503//_____________________________________________________________________________________________
fb8eab96 504Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
ba19aaf1 505{
506 // Look to cell neighbourhood and reject if it seems exotic
507 // Do before recalibrating the cells
508
8517fbd8 509 if (!fRejectExoticCells) return kFALSE;
ba19aaf1 510
511 Float_t ecell = 0;
512 Double_t tcell = 0;
513 Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
514
8517fbd8 515 if (!accept) return kTRUE; // reject this cell
ba19aaf1 516
8517fbd8 517 if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
a7e5a381 518
ba19aaf1 519 Float_t eCross = GetECross(absID,tcell,cells,bc);
a7e5a381 520
8517fbd8 521 if (1-eCross/ecell > fExoticCellFraction) {
841dbf60 522 AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
523 absID,ecell,eCross,1-eCross/ecell));
a7e5a381 524 return kTRUE;
525 }
841dbf60 526
a7e5a381 527 return kFALSE;
a7e5a381 528}
529
a520bcd0 530//___________________________________________________________________
531Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
532 AliVCaloCells *cells,
fb8eab96 533 Int_t bc)
841dbf60 534{
a7e5a381 535 // Check if the cluster highest energy tower is exotic
2aeb4226 536
8517fbd8 537 if (!cluster) {
2aeb4226 538 AliInfo("Cluster pointer null!");
539 return kFALSE;
540 }
45516c1f 541
8517fbd8 542 if (!fRejectExoticCluster) return kFALSE;
45516c1f 543
a7e5a381 544 // Get highest energy tower
545 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
546 Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
547 Bool_t shared = kFALSE;
548 GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
7f5392da 549
a7e5a381 550 return IsExoticCell(absId,cells,bc);
b5078f5d 551}
552
a520bcd0 553//_______________________________________________________________________
841dbf60 554Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
555{
01d44f1f 556 //In case of MC analysis, smear energy to match resolution/calibration in real data
557
8517fbd8 558 if (!cluster) {
01d44f1f 559 AliInfo("Cluster pointer null!");
560 return 0;
561 }
562
563 Float_t energy = cluster->E() ;
564 Float_t rdmEnergy = energy ;
8517fbd8 565 if (fSmearClusterEnergy) {
01d44f1f 566 rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
567 fSmearClusterParam[1] * energy +
568 fSmearClusterParam[2] );
569 AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
570 }
571
841dbf60 572 return rdmEnergy;
01d44f1f 573}
574
a520bcd0 575//____________________________________________________________________________
841dbf60 576Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
577{
578 // Correct cluster energy from non linearity functions
2aeb4226 579
8517fbd8 580 if (!cluster) {
2aeb4226 581 AliInfo("Cluster pointer null!");
582 return 0;
583 }
584
d9b3567c 585 Float_t energy = cluster->E();
57131575 586
8517fbd8 587 if (energy < 0.05) {
24f680c7 588 // Clusters with less than 50 MeV or negative are not possible
589 AliInfo(Form("Too Low Cluster energy!, E = %f < 0.05 GeV",energy));
5890e234 590 return 0;
591 }
592
7f5392da 593 switch (fNonLinearityFunction)
594 {
d9b3567c 595 case kPi0MC:
871aee7a 596 {
d9b3567c 597 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
4439cfaf 598 //fNonLinearityParams[0] = 1.014;
599 //fNonLinearityParams[1] =-0.03329;
600 //fNonLinearityParams[2] =-0.3853;
601 //fNonLinearityParams[3] = 0.5423;
602 //fNonLinearityParams[4] =-0.4335;
8cdd1f1f 603 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
d9b3567c 604 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
605 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
606 break;
871aee7a 607 }
dff9e2e3 608
4439cfaf 609 case kPi0MCv2:
610 {
611 //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
612 //fNonLinearityParams[0] = 3.11111e-02;
613 //fNonLinearityParams[1] =-5.71666e-02;
614 //fNonLinearityParams[2] = 5.67995e-01;
615
616 energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
617 break;
618 }
e33e0c4a 619
620 case kPi0MCv3:
621 {
622 //Same as beam test corrected, change parameters
623 //fNonLinearityParams[0] = 9.81039e-01
624 //fNonLinearityParams[1] = 1.13508e-01;
625 //fNonLinearityParams[2] = 1.00173e+00;
626 //fNonLinearityParams[3] = 9.67998e-02;
627 //fNonLinearityParams[4] = 2.19381e+02;
628 //fNonLinearityParams[5] = 6.31604e+01;
629 //fNonLinearityParams[6] = 1;
630 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
631
632 break;
633 }
634
4439cfaf 635
d9b3567c 636 case kPi0GammaGamma:
871aee7a 637 {
d9b3567c 638 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
4439cfaf 639 //fNonLinearityParams[0] = 1.04;
640 //fNonLinearityParams[1] = -0.1445;
641 //fNonLinearityParams[2] = 1.046;
d9b3567c 642 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
643 break;
871aee7a 644 }
d9b3567c 645
646 case kPi0GammaConversion:
871aee7a 647 {
d9b3567c 648 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
871aee7a 649 //fNonLinearityParams[0] = 0.139393/0.1349766;
650 //fNonLinearityParams[1] = 0.0566186;
651 //fNonLinearityParams[2] = 0.982133;
d9b3567c 652 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
653
654 break;
871aee7a 655 }
656
657 case kBeamTest:
658 {
659 //From beam test, Alexei's results, for different ZS thresholds
660 // th=30 MeV; th = 45 MeV; th = 75 MeV
96957075 661 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
871aee7a 662 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
663 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
96957075 664 //Rescale the param[0] with 1.03
871aee7a 665 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
666
667 break;
668 }
dff9e2e3 669
4b58ac4f 670 case kBeamTestCorrected:
671 {
672 //From beam test, corrected for material between beam and EMCAL
dff9e2e3 673 //fNonLinearityParams[0] = 0.99078
674 //fNonLinearityParams[1] = 0.161499;
675 //fNonLinearityParams[2] = 0.655166;
676 //fNonLinearityParams[3] = 0.134101;
677 //fNonLinearityParams[4] = 163.282;
678 //fNonLinearityParams[5] = 23.6904;
679 //fNonLinearityParams[6] = 0.978;
680 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
681
4b58ac4f 682 break;
683 }
a31af82c 684
685 case kBeamTestCorrectedv2:
686 {
687 //From beam test, corrected for material between beam and EMCAL
688 //fNonLinearityParams[0] = 0.983504;
689 //fNonLinearityParams[1] = 0.210106;
690 //fNonLinearityParams[2] = 0.897274;
691 //fNonLinearityParams[3] = 0.0829064;
692 //fNonLinearityParams[4] = 152.299;
693 //fNonLinearityParams[5] = 31.5028;
694 //fNonLinearityParams[6] = 0.968;
695 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
696
697 break;
698 }
d9b3567c 699
d24a0a76 700 case kSDMv5:
701 {
702 //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
703 //fNonLinearityParams[0] = 1.0;
704 //fNonLinearityParams[1] = 6.64778e-02;
705 //fNonLinearityParams[2] = 1.570;
706 //fNonLinearityParams[3] = 9.67998e-02;
707 //fNonLinearityParams[4] = 2.19381e+02;
708 //fNonLinearityParams[5] = 6.31604e+01;
709 //fNonLinearityParams[6] = 1.01286;
710 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));
711
712 break;
713 }
714
715 case kPi0MCv5:
716 {
717 //Based on comparing MC truth information to the reconstructed energy of clusters.
718 //fNonLinearityParams[0] = 1.0;
719 //fNonLinearityParams[1] = 6.64778e-02;
720 //fNonLinearityParams[2] = 1.570;
721 //fNonLinearityParams[3] = 9.67998e-02;
722 //fNonLinearityParams[4] = 2.19381e+02;
723 //fNonLinearityParams[5] = 6.31604e+01;
724 //fNonLinearityParams[6] = 1.01286;
725 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
726
727 break;
728 }
729
4c489e4e 730 case kSDMv6:
731 {
732 //Based on fit to the MC/data using kNoCorrection on the data
733 // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
734 // - parameters constrained by the test beam data as well
735 // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
736 //fNonLinearityParams[0] = 1.0;
737 //fNonLinearityParams[1] = 0.237767;
738 //fNonLinearityParams[2] = 0.651203;
739 //fNonLinearityParams[3] = 0.183741;
740 //fNonLinearityParams[4] = 155.427;
741 //fNonLinearityParams[5] = 17.0335;
742 //fNonLinearityParams[6] = 0.987054;
743 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
744
745 break;
746 }
747
748 case kPi0MCv6:
749 {
750 //Based on comparing MC truth information to the reconstructed energy of clusters.
751 // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
752 //fNonLinearityParams[0] = 1.0;
753 //fNonLinearityParams[1] = 0.0797873;
754 //fNonLinearityParams[2] = 1.68322;
755 //fNonLinearityParams[3] = 0.0806098;
756 //fNonLinearityParams[4] = 244.586;
757 //fNonLinearityParams[5] = 116.938;
758 //fNonLinearityParams[6] = 1.00437;
759 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
760
761 break;
762 }
763
d9b3567c 764 case kNoCorrection:
765 AliDebug(2,"No correction on the energy\n");
766 break;
767
768 }
57131575 769
d9b3567c 770 return energy;
d9b3567c 771}
841dbf60 772
7e0ecb89 773//__________________________________________________
774void AliEMCALRecoUtils::InitNonLinearityParam()
775{
841dbf60 776 //Initialising Non Linearity Parameters
fbddd006 777
8517fbd8 778 if (fNonLinearityFunction == kPi0MC) {
841dbf60 779 fNonLinearityParams[0] = 1.014;
780 fNonLinearityParams[1] = -0.03329;
781 fNonLinearityParams[2] = -0.3853;
782 fNonLinearityParams[3] = 0.5423;
783 fNonLinearityParams[4] = -0.4335;
784 }
785
8517fbd8 786 if (fNonLinearityFunction == kPi0MCv2) {
e33e0c4a 787 fNonLinearityParams[0] = 3.11111e-02;
788 fNonLinearityParams[1] =-5.71666e-02;
789 fNonLinearityParams[2] = 5.67995e-01;
790 }
791
8517fbd8 792 if (fNonLinearityFunction == kPi0MCv3) {
e33e0c4a 793 fNonLinearityParams[0] = 9.81039e-01;
794 fNonLinearityParams[1] = 1.13508e-01;
795 fNonLinearityParams[2] = 1.00173e+00;
796 fNonLinearityParams[3] = 9.67998e-02;
797 fNonLinearityParams[4] = 2.19381e+02;
798 fNonLinearityParams[5] = 6.31604e+01;
799 fNonLinearityParams[6] = 1;
800 }
801
8517fbd8 802 if (fNonLinearityFunction == kPi0GammaGamma) {
841dbf60 803 fNonLinearityParams[0] = 1.04;
804 fNonLinearityParams[1] = -0.1445;
805 fNonLinearityParams[2] = 1.046;
fbddd006 806 }
841dbf60 807
8517fbd8 808 if (fNonLinearityFunction == kPi0GammaConversion) {
841dbf60 809 fNonLinearityParams[0] = 0.139393;
810 fNonLinearityParams[1] = 0.0566186;
811 fNonLinearityParams[2] = 0.982133;
fbddd006 812 }
841dbf60 813
8517fbd8 814 if (fNonLinearityFunction == kBeamTest) {
815 if (fNonLinearThreshold == 30) {
841dbf60 816 fNonLinearityParams[0] = 1.007;
817 fNonLinearityParams[1] = 0.894;
818 fNonLinearityParams[2] = 0.246;
819 }
8517fbd8 820 if (fNonLinearThreshold == 45) {
841dbf60 821 fNonLinearityParams[0] = 1.003;
822 fNonLinearityParams[1] = 0.719;
823 fNonLinearityParams[2] = 0.334;
824 }
8517fbd8 825 if (fNonLinearThreshold == 75) {
841dbf60 826 fNonLinearityParams[0] = 1.002;
827 fNonLinearityParams[1] = 0.797;
828 fNonLinearityParams[2] = 0.358;
829 }
830 }
831
8517fbd8 832 if (fNonLinearityFunction == kBeamTestCorrected) {
841dbf60 833 fNonLinearityParams[0] = 0.99078;
834 fNonLinearityParams[1] = 0.161499;
835 fNonLinearityParams[2] = 0.655166;
836 fNonLinearityParams[3] = 0.134101;
837 fNonLinearityParams[4] = 163.282;
838 fNonLinearityParams[5] = 23.6904;
839 fNonLinearityParams[6] = 0.978;
840 }
a31af82c 841
8517fbd8 842 if (fNonLinearityFunction == kBeamTestCorrectedv2) {
a31af82c 843 fNonLinearityParams[0] = 0.983504;
844 fNonLinearityParams[1] = 0.210106;
845 fNonLinearityParams[2] = 0.897274;
846 fNonLinearityParams[3] = 0.0829064;
847 fNonLinearityParams[4] = 152.299;
848 fNonLinearityParams[5] = 31.5028;
849 fNonLinearityParams[6] = 0.968;
850 }
d24a0a76 851
852 if (fNonLinearityFunction == kSDMv5) {
853 fNonLinearityParams[0] = 1.0;
854 fNonLinearityParams[1] = 6.64778e-02;
855 fNonLinearityParams[2] = 1.570;
856 fNonLinearityParams[3] = 9.67998e-02;
857 fNonLinearityParams[4] = 2.19381e+02;
858 fNonLinearityParams[5] = 6.31604e+01;
859 fNonLinearityParams[6] = 1.01286;
860 }
861
862 if (fNonLinearityFunction == kPi0MCv5) {
863 fNonLinearityParams[0] = 1.0;
864 fNonLinearityParams[1] = 6.64778e-02;
865 fNonLinearityParams[2] = 1.570;
866 fNonLinearityParams[3] = 9.67998e-02;
867 fNonLinearityParams[4] = 2.19381e+02;
868 fNonLinearityParams[5] = 6.31604e+01;
869 fNonLinearityParams[6] = 1.01286;
870 }
871
4c489e4e 872 if (fNonLinearityFunction == kSDMv6) {
873 fNonLinearityParams[0] = 1.0;
874 fNonLinearityParams[1] = 0.237767;
875 fNonLinearityParams[2] = 0.651203;
876 fNonLinearityParams[3] = 0.183741;
877 fNonLinearityParams[4] = 155.427;
878 fNonLinearityParams[5] = 17.0335;
879 fNonLinearityParams[6] = 0.987054;
880 }
881
882 if (fNonLinearityFunction == kPi0MCv6) {
883 fNonLinearityParams[0] = 1.0;
884 fNonLinearityParams[1] = 0.0797873;
885 fNonLinearityParams[2] = 1.68322;
886 fNonLinearityParams[3] = 0.0806098;
887 fNonLinearityParams[4] = 244.586;
888 fNonLinearityParams[5] = 116.938;
889 fNonLinearityParams[6] = 1.00437;
890 }
891
7e0ecb89 892}
893
a520bcd0 894//_________________________________________________________
fb8eab96 895Float_t AliEMCALRecoUtils::GetDepth(Float_t energy,
896 Int_t iParticle,
897 Int_t iSM) const
094786cc 898{
899 //Calculate shower depth for a given cluster energy and particle type
900
901 // parameters
cb231979 902 Float_t x0 = 1.31;
094786cc 903 Float_t ecr = 8;
904 Float_t depth = 0;
a840d589 905 Float_t arg = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV
094786cc 906
907 switch ( iParticle )
908 {
909 case kPhoton:
a840d589 910 if (arg < 1)
911 depth = 0;
912 else
913 depth = x0 * (TMath::Log(arg) + 0.5);
094786cc 914 break;
915
916 case kElectron:
a840d589 917 if (arg < 1)
918 depth = 0;
919 else
920 depth = x0 * (TMath::Log(arg) - 0.5);
094786cc 921 break;
922
923 case kHadron:
924 // hadron
925 // boxes anc. here
8517fbd8 926 if (gGeoManager) {
094786cc 927 gGeoManager->cd("ALIC_1/XEN1_1");
928 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
929 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
8517fbd8 930 if (geoSM) {
fd6df01c 931 TGeoVolume *geoSMVol = geoSM->GetVolume();
932 TGeoShape *geoSMShape = geoSMVol->GetShape();
933 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
8517fbd8 934 if (geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
fd6df01c 935 else AliFatal("Null GEANT box");
7f5392da 936 }
937 else AliFatal("NULL GEANT node matrix");
094786cc 938 }
7f5392da 939 else
940 {//electron
a840d589 941 if (arg < 1)
942 depth = 0;
943 else
944 depth = x0 * (TMath::Log(arg) - 0.5);
094786cc 945 }
946
947 break;
948
949 default://photon
a840d589 950 if (arg < 1)
951 depth = 0;
952 else
953 depth = x0 * (TMath::Log(arg) + 0.5);
094786cc 954 }
955
956 return depth;
094786cc 957}
958
88b96ad8 959//____________________________________________________________________
960void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
961 AliVCaloCells* cells,
962 const AliVCluster* clu,
a520bcd0 963 Int_t & absId,
964 Int_t & iSupMod,
965 Int_t & ieta,
966 Int_t & iphi,
967 Bool_t & shared)
d9b3567c 968{
969 //For a given CaloCluster gets the absId of the cell
970 //with maximum energy deposit.
971
972 Double_t eMax = -1.;
973 Double_t eCell = -1.;
094786cc 974 Float_t fraction = 1.;
975 Float_t recalFactor = 1.;
d9b3567c 976 Int_t cellAbsId = -1 ;
094786cc 977
d9b3567c 978 Int_t iTower = -1;
979 Int_t iIphi = -1;
980 Int_t iIeta = -1;
cb231979 981 Int_t iSupMod0= -1;
2aeb4226 982
8517fbd8 983 if (!clu) {
2aeb4226 984 AliInfo("Cluster pointer null!");
985 absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
986 return;
987 }
988
8517fbd8 989 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 990 cellAbsId = clu->GetCellAbsId(iDig);
991 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 992 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
8517fbd8 993 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 994 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
995 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
8517fbd8 996 if (iDig==0) {
7f5392da 997 iSupMod0=iSupMod;
8517fbd8 998 } else if (iSupMod0!=iSupMod) {
cb231979 999 shared = kTRUE;
1000 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1001 }
8517fbd8 1002 if (!fCellsRecalibrated && IsRecalibrationOn()) {
094786cc 1003 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1004 }
1005 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 1006 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
8517fbd8 1007 if (eCell > eMax) {
d9b3567c 1008 eMax = eCell;
1009 absId = cellAbsId;
1010 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1011 }
1012 }// cell loop
1013
1014 //Get from the absid the supermodule, tower and eta/phi numbers
1015 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1016 //Gives SuperModule and Tower numbers
1017 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 1018 iIphi, iIeta,iphi,ieta);
1019 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1020 //printf("Max end---\n");
d9b3567c 1021}
1022
88b96ad8 1023//______________________________________
1024void AliEMCALRecoUtils::InitParameters()
1025{
1026 // Initialize data members with default values
1027
1028 fParticleType = kPhoton;
1029 fPosAlgo = kUnchanged;
1030 fW0 = 4.5;
1031
1032 fNonLinearityFunction = kNoCorrection;
1033 fNonLinearThreshold = 30;
1034
1035 fExoticCellFraction = 0.97;
1036 fExoticCellDiffTime = 1e6;
31ba2b24 1037 fExoticCellMinAmplitude = 4.0;
88b96ad8 1038
a6a1e3ab 1039 fAODFilterMask = 128;
1040 fAODHybridTracks = kFALSE;
1041 fAODTPCOnlyTracks = kTRUE;
88b96ad8 1042
1043 fCutEtaPhiSum = kTRUE;
1044 fCutEtaPhiSeparate = kFALSE;
1045
1046 fCutR = 0.05;
1047 fCutEta = 0.025;
1048 fCutPhi = 0.05;
1049
1050 fClusterWindow = 100;
1051 fMass = 0.139;
1052
1053 fStepSurface = 20.;
1054 fStepCluster = 5.;
1055 fTrackCutsType = kLooseCut;
1056
1057 fCutMinTrackPt = 0;
1058 fCutMinNClusterTPC = -1;
1059 fCutMinNClusterITS = -1;
1060
1061 fCutMaxChi2PerClusterTPC = 1e10;
1062 fCutMaxChi2PerClusterITS = 1e10;
1063
1064 fCutRequireTPCRefit = kFALSE;
1065 fCutRequireITSRefit = kFALSE;
1066 fCutAcceptKinkDaughters = kFALSE;
1067
1068 fCutMaxDCAToVertexXY = 1e10;
1069 fCutMaxDCAToVertexZ = 1e10;
1070 fCutDCAToVertex2D = kFALSE;
1071
42ceff04 1072 fCutRequireITSStandAlone = kFALSE; //MARCEL
1073 fCutRequireITSpureSA = kFALSE; //Marcel
88b96ad8 1074
1075 //Misalignment matrices
8517fbd8 1076 for (Int_t i = 0; i < 15 ; i++)
7f5392da 1077 {
88b96ad8 1078 fMisalTransShift[i] = 0.;
1079 fMisalRotShift[i] = 0.;
1080 }
1081
1082 //Non linearity
8517fbd8 1083 for (Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
88b96ad8 1084
a31af82c 1085 //For kBeamTestCorrectedv2 case, but default is no correction
1086 fNonLinearityParams[0] = 0.983504;
1087 fNonLinearityParams[1] = 0.210106;
1088 fNonLinearityParams[2] = 0.897274;
1089 fNonLinearityParams[3] = 0.0829064;
1090 fNonLinearityParams[4] = 152.299;
1091 fNonLinearityParams[5] = 31.5028;
1092 fNonLinearityParams[6] = 0.968;
88b96ad8 1093
1094 //Cluster energy smearing
1095 fSmearClusterEnergy = kFALSE;
1096 fSmearClusterParam[0] = 0.07; // * sqrt E term
1097 fSmearClusterParam[1] = 0.00; // * E term
1098 fSmearClusterParam[2] = 0.00; // constant
88b96ad8 1099}
1100
1101//_____________________________________________________
1102void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
1103{
841dbf60 1104 //Init EMCAL recalibration factors
1105 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
fbddd006 1106 //In order to avoid rewriting the same histograms
841dbf60 1107 Bool_t oldStatus = TH1::AddDirectoryStatus();
1108 TH1::AddDirectory(kFALSE);
1109
c3d9f926 1110 fEMCALRecalibrationFactors = new TObjArray(12);
1111 for (int i = 0; i < 12; i++)
841dbf60 1112 fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1113 Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1114 //Init the histograms with 1
7f5392da 1115 for (Int_t sm = 0; sm < 12; sm++)
1116 {
1117 for (Int_t i = 0; i < 48; i++)
1118 {
1119 for (Int_t j = 0; j < 24; j++)
1120 {
841dbf60 1121 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1122 }
1123 }
1124 }
7f5392da 1125
841dbf60 1126 fEMCALRecalibrationFactors->SetOwner(kTRUE);
1127 fEMCALRecalibrationFactors->Compress();
fbddd006 1128
841dbf60 1129 //In order to avoid rewriting the same histograms
fbddd006 1130 TH1::AddDirectory(oldStatus);
094786cc 1131}
1132
a520bcd0 1133//_________________________________________________________
841dbf60 1134void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1135{
1136 //Init EMCAL recalibration factors
1137 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1138 //In order to avoid rewriting the same histograms
1139 Bool_t oldStatus = TH1::AddDirectoryStatus();
1140 TH1::AddDirectory(kFALSE);
1141
1142 fEMCALTimeRecalibrationFactors = new TObjArray(4);
1143 for (int i = 0; i < 4; i++)
3bfc4732 1144 fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1145 Form("hAllTimeAvBC%d",i),
c3d9f926 1146 48*24*12,0.,48*24*12) );
841dbf60 1147 //Init the histograms with 1
7f5392da 1148 for (Int_t bc = 0; bc < 4; bc++)
1149 {
c3d9f926 1150 for (Int_t i = 0; i < 48*24*12; i++)
841dbf60 1151 SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
3bfc4732 1152 }
841dbf60 1153
1154 fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1155 fEMCALTimeRecalibrationFactors->Compress();
1156
1157 //In order to avoid rewriting the same histograms
fbddd006 1158 TH1::AddDirectory(oldStatus);
3bfc4732 1159}
094786cc 1160
a520bcd0 1161//____________________________________________________
841dbf60 1162void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1163{
1164 //Init EMCAL bad channels map
1165 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1166 //In order to avoid rewriting the same histograms
1167 Bool_t oldStatus = TH1::AddDirectoryStatus();
608c80a3 1168 TH1::AddDirectory(kFALSE);
1169
c3d9f926 1170 fEMCALBadChannelMap = new TObjArray(12);
841dbf60 1171 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
7f5392da 1172 for (int i = 0; i < 12; i++)
1173 {
841dbf60 1174 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1175 }
1176
1177 fEMCALBadChannelMap->SetOwner(kTRUE);
1178 fEMCALBadChannelMap->Compress();
1179
1180 //In order to avoid rewriting the same histograms
fbddd006 1181 TH1::AddDirectory(oldStatus);
fd6df01c 1182}
1183
88b96ad8 1184//____________________________________________________________________________
1185void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1186 AliVCluster * cluster,
1187 AliVCaloCells * cells,
fb8eab96 1188 Int_t bc)
88b96ad8 1189{
841dbf60 1190 // Recalibrate the cluster energy and Time, considering the recalibration map
3bfc4732 1191 // and the energy of the cells and time that compose the cluster.
1192 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
fbddd006 1193
8517fbd8 1194 if (!cluster) {
2aeb4226 1195 AliInfo("Cluster pointer null!");
1196 return;
1197 }
1198
841dbf60 1199 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1200 UShort_t * index = cluster->GetCellsAbsId() ;
1201 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1202 Int_t ncells = cluster->GetNCells();
fbddd006 1203
841dbf60 1204 //Initialize some used variables
1205 Float_t energy = 0;
1206 Int_t absId =-1;
3bfc4732 1207 Int_t icol =-1, irow =-1, imod=1;
841dbf60 1208 Float_t factor = 1, frac = 0;
3bfc4732 1209 Int_t absIdMax = -1;
1210 Float_t emax = 0;
1211
841dbf60 1212 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
8517fbd8 1213 for (Int_t icell = 0; icell < ncells; icell++)
7f5392da 1214 {
841dbf60 1215 absId = index[icell];
1216 frac = fraction[icell];
8517fbd8 1217 if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
3bfc4732 1218
8517fbd8 1219 if (!fCellsRecalibrated && IsRecalibrationOn()) {
3bfc4732 1220 // Energy
1221 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1222 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
8517fbd8 1223 if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1224 continue;
fbddd006 1225 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
3bfc4732 1226 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1227
1228 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1229 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1230
1231 }
1232
841dbf60 1233 energy += cells->GetCellAmplitude(absId)*factor*frac;
3bfc4732 1234
8517fbd8 1235 if (emax < cells->GetCellAmplitude(absId)*factor*frac) {
3bfc4732 1236 emax = cells->GetCellAmplitude(absId)*factor*frac;
1237 absIdMax = absId;
1238 }
841dbf60 1239 }
fbddd006 1240
3a5708cd 1241 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
3bfc4732 1242
3a5708cd 1243 cluster->SetE(energy);
3bfc4732 1244
bc8330e9 1245 // Recalculate time of cluster
3a5708cd 1246 Double_t timeorg = cluster->GetTOF();
bc8330e9 1247
1248 Double_t time = cells->GetCellTime(absIdMax);
8517fbd8 1249 if (!fCellsRecalibrated && IsTimeRecalibrationOn())
3a5708cd 1250 RecalibrateCellTime(absIdMax,bc,time);
3bfc4732 1251
bc8330e9 1252 cluster->SetTOF(time);
3bfc4732 1253
bc8330e9 1254 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
3bfc4732 1255}
1256
a520bcd0 1257//_____________________________________________________________
1258void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells,
fb8eab96 1259 Int_t bc)
841dbf60 1260{
1261 // Recalibrate the cells time and energy, considering the recalibration map and the energy
3bfc4732 1262 // of the cells that compose the cluster.
1263 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1264
8517fbd8 1265 if (!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn())
1266 return;
3bfc4732 1267
8517fbd8 1268 if (!cells) {
3bfc4732 1269 AliInfo("Cells pointer null!");
1270 return;
1271 }
1272
fbddd006 1273 Short_t absId =-1;
a7e5a381 1274 Bool_t accept = kFALSE;
1275 Float_t ecell = 0;
1276 Double_t tcell = 0;
fbddd006 1277 Double_t ecellin = 0;
1278 Double_t tcellin = 0;
60d77596 1279 Int_t mclabel = -1;
fbddd006 1280 Double_t efrac = 0;
3bfc4732 1281
841dbf60 1282 Int_t nEMcell = cells->GetNumberOfCells() ;
7f5392da 1283 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1284 {
fbddd006 1285 cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
3bfc4732 1286
a7e5a381 1287 accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
8517fbd8 1288 if (!accept) {
a7e5a381 1289 ecell = 0;
fbddd006 1290 tcell = -1;
a7e5a381 1291 }
3bfc4732 1292
1293 //Set new values
fbddd006 1294 cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
3bfc4732 1295 }
841dbf60 1296
1297 fCellsRecalibrated = kTRUE;
094786cc 1298}
1299
88b96ad8 1300//_______________________________________________________________________________________________________
fb8eab96 1301void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime) const
7d692da6 1302{
841dbf60 1303 // Recalibrate time of cell with absID considering the recalibration map
3bfc4732 1304 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
7d692da6 1305
8517fbd8 1306 if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
7d692da6 1307 celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; ;
3bfc4732 1308 }
3bfc4732 1309}
1310
a520bcd0 1311//______________________________________________________________________________
1312void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1313 AliVCaloCells* cells,
1314 AliVCluster* clu)
d9b3567c 1315{
1316 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1317
8517fbd8 1318 if (!clu) {
2aeb4226 1319 AliInfo("Cluster pointer null!");
1320 return;
1321 }
1322
8517fbd8 1323 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1324 else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1325 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 1326}
1327
a520bcd0 1328//_____________________________________________________________________________________________
1329void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom,
1330 AliVCaloCells* cells,
1331 AliVCluster* clu)
094786cc 1332{
1333 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1334 // The algorithm is a copy of what is done in AliEMCALRecPoint
1335
1336 Double_t eCell = 0.;
1337 Float_t fraction = 1.;
1338 Float_t recalFactor = 1.;
1339
1340 Int_t absId = -1;
1341 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1342 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1343 Float_t weight = 0., totalWeight=0.;
1344 Float_t newPos[3] = {0,0,0};
1345 Double_t pLocal[3], pGlobal[3];
cb231979 1346 Bool_t shared = kFALSE;
1347
094786cc 1348 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
a840d589 1349 if (clEnergy <= 0)
1350 return;
cb231979 1351 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1352 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1353
83bfd77a 1354 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1355
7f5392da 1356 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1357 {
841dbf60 1358 absId = clu->GetCellAbsId(iDig);
1359 fraction = clu->GetCellAmplitudeFraction(iDig);
8517fbd8 1360 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1361
8517fbd8 1362 if (!fCellsRecalibrated) {
3bfc4732 1363 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
fbddd006 1364 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
8517fbd8 1365 if (IsRecalibrationOn()) {
3bfc4732 1366 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1367 }
094786cc 1368 }
3bfc4732 1369
094786cc 1370 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1371
1372 weight = GetCellWeight(eCell,clEnergy);
1373 totalWeight += weight;
3bfc4732 1374
094786cc 1375 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 1376 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 1377 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 1378 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1379
8517fbd8 1380 for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
094786cc 1381 }// cell loop
1382
8517fbd8 1383 if (totalWeight>0) {
1384 for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
094786cc 1385 }
1386
094786cc 1387 //Float_t pos[]={0,0,0};
1388 //clu->GetPosition(pos);
1389 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 1390 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 1391
8517fbd8 1392 if (iSupModMax > 1) { //sector 1
841dbf60 1393 newPos[0] +=fMisalTransShift[3];//-=3.093;
1394 newPos[1] +=fMisalTransShift[4];//+=6.82;
1395 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 1396 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
8517fbd8 1397 } else { //sector 0
841dbf60 1398 newPos[0] +=fMisalTransShift[0];//+=1.134;
1399 newPos[1] +=fMisalTransShift[1];//+=8.2;
1400 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 1401 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
841dbf60 1402 }
83bfd77a 1403 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1404
094786cc 1405 clu->SetPosition(newPos);
094786cc 1406}
1407
a520bcd0 1408//____________________________________________________________________________________________
1409void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom,
1410 AliVCaloCells* cells,
1411 AliVCluster* clu)
094786cc 1412{
1413 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1414 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1415
1416 Double_t eCell = 1.;
1417 Float_t fraction = 1.;
1418 Float_t recalFactor = 1.;
1419
1420 Int_t absId = -1;
d9b3567c 1421 Int_t iTower = -1;
094786cc 1422 Int_t iIphi = -1, iIeta = -1;
841dbf60 1423 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 1424 Int_t iphi = -1, ieta =-1;
cb231979 1425 Bool_t shared = kFALSE;
1426
d9b3567c 1427 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
a1c0f640 1428
a840d589 1429 if (clEnergy <= 0)
1430 return;
cb231979 1431 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1432 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1433
d9b3567c 1434 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 1435 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1436 Int_t startingSM = -1;
d9b3567c 1437
7f5392da 1438 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1439 {
094786cc 1440 absId = clu->GetCellAbsId(iDig);
1441 fraction = clu->GetCellAmplitudeFraction(iDig);
8517fbd8 1442 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1443
8517fbd8 1444 if (iDig==0) startingSM = iSupMod;
1445 else if (iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 1446
1447 eCell = cells->GetCellAmplitude(absId);
d9b3567c 1448
3bfc4732 1449 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
fbddd006 1450 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
3bfc4732 1451
7f5392da 1452 if (!fCellsRecalibrated)
1453 {
8517fbd8 1454 if (IsRecalibrationOn()) {
3bfc4732 1455 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3bfc4732 1456 }
094786cc 1457 }
3bfc4732 1458
094786cc 1459 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 1460
094786cc 1461 weight = GetCellWeight(eCell,clEnergy);
8517fbd8 1462 if (weight < 0) weight = 0;
d9b3567c 1463 totalWeight += weight;
1464 weightedCol += ieta*weight;
1465 weightedRow += iphi*weight;
1466
1467 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
841dbf60 1468 }// cell loop
094786cc 1469
d9b3567c 1470 Float_t xyzNew[]={0.,0.,0.};
8517fbd8 1471 if (areInSameSM == kTRUE) {
d9b3567c 1472 //printf("In Same SM\n");
1473 weightedCol = weightedCol/totalWeight;
1474 weightedRow = weightedRow/totalWeight;
094786cc 1475 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
7f5392da 1476 }
1477 else
1478 {
d9b3567c 1479 //printf("In Different SM\n");
094786cc 1480 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 1481 }
d9b3567c 1482
094786cc 1483 clu->SetPosition(xyzNew);
d9b3567c 1484}
1485
a520bcd0 1486//___________________________________________________________________________________________
1487void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
1488 AliVCaloCells* cells,
1489 AliVCluster * cluster)
841dbf60 1490{
cb231979 1491 //re-evaluate distance to bad channel with updated bad map
1492
8517fbd8 1493 if (!fRecalDistToBadChannels) return;
cb231979 1494
8517fbd8 1495 if (!cluster)
7f5392da 1496 {
2aeb4226 1497 AliInfo("Cluster pointer null!");
1498 return;
1499 }
1500
841dbf60 1501 //Get channels map of the supermodule where the cluster is.
fbddd006 1502 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
cb231979 1503 Bool_t shared = kFALSE;
1504 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1505 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1506
fbddd006 1507 Int_t dRrow, dRcol;
841dbf60 1508 Float_t minDist = 10000.;
1509 Float_t dist = 0.;
cb231979 1510
1511 //Loop on tower status map
8517fbd8 1512 for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
7f5392da 1513 {
8517fbd8 1514 for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
7f5392da 1515 {
841dbf60 1516 //Check if tower is bad.
8517fbd8 1517 if (hMap->GetBinContent(icol,irow)==0) continue;
cb231979 1518 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
6fe0e6d0 1519 // iSupMod,icol, irow, icolM,irowM);
cb231979 1520
1521 dRrow=TMath::Abs(irowM-irow);
1522 dRcol=TMath::Abs(icolM-icol);
1523 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
8517fbd8 1524 if (dist < minDist)
7f5392da 1525 {
cb231979 1526 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1527 minDist = dist;
1528 }
841dbf60 1529 }
1530 }
cb231979 1531
841dbf60 1532 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
7f5392da 1533 if (shared)
1534 {
841dbf60 1535 TH2D* hMap2 = 0;
1536 Int_t iSupMod2 = -1;
cb231979 1537
841dbf60 1538 //The only possible combinations are (0,1), (2,3) ... (8,9)
8517fbd8 1539 if (iSupMod%2) iSupMod2 = iSupMod-1;
1540 else iSupMod2 = iSupMod+1;
841dbf60 1541 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
cb231979 1542
841dbf60 1543 //Loop on tower status map of second super module
8517fbd8 1544 for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
7f5392da 1545 {
8517fbd8 1546 for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
7f5392da 1547 {
841dbf60 1548 //Check if tower is bad.
8517fbd8 1549 if (hMap2->GetBinContent(icol,irow)==0)
1550 continue;
841dbf60 1551 //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",
1552 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
cb231979 1553 dRrow=TMath::Abs(irow-irowM);
1554
8517fbd8 1555 if (iSupMod%2) {
841dbf60 1556 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
8517fbd8 1557 } else {
cb231979 1558 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
841dbf60 1559 }
cb231979 1560
841dbf60 1561 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
8517fbd8 1562 if (dist < minDist) minDist = dist;
841dbf60 1563 }
1564 }
1565 }// shared cluster in 2 SuperModules
78467229 1566
6fe0e6d0 1567 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1568 cluster->SetDistanceToBadChannel(minDist);
cb231979 1569}
1570
a520bcd0 1571//__________________________________________________________________
841dbf60 1572void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1573{
83bfd77a 1574 //re-evaluate identification parameters with bayesian
2aeb4226 1575
8517fbd8 1576 if (!cluster) {
2aeb4226 1577 AliInfo("Cluster pointer null!");
1578 return;
1579 }
1580
8517fbd8 1581 if (cluster->GetM02() != 0)
83bfd77a 1582 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1583
00a38d07 1584 Float_t pidlist[AliPID::kSPECIESCN+1];
8517fbd8 1585 for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
841dbf60 1586
83bfd77a 1587 cluster->SetPID(pidlist);
83bfd77a 1588}
1589
f0e9e976 1590//___________________________________________________________________________________________________________________
a520bcd0 1591void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1592 AliVCaloCells* cells,
f0e9e976 1593 AliVCluster * cluster,
1594 Float_t & l0, Float_t & l1,
1595 Float_t & disp, Float_t & dEta, Float_t & dPhi,
1596 Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
83bfd77a 1597{
1598 // Calculates new center of gravity in the local EMCAL-module coordinates
1599 // and tranfers into global ALICE coordinates
1600 // Calculates Dispersion and main axis
1601
8517fbd8 1602 if (!cluster) {
2aeb4226 1603 AliInfo("Cluster pointer null!");
1604 return;
1605 }
1606
83bfd77a 1607 Double_t eCell = 0.;
1608 Float_t fraction = 1.;
1609 Float_t recalFactor = 1.;
1610
f0e9e976 1611 Int_t iSupMod = -1;
1612 Int_t iTower = -1;
1613 Int_t iIphi = -1;
1614 Int_t iIeta = -1;
1615 Int_t iphi = -1;
1616 Int_t ieta = -1;
1617 Double_t etai = -1.;
1618 Double_t phii = -1.;
1619
1620 Int_t nstat = 0 ;
1621 Float_t wtot = 0.;
1622 Double_t w = 0.;
1623 Double_t etaMean = 0.;
1624 Double_t phiMean = 0.;
a1c0f640 1625
1626 //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added
1627 // and to check if the cluster is between 2 SM in eta
1628 Int_t iSM0 = -1;
1629 Bool_t shared = kFALSE;
1630 Float_t energy = 0;
1631
8517fbd8 1632 for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
a1c0f640 1633 {
1634 //Get from the absid the supermodule, tower and eta/phi numbers
1635 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1636 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1637
1638 //Check if there are cells of different SM
8517fbd8 1639 if (iDigit == 0 ) iSM0 = iSupMod;
1640 else if (iSupMod!= iSM0) shared = kTRUE;
a1c0f640 1641
1642 //Get the cell energy, if recalibration is on, apply factors
1643 fraction = cluster->GetCellAmplitudeFraction(iDigit);
8517fbd8 1644 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
83bfd77a 1645
8517fbd8 1646 if (IsRecalibrationOn()) {
a1c0f640 1647 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1648 }
1649
1650 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1651
1652 energy += eCell;
1653
1654 }//cell loop
1655
83bfd77a 1656 //Loop on cells
8517fbd8 1657 for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
7f5392da 1658 {
83bfd77a 1659 //Get from the absid the supermodule, tower and eta/phi numbers
1660 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1661 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1662
1663 //Get the cell energy, if recalibration is on, apply factors
1664 fraction = cluster->GetCellAmplitudeFraction(iDigit);
8517fbd8 1665 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1666
8517fbd8 1667 if (!fCellsRecalibrated) {
1668 if (IsRecalibrationOn()) {
3bfc4732 1669 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1670 }
83bfd77a 1671 }
3bfc4732 1672
83bfd77a 1673 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1674
a1c0f640 1675 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1676 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
8517fbd8 1677 if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
a1c0f640 1678
8517fbd8 1679 if (cluster->E() > 0 && eCell > 0) {
83bfd77a 1680 w = GetCellWeight(eCell,cluster->E());
1681
1682 etai=(Double_t)ieta;
fbddd006 1683 phii=(Double_t)iphi;
f0e9e976 1684
8517fbd8 1685 if (w > 0.0) {
83bfd77a 1686 wtot += w ;
fbddd006 1687 nstat++;
83bfd77a 1688 //Shower shape
f0e9e976 1689 sEta += w * etai * etai ;
1690 etaMean += w * etai ;
1691 sPhi += w * phii * phii ;
1692 phiMean += w * phii ;
1693 sEtaPhi += w * etai * phii ;
83bfd77a 1694 }
8517fbd8 1695 } else
83bfd77a 1696 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1697 }//cell loop
1698
fbddd006 1699 //Normalize to the weight
8517fbd8 1700 if (wtot > 0) {
f0e9e976 1701 etaMean /= wtot ;
1702 phiMean /= wtot ;
8517fbd8 1703 } else
83bfd77a 1704 AliError(Form("Wrong weight %f\n", wtot));
1705
fbddd006 1706 //Calculate dispersion
8517fbd8 1707 for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
7f5392da 1708 {
83bfd77a 1709 //Get from the absid the supermodule, tower and eta/phi numbers
1710 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1711 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1712
1713 //Get the cell energy, if recalibration is on, apply factors
1714 fraction = cluster->GetCellAmplitudeFraction(iDigit);
8517fbd8 1715 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1716 if (IsRecalibrationOn()) {
83bfd77a 1717 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1718 }
1719 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1720
a1c0f640 1721 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1722 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
8517fbd8 1723 if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
a1c0f640 1724
8517fbd8 1725 if (cluster->E() > 0 && eCell > 0) {
83bfd77a 1726 w = GetCellWeight(eCell,cluster->E());
1727
1728 etai=(Double_t)ieta;
fbddd006 1729 phii=(Double_t)iphi;
8517fbd8 1730 if (w > 0.0) {
f0e9e976 1731 disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
1732 dEta += w * (etai-etaMean)*(etai-etaMean) ;
1733 dPhi += w * (phii-phiMean)*(phii-phiMean) ;
1734 }
8517fbd8 1735 } else
83bfd77a 1736 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1737 }// cell loop
1738
1739 //Normalize to the weigth and set shower shape parameters
8517fbd8 1740 if (wtot > 0 && nstat > 1) {
f0e9e976 1741 disp /= wtot ;
1742 dEta /= wtot ;
1743 dPhi /= wtot ;
1744 sEta /= wtot ;
1745 sPhi /= wtot ;
1746 sEtaPhi /= wtot ;
1747
1748 sEta -= etaMean * etaMean ;
1749 sPhi -= phiMean * phiMean ;
1750 sEtaPhi -= etaMean * phiMean ;
1751
1752 l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1753 l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
8517fbd8 1754 } else {
f0e9e976 1755 l0 = 0. ;
1756 l1 = 0. ;
1757 dEta = 0. ; dPhi = 0. ; disp = 0. ;
1758 sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
fbddd006 1759 }
83bfd77a 1760}
1761
f0e9e976 1762//____________________________________________________________________________________________
1763void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1764 AliVCaloCells* cells,
1765 AliVCluster * cluster)
1766{
1767 // Calculates new center of gravity in the local EMCAL-module coordinates
1768 // and tranfers into global ALICE coordinates
1769 // Calculates Dispersion and main axis and puts them into the cluster
1770
1771 Float_t l0 = 0., l1 = 0.;
1772 Float_t disp = 0., dEta = 0., dPhi = 0.;
1773 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1774
1775 AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
1776 dEta, dPhi, sEta, sPhi, sEtaPhi);
1777
1778 cluster->SetM02(l0);
1779 cluster->SetM20(l1);
8517fbd8 1780 if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
f0e9e976 1781
1782}
1783
b540d03f 1784//____________________________________________________________________________
a520bcd0 1785void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
1786 TObjArray * clusterArr,
1787 const AliEMCALGeometry *geom)
bd8c7aef 1788{
1789 //This function should be called before the cluster loop
1790 //Before call this function, please recalculate the cluster positions
1791 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1792 //Store matched cluster indexes and residuals
61160f1f 1793
88b96ad8 1794 fMatchedTrackIndex ->Reset();
bd8c7aef 1795 fMatchedClusterIndex->Reset();
fa4287a2 1796 fResidualPhi->Reset();
1797 fResidualEta->Reset();
bd8c7aef 1798
841dbf60 1799 fMatchedTrackIndex ->Set(1000);
1800 fMatchedClusterIndex->Set(1000);
1801 fResidualPhi->Set(1000);
1802 fResidualEta->Set(1000);
bd8c7aef 1803
1c7a2bf4 1804 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1805 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
57131575 1806
88b96ad8 1807 // init the magnetic field if not already on
8517fbd8 1808 if (!TGeoGlobalMagField::Instance()->GetField()) {
092f179b 1809 if (!event->InitMagneticField())
88b96ad8 1810 {
1811 AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1812 }
88b96ad8 1813 } // Init mag field
1814
42ceff04 1815 if (esdevent) {
1816 UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
1817 UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
1818 Bool_t desc1 = (mask1 >> 3) & 0x1;
1819 Bool_t desc2 = (mask2 >> 3) & 0x1;
1820 if (desc1==0 || desc2==0) {
30e29b2a 1821// AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
1822// mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
1823// mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
42ceff04 1824 fITSTrackSA=kTRUE;
1825 }
1826 }
1827
8fc351e3 1828 TObjArray *clusterArray = 0x0;
8517fbd8 1829 if (!clusterArr) {
092f179b 1830 clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
8517fbd8 1831 for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
092f179b
CL		 1832 {
1833 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8517fbd8 1834 if (geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
1835 continue;
092f179b 1836 clusterArray->AddAt(cluster,icl);
8fc351e3 1837 }
092f179b 1838 }
61160f1f 1839
bd8c7aef 1840 Int_t matched=0;
bb6f5f0b 1841 Double_t cv[21];
1842 for (Int_t i=0; i<21;i++) cv[i]=0;
8517fbd8 1843 for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
bd8c7aef 1844 {
456126ad 1845 AliExternalTrackParam *trackParam = 0;
61160f1f 1846
bb6f5f0b 1847 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
8fc351e3 1848 AliESDtrack *esdTrack = 0;
1849 AliAODTrack *aodTrack = 0;
8517fbd8 1850 if (esdevent) {
8fc351e3 1851 esdTrack = esdevent->GetTrack(itr);
8517fbd8 1852 if (!esdTrack) continue;
1853 if (!IsAccepted(esdTrack)) continue;
1854 if (esdTrack->Pt()<fCutMinTrackPt) continue;
8fc351e3 1855 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
8517fbd8 1856 if (TMath::Abs(esdTrack->Eta())>0.9 || phi <= 10 || phi >= 250 ) continue;
1857 if (!fITSTrackSA)
42ceff04 1858 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
1859 else
1860 trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
61160f1f 1861 }
bb6f5f0b 1862
1863 //If the input event is AOD, the starting point for extrapolation is at vertex
8fc351e3 1864 //AOD tracks are selected according to its filterbit.
8517fbd8 1865 else if (aodevent) {
f15c1f69 1866 aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
1867 if(!aodTrack) AliFatal("Not a standard AOD");
8517fbd8 1868 if (!aodTrack) continue;
a6a1e3ab 1869
8517fbd8 1870 if (fAODTPCOnlyTracks) { // Match with TPC only tracks, default from May 2013, before filter bit 32
a6a1e3ab 1871 //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
c860bb71 1872 if (!aodTrack->IsTPCConstrained()) continue ;
8517fbd8 1873 } else if (fAODHybridTracks) { // Match with hybrid tracks
a6a1e3ab 1874 //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
8517fbd8 1875 if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
1876 } else { // Match with tracks on a mask
a6a1e3ab 1877 //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
8517fbd8 1878 if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
a6a1e3ab 1879 }
1af378e6 1880
8517fbd8 1881 if (aodTrack->Pt()<fCutMinTrackPt) continue;
1af378e6 1882
8fc351e3 1883 Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
8517fbd8 1884 if (TMath::Abs(aodTrack->Eta())>0.9 || phi <= 10 || phi >= 250 )
1885 continue;
61160f1f 1886 Double_t pos[3],mom[3];
1887 aodTrack->GetXYZ(pos);
1888 aodTrack->GetPxPyPz(mom);
1889 AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
1af378e6 1890
61160f1f 1891 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1892 }
bd8c7aef 1893
bb6f5f0b 1894 //Return if the input data is not "AOD" or "ESD"
8517fbd8 1895 else {
61160f1f 1896 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
8517fbd8 1897 if (clusterArray) {
092f179b
CL		 1898 clusterArray->Clear();
1899 delete clusterArray;
1900 }
61160f1f 1901 return;
1902 }
1903
8517fbd8 1904 if (!trackParam) continue;
8fc351e3 1905
1906 //Extrapolate the track to EMCal surface
1907 AliExternalTrackParam emcalParam(*trackParam);
a29b2a8a 1908 Float_t eta, phi, pt;
8517fbd8 1909 if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) {
1910 if (aodevent && trackParam) delete trackParam;
1911 if (fITSTrackSA && trackParam) delete trackParam;
092f179b
CL		 1912 continue;
1913 }
8fc351e3 1914
8517fbd8 1915 if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) {
1916 if (aodevent && trackParam) delete trackParam;
1917 if (fITSTrackSA && trackParam) delete trackParam;
092f179b
CL		 1918 continue;
1919 }
8fc351e3 1920
1921 //Find matched clusters
bd8c7aef 1922 Int_t index = -1;
8fc351e3 1923 Float_t dEta = -999, dPhi = -999;
8517fbd8 1924 if (!clusterArr) {
092f179b 1925 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
8517fbd8 1926 } else {
092f179b
CL		 1927 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1928 }
61160f1f 1929
8517fbd8 1930 if (index>-1) {
8fc351e3 1931 fMatchedTrackIndex ->AddAt(itr,matched);
1932 fMatchedClusterIndex ->AddAt(index,matched);
1933 fResidualEta ->AddAt(dEta,matched);
1934 fResidualPhi ->AddAt(dPhi,matched);
bd8c7aef 1935 matched++;
1936 }
8517fbd8 1937 if (aodevent && trackParam) delete trackParam;
1938 if (fITSTrackSA && trackParam) delete trackParam;
bd8c7aef 1939 }//track loop
8fc351e3 1940
8517fbd8 1941 if (clusterArray) {
092f179b
CL		 1942 clusterArray->Clear();
1943 delete clusterArray;
1944 }
b540d03f 1945
1946 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1947
8fc351e3 1948 fMatchedTrackIndex ->Set(matched);
1949 fMatchedClusterIndex ->Set(matched);
1950 fResidualPhi ->Set(matched);
1951 fResidualEta ->Set(matched);
bd8c7aef 1952}
1953
b540d03f 1954//________________________________________________________________________________
a520bcd0 1955Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
1956 const AliVEvent *event,
1957 const AliEMCALGeometry *geom,
1958 Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1959{
1960 //
1961 // This function returns the index of matched cluster to input track
fa4287a2 1962 // Returns -1 if no match is found
bb6f5f0b 1963 Int_t index = -1;
8fc351e3 1964 Double_t phiV = track->Phi()*TMath::RadToDeg();
8517fbd8 1965 if (TMath::Abs(track->Eta())>0.9 || phiV <= 10 || phiV >= 250 ) return index;
42ceff04 1966 AliExternalTrackParam *trackParam = 0;
8517fbd8 1967 if (!fITSTrackSA)
42ceff04 1968 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
1969 else
1970 trackParam = new AliExternalTrackParam(*track);
1971
8517fbd8 1972 if (!trackParam) return index;
8fc351e3 1973 AliExternalTrackParam emcalParam(*trackParam);
a29b2a8a 1974 Float_t eta, phi, pt;
8fc351e3 1975
8517fbd8 1976 if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) {
1977 if (fITSTrackSA) delete trackParam;
1978 return index;
30e29b2a 1979 }
8517fbd8 1980 if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) {
1981 if (fITSTrackSA) delete trackParam;
1982 return index;
30e29b2a 1983 }
1984
8fc351e3 1985 TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
1986
8517fbd8 1987 for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 1988 {
1989 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8517fbd8 1990 if (geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
8fc351e3 1991 clusterArr->AddAt(cluster,icl);
1992 }
1993
1994 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1995 clusterArr->Clear();
1996 delete clusterArr;
8517fbd8 1997 if (fITSTrackSA) delete trackParam;
30e29b2a 1998
8fc351e3 1999 return index;
2000}
2001
7f5392da 2002//_______________________________________________________________________________________________
2003Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
a520bcd0 2004 AliExternalTrackParam *trkParam,
7f5392da 2005 const TObjArray * clusterArr,
a520bcd0 2006 Float_t &dEta, Float_t &dPhi)
8fc351e3 2007{
7f5392da 2008 // Find matched cluster in array
2009
8fc351e3 2010 dEta=-999, dPhi=-999;
2011 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2012 Int_t index = -1;
ee602376 2013 Float_t tmpEta=-999, tmpPhi=-999;
8fc351e3 2014
2015 Double_t exPos[3] = {0.,0.,0.};
8517fbd8 2016 if (!emcalParam->GetXYZ(exPos)) return index;
8fc351e3 2017
2018 Float_t clsPos[3] = {0.,0.,0.};
8517fbd8 2019 for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
092f179b
CL		 2020 {
2021 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
8517fbd8 2022 if (!cluster || !cluster->IsEMCAL()) continue;
092f179b
CL		 2023 cluster->GetPosition(clsPos);
2024 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));
8517fbd8 2025 if (dR > fClusterWindow) continue;
092f179b
CL		 2026
2027 AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
8517fbd8 2028 if (!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2029 if (fCutEtaPhiSum) {
092f179b 2030 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
8517fbd8 2031 if (tmpR<dRMax) {
fbddd006 2032 dRMax=tmpR;
2033 dEtaMax=tmpEta;
2034 dPhiMax=tmpPhi;
2035 index=icl;
2036 }
8517fbd8 2037 } else if (fCutEtaPhiSeparate) {
2038 if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax)) {
fbddd006 2039 dEtaMax = tmpEta;
2040 dPhiMax = tmpPhi;
2041 index=icl;
2042 }
8517fbd8 2043 } else {
092f179b
CL		 2044 printf("Error: please specify your cut criteria\n");
2045 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2046 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2047 return index;
2048 }
2049 }
8fc351e3 2050
2051 dEta=dEtaMax;
2052 dPhi=dPhiMax;
2053
bb6f5f0b 2054 return index;
2055}
2056
40891fa2
CL		 2057//------------------------------------------------------------------------------------
2058Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliVTrack *track,
18a566cc 2059 Double_t emcalR, Double_t mass,
2060 Double_t step, Double_t minpt,
2061 Bool_t useMassForTracking)
40891fa2 2062{
8517fbd8 2063 // Extrapolate track to EMCAL surface
40891fa2
CL		 2064
2065 track->SetTrackPhiEtaPtOnEMCal(-999, -999, -999);
2066
6e0ae02d 2067 if (track->Pt()<minpt)
40891fa2
CL		 2068 return kFALSE;
2069
2070 Double_t phi = track->Phi()*TMath::RadToDeg();
2071 if (TMath::Abs(track->Eta())>0.9 || phi <= 10 || phi >= 250)
2072 return kFALSE;
2073
2074 AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(track);
2075 AliAODTrack *aodt = 0;
2076 if (!esdt) {
2077 aodt = dynamic_cast<AliAODTrack*>(track);
2078 if (!aodt)
2079 return kFALSE;
2080 }
2081
26664f75 2082 // Select the mass hypothesis
2083 if ( mass < 0 )
2084 {
40891fa2 2085 Bool_t onlyTPC = kFALSE;
26664f75 2086 if ( mass == -99 ) onlyTPC=kTRUE;
2087
2088 if (esdt)
2089 {
18a566cc 2090 if ( useMassForTracking ) mass = esdt->GetMassForTracking();
2091 else mass = esdt->GetMass(onlyTPC);
2092 }
2093 else
26664f75 2094 {
2095 if ( useMassForTracking ) mass = aodt->GetMassForTracking();
2096 else mass = aodt->M();
2097 }
40891fa2
CL		 2098 }
2099
2100 AliExternalTrackParam *trackParam = 0;
2101 if (esdt) {
0e1c2338 2102 const AliExternalTrackParam *in = esdt->GetInnerParam();
7af66d73
MP		 2103 if (!in)
2104 return kFALSE;
2105 trackParam = new AliExternalTrackParam(*in);
40891fa2
CL		 2106 } else {
2107 Double_t xyz[3] = {0}, pxpypz[3] = {0}, cv[21] = {0};
2108 aodt->PxPyPz(pxpypz);
2109 aodt->XvYvZv(xyz);
2110 aodt->GetCovarianceXYZPxPyPz(cv);
2111 trackParam = new AliExternalTrackParam(xyz,pxpypz,cv,aodt->Charge());
2112 }
2113 if (!trackParam)
2114 return kFALSE;
2115
2116 Float_t etaout=-999, phiout=-999, ptout=-999;
2117 Bool_t ret = ExtrapolateTrackToEMCalSurface(trackParam,
2118 emcalR,
2119 mass,
2120 step,
2121 etaout,
2122 phiout,
2123 ptout);
3ad1e868 2124 delete trackParam;
40891fa2
CL		 2125 if (!ret)
2126 return kFALSE;
40891fa2
CL		 2127 if (TMath::Abs(etaout)>0.75 || (phiout<70*TMath::DegToRad()) || (phiout>190*TMath::DegToRad()))
2128 return kFALSE;
40891fa2
CL		 2129 track->SetTrackPhiEtaPtOnEMCal(phiout, etaout, ptout);
2130 return kTRUE;
2131}
2132
2133
88b96ad8 2134//------------------------------------------------------------------------------------
2135Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
fb8eab96 2136 Double_t emcalR,
2137 Double_t mass,
2138 Double_t step,
88b96ad8 2139 Float_t &eta,
a29b2a8a 2140 Float_t &phi,
2141 Float_t &pt)
ee602376 2142{
88b96ad8 2143 //Extrapolate track to EMCAL surface
2144
a29b2a8a 2145 eta = -999, phi = -999, pt = -999;
8517fbd8 2146 if (!trkParam) return kFALSE;
2147 if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
ee602376 2148 Double_t trkPos[3] = {0.,0.,0.};
8517fbd8 2149 if (!trkParam->GetXYZ(trkPos)) return kFALSE;
ee602376 2150 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2151 eta = trkPosVec.Eta();
2152 phi = trkPosVec.Phi();
a29b2a8a 2153 pt = trkParam->Pt();
8517fbd8 2154 if (phi<0)
ee602376 2155 phi += 2*TMath::Pi();
2156
2157 return kTRUE;
2158}
2159
88b96ad8 2160//-----------------------------------------------------------------------------------
2161Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2162 const Float_t *clsPos,
2163 Double_t mass,
2164 Double_t step,
2165 Float_t &tmpEta,
2166 Float_t &tmpPhi)
ee602376 2167{
2168 //
2169 //Return the residual by extrapolating a track param to a global position
2170 //
2171 tmpEta = -999;
2172 tmpPhi = -999;
8517fbd8 2173 if (!trkParam) return kFALSE;
ee602376 2174 Double_t trkPos[3] = {0.,0.,0.};
2175 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2176 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2177 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
8517fbd8 2178 if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
2179 if (!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
ee602376 2180
2181 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2182 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2183
2184 // track cluster matching
2185 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
2186 tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2187
2188 return kTRUE;
2189}
2190
88b96ad8 2191//----------------------------------------------------------------------------------
2192Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2193 const AliVCluster *cluster,
fb8eab96 2194 Double_t mass,
2195 Double_t step,
88b96ad8 2196 Float_t &tmpEta,
2197 Float_t &tmpPhi)
ee602376 2198{
2199 //
2200 //Return the residual by extrapolating a track param to a cluster
2201 //
2202 tmpEta = -999;
2203 tmpPhi = -999;
8517fbd8 2204 if (!cluster || !trkParam)
2205 return kFALSE;
ee602376 2206
2207 Float_t clsPos[3] = {0.,0.,0.};
2208 cluster->GetPosition(clsPos);
2209
2210 return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2211}
2212
88b96ad8 2213//---------------------------------------------------------------------------------
2214Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2215 const AliVCluster *cluster,
88b96ad8 2216 Float_t &tmpEta,
2217 Float_t &tmpPhi)
bb6f5f0b 2218{
2219 //
ee602376 2220 //Return the residual by extrapolating a track param to a clusterfStepCluster
bb6f5f0b 2221 //
8fc351e3 2222
ee602376 2223 return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
bb6f5f0b 2224}
2225
a520bcd0 2226//_______________________________________________________________________
fb8eab96 2227void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex,
a520bcd0 2228 Float_t &dEta, Float_t &dPhi)
bd8c7aef 2229{
bb6f5f0b 2230 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 2231 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 2232 //Works with ESDs and AODs
bd8c7aef 2233
8517fbd8 2234 if (FindMatchedPosForCluster(clsIndex) >= 999) {
bd8c7aef 2235 AliDebug(2,"No matched tracks found!\n");
fa4287a2 2236 dEta=999.;
2237 dPhi=999.;
bd8c7aef 2238 return;
2239 }
fa4287a2 2240 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2241 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 2242}
841dbf60 2243
88b96ad8 2244//______________________________________________________________________________________________
fa4287a2 2245void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 2246{
2247 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 2248 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 2249 //Works with ESDs and AODs
2250
8517fbd8 2251 if (FindMatchedPosForTrack(trkIndex) >= 999) {
bb6f5f0b 2252 AliDebug(2,"No matched cluster found!\n");
fa4287a2 2253 dEta=999.;
2254 dPhi=999.;
bb6f5f0b 2255 return;
2256 }
fa4287a2 2257 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2258 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 2259}
2260
2261//__________________________________________________________
2262Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2263{
2264 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2265 //Get the index of matched track to this cluster
2266 //Works with ESDs and AODs
2267
8517fbd8 2268 if (IsClusterMatched(clsIndex))
bb6f5f0b 2269 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2270 else
2271 return -1;
bd8c7aef 2272}
2273
b540d03f 2274//__________________________________________________________
bb6f5f0b 2275Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 2276{
bb6f5f0b 2277 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2278 //Get the index of matched cluster to this track
2279 //Works with ESDs and AODs
b540d03f 2280
8517fbd8 2281 if (IsTrackMatched(trkIndex))
bb6f5f0b 2282 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 2283 else
2284 return -1;
2285}
2286
7f5392da 2287//______________________________________________________________
7cdec71f 2288Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 2289{
2290 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2291 //Returns if the cluster has a match
8517fbd8 2292 if (FindMatchedPosForCluster(clsIndex) < 999)
bb6f5f0b 2293 return kTRUE;
2294 else
2295 return kFALSE;
2296}
b540d03f 2297
7f5392da 2298//____________________________________________________________
7cdec71f 2299Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 2300{
bb6f5f0b 2301 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2302 //Returns if the track has a match
8517fbd8 2303 if (FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 2304 return kTRUE;
bd8c7aef 2305 else
2306 return kFALSE;
2307}
bb6f5f0b 2308
7f5392da 2309//______________________________________________________________________
bb6f5f0b 2310UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 2311{
bb6f5f0b 2312 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 2313 //Returns the position of the match in the fMatchedClusterIndex array
2314 Float_t tmpR = fCutR;
81efb149 2315 UInt_t pos = 999;
b540d03f 2316
8517fbd8 2317 for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
7f5392da 2318 {
8517fbd8 2319 if (fMatchedClusterIndex->At(i)==clsIndex) {
841dbf60 2320 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
8517fbd8 2321 if (r<tmpR) {
841dbf60 2322 pos=i;
2323 tmpR=r;
7f5392da 2324 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2325 fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2326 }
841dbf60 2327 }
bb6f5f0b 2328 }
2329 return pos;
2330}
2331
7f5392da 2332//____________________________________________________________________
bb6f5f0b 2333UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2334{
2335 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2336 //Returns the position of the match in the fMatchedTrackIndex array
2337 Float_t tmpR = fCutR;
2338 UInt_t pos = 999;
2339
8517fbd8 2340 for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
7f5392da 2341 {
8517fbd8 2342 if (fMatchedTrackIndex->At(i)==trkIndex) {
841dbf60 2343 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
8517fbd8 2344 if (r<tmpR) {
841dbf60 2345 pos=i;
2346 tmpR=r;
7f5392da 2347 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2348 fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2349 }
841dbf60 2350 }
b540d03f 2351 }
bd8c7aef 2352 return pos;
2353}
2354
a520bcd0 2355//__________________________________________________________________________
2356Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
2357 const AliEMCALGeometry *geom,
fb8eab96 2358 AliVCaloCells* cells, Int_t bc)
b5078f5d 2359{
2360 // check if the cluster survives some quality cut
2361 //
2362 //
2363 Bool_t isGood=kTRUE;
7f5392da 2364
8517fbd8 2365 if (!cluster || !cluster->IsEMCAL()) return kFALSE;
2366 if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
2367 if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
2368 if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
b5078f5d 2369
2370 return isGood;
2371}
2372
b540d03f 2373//__________________________________________________________
bd8c7aef 2374Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2375{
2376 // Given a esd track, return whether the track survive all the cuts
2377
2378 // The different quality parameter are first
2379 // retrieved from the track. then it is found out what cuts the
2380 // track did not survive and finally the cuts are imposed.
2381
2382 UInt_t status = esdTrack->GetStatus();
2383
2384 Int_t nClustersITS = esdTrack->GetITSclusters(0);
2385 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2386
2387 Float_t chi2PerClusterITS = -1;
2388 Float_t chi2PerClusterTPC = -1;
2389 if (nClustersITS!=0)
2390 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2391 if (nClustersTPC!=0)
2392 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 2393
2394
2395 //DCA cuts
8517fbd8 2396 if (fTrackCutsType==kGlobalCut) {
2397 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2398 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2399 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2400 }
82d09e74 2401
bd8c7aef 2402 Float_t b[2];
2403 Float_t bCov[3];
2404 esdTrack->GetImpactParameters(b,bCov);
8517fbd8 2405 if (bCov[0]<=0 || bCov[2]<=0) {
bd8c7aef 2406 AliDebug(1, "Estimated b resolution lower or equal zero!");
2407 bCov[0]=0; bCov[2]=0;
2408 }
2409
2410 Float_t dcaToVertexXY = b[0];
2411 Float_t dcaToVertexZ = b[1];
2412 Float_t dcaToVertex = -1;
2413
2414 if (fCutDCAToVertex2D)
2415 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2416 else
2417 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2418
2419 // cut the track?
2420
2421 Bool_t cuts[kNCuts];
2422 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2423
2424 // track quality cuts
2425 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2426 cuts[0]=kTRUE;
2427 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2428 cuts[1]=kTRUE;
2429 if (nClustersTPC<fCutMinNClusterTPC)
2430 cuts[2]=kTRUE;
2431 if (nClustersITS<fCutMinNClusterITS)
2432 cuts[3]=kTRUE;
2433 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
2434 cuts[4]=kTRUE;
2435 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
2436 cuts[5]=kTRUE;
2437 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2438 cuts[6]=kTRUE;
2439 if (fCutDCAToVertex2D && dcaToVertex > 1)
2440 cuts[7] = kTRUE;
2441 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2442 cuts[8] = kTRUE;
2443 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2444 cuts[9] = kTRUE;
2445
8517fbd8 2446 if (fTrackCutsType==kGlobalCut) {
2447 //Require at least one SPD point + anything else in ITS
2448 if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
2449 cuts[10] = kTRUE;
2450 }
82d09e74 2451
8517fbd8 2452 // ITS
2453 if (fCutRequireITSStandAlone || fCutRequireITSpureSA) {
2454 if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)) {
42ceff04 2455 // TPC tracks
2456 cuts[11] = kTRUE;
8517fbd8 2457 } else {
42ceff04 2458 // ITS standalone tracks
8517fbd8 2459 if (fCutRequireITSStandAlone && !fCutRequireITSpureSA) {
2460 if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
2461 } else if (fCutRequireITSpureSA) {
2462 if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
42ceff04 2463 }
2464 }
2465 }
2466
bd8c7aef 2467 Bool_t cut=kFALSE;
827f9f23 2468 for (Int_t i=0; i<kNCuts; i++)
7f5392da 2469 if (cuts[i]) { cut = kTRUE ; }
bd8c7aef 2470
2471 // cut the track
2472 if (cut)
2473 return kFALSE;
2474 else
2475 return kTRUE;
2476}
827f9f23 2477
88b96ad8 2478//_____________________________________
bd8c7aef 2479void AliEMCALRecoUtils::InitTrackCuts()
2480{
2481 //Intilize the track cut criteria
5f7714ad 2482 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
bd8c7aef 2483 //Also you can customize the cuts using the setters
82d09e74 2484
5f7714ad 2485 switch (fTrackCutsType)
88b96ad8 2486 {
5f7714ad 2487 case kTPCOnlyCut:
88b96ad8 2488 {
2489 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2490 //TPC
2491 SetMinNClustersTPC(70);
2492 SetMaxChi2PerClusterTPC(4);
2493 SetAcceptKinkDaughters(kFALSE);
2494 SetRequireTPCRefit(kFALSE);
2495
2496 //ITS
2497 SetRequireITSRefit(kFALSE);
2498 SetMaxDCAToVertexZ(3.2);
2499 SetMaxDCAToVertexXY(2.4);
2500 SetDCAToVertex2D(kTRUE);
2501
2502 break;
2503 }
2504
5f7714ad 2505 case kGlobalCut:
88b96ad8 2506 {
2507 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2508 //TPC
2509 SetMinNClustersTPC(70);
2510 SetMaxChi2PerClusterTPC(4);
2511 SetAcceptKinkDaughters(kFALSE);
2512 SetRequireTPCRefit(kTRUE);
2513
2514 //ITS
2515 SetRequireITSRefit(kTRUE);
2516 SetMaxDCAToVertexZ(2);
2517 SetMaxDCAToVertexXY();
2518 SetDCAToVertex2D(kFALSE);
2519
2520 break;
2521 }
2522
0e7de35b 2523 case kLooseCut:
88b96ad8 2524 {
2525 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2526 SetMinNClustersTPC(50);
2527 SetAcceptKinkDaughters(kTRUE);
2528
2529 break;
5f7714ad 2530 }
42ceff04 2531
2532 case kITSStandAlone:
2533 {
2534 AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
2535 SetRequireITSRefit(kTRUE);
2536 SetRequireITSStandAlone(kTRUE);
2537 SetITSTrackSA(kTRUE);
2538 break;
2539 }
2540
88b96ad8 2541 }
bd8c7aef 2542}
83bfd77a 2543
57131575 2544
88b96ad8 2545//________________________________________________________________________
dda65b42 2546void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
57131575 2547{
2548 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
dda65b42 2549
57131575 2550 Int_t nTracks = event->GetNumberOfTracks();
7f5392da 2551 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
2552 {
dda65b42 2553 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
7f5392da 2554 if (!track)
2555 {
57131575 2556 AliWarning(Form("Could not receive track %d", iTrack));
2557 continue;
2558 }
7f5392da 2559
fbddd006 2560 Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
57131575 2561 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
dda65b42 2562 /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2563 AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2564 if (esdtrack) {
8517fbd8 2565 if (matchClusIndex != -1)
dda65b42 2566 esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2567 else
2568 esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2569 } else {
2570 AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
8517fbd8 2571 if (matchClusIndex != -1)
dda65b42 2572 aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2573 else
2574 aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2575 }
57131575 2576 }
fbddd006 2577 AliDebug(2,"Track matched to closest cluster");
57131575 2578}
2579
88b96ad8 2580//_________________________________________________________________________
dda65b42 2581void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
57131575 2582{
2583 // Checks if EMC clusters are matched to ESD track.
2584 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2585
7f5392da 2586 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
2587 {
dda65b42 2588 AliVCluster *cluster = event->GetCaloCluster(iClus);
57131575 2589 if (!cluster->IsEMCAL())
2590 continue;
2591
2592 Int_t nTracks = event->GetNumberOfTracks();
2593 TArrayI arrayTrackMatched(nTracks);
2594
2595 // Get the closest track matched to the cluster
2596 Int_t nMatched = 0;
2597 Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
7f5392da 2598 if (matchTrackIndex != -1)
2599 {
57131575 2600 arrayTrackMatched[nMatched] = matchTrackIndex;
2601 nMatched++;
2602 }
2603
2604 // Get all other tracks matched to the cluster
8517fbd8 2605 for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
7f5392da 2606 {
dda65b42 2607 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
0ddff478 2608
2609 if( !track ) continue;
2610
2611 if ( iTrk == matchTrackIndex ) continue;
2612
2613 if ( track->GetEMCALcluster() == iClus )
2614 {
57131575 2615 arrayTrackMatched[nMatched] = iTrk;
2616 ++nMatched;
2617 }
2618 }
2619
2620 //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2621
2622 arrayTrackMatched.Set(nMatched);
dda65b42 2623 AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2624 if (esdcluster)
2625 esdcluster->AddTracksMatched(arrayTrackMatched);
2626 else if (nMatched>0) {
2627 AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2628 if (aodcluster)
2629 aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2630 }
57131575 2631
2632 Float_t eta= -999, phi = -999;
2633 if (matchTrackIndex != -1)
2634 GetMatchedResiduals(iClus, eta, phi);
2635 cluster->SetTrackDistance(phi, eta);
2636 }
2637
fbddd006 2638 AliDebug(2,"Cluster matched to tracks");
57131575 2639}
2640
b540d03f 2641//___________________________________________________
d9b3567c 2642void AliEMCALRecoUtils::Print(const Option_t *) const
2643{
2644 // Print Parameters
2645
2646 printf("AliEMCALRecoUtils Settings: \n");
2647 printf("Misalignment shifts\n");
8517fbd8 2648 for (Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
2a71e873 2649 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2650 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 2651 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
8517fbd8 2652 for (Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 2653
2654 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 2655
fa4287a2 2656 printf("Matching criteria: ");
8517fbd8 2657 if (fCutEtaPhiSum) {
2658 printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
2659 } else if (fCutEtaPhiSeparate) {
2660 printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
2661 } else {
2662 printf("Error\n");
2663 printf("please specify your cut criteria\n");
2664 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2665 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2666 }
fa4287a2 2667
8fc351e3 2668 printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
2669 printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
bd8c7aef 2670
2671 printf("Track cuts: \n");
fa4287a2 2672 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
a6a1e3ab 2673 printf("AOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
bd8c7aef 2674 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2675 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2676 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2677 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2678 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
d9b3567c 2679}