]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALRecoUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Run initialization fixed
[u/mrichter/AliRoot.git] / EMCAL / 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
d9b3567c 730 case kNoCorrection:
731 AliDebug(2,"No correction on the energy\n");
732 break;
733
734 }
57131575 735
d9b3567c 736 return energy;
d9b3567c 737}
841dbf60 738
7e0ecb89 739//__________________________________________________
740void AliEMCALRecoUtils::InitNonLinearityParam()
741{
841dbf60 742 //Initialising Non Linearity Parameters
fbddd006 743
8517fbd8 744 if (fNonLinearityFunction == kPi0MC) {
841dbf60 745 fNonLinearityParams[0] = 1.014;
746 fNonLinearityParams[1] = -0.03329;
747 fNonLinearityParams[2] = -0.3853;
748 fNonLinearityParams[3] = 0.5423;
749 fNonLinearityParams[4] = -0.4335;
750 }
751
8517fbd8 752 if (fNonLinearityFunction == kPi0MCv2) {
e33e0c4a 753 fNonLinearityParams[0] = 3.11111e-02;
754 fNonLinearityParams[1] =-5.71666e-02;
755 fNonLinearityParams[2] = 5.67995e-01;
756 }
757
8517fbd8 758 if (fNonLinearityFunction == kPi0MCv3) {
e33e0c4a 759 fNonLinearityParams[0] = 9.81039e-01;
760 fNonLinearityParams[1] = 1.13508e-01;
761 fNonLinearityParams[2] = 1.00173e+00;
762 fNonLinearityParams[3] = 9.67998e-02;
763 fNonLinearityParams[4] = 2.19381e+02;
764 fNonLinearityParams[5] = 6.31604e+01;
765 fNonLinearityParams[6] = 1;
766 }
767
8517fbd8 768 if (fNonLinearityFunction == kPi0GammaGamma) {
841dbf60 769 fNonLinearityParams[0] = 1.04;
770 fNonLinearityParams[1] = -0.1445;
771 fNonLinearityParams[2] = 1.046;
fbddd006 772 }
841dbf60 773
8517fbd8 774 if (fNonLinearityFunction == kPi0GammaConversion) {
841dbf60 775 fNonLinearityParams[0] = 0.139393;
776 fNonLinearityParams[1] = 0.0566186;
777 fNonLinearityParams[2] = 0.982133;
fbddd006 778 }
841dbf60 779
8517fbd8 780 if (fNonLinearityFunction == kBeamTest) {
781 if (fNonLinearThreshold == 30) {
841dbf60 782 fNonLinearityParams[0] = 1.007;
783 fNonLinearityParams[1] = 0.894;
784 fNonLinearityParams[2] = 0.246;
785 }
8517fbd8 786 if (fNonLinearThreshold == 45) {
841dbf60 787 fNonLinearityParams[0] = 1.003;
788 fNonLinearityParams[1] = 0.719;
789 fNonLinearityParams[2] = 0.334;
790 }
8517fbd8 791 if (fNonLinearThreshold == 75) {
841dbf60 792 fNonLinearityParams[0] = 1.002;
793 fNonLinearityParams[1] = 0.797;
794 fNonLinearityParams[2] = 0.358;
795 }
796 }
797
8517fbd8 798 if (fNonLinearityFunction == kBeamTestCorrected) {
841dbf60 799 fNonLinearityParams[0] = 0.99078;
800 fNonLinearityParams[1] = 0.161499;
801 fNonLinearityParams[2] = 0.655166;
802 fNonLinearityParams[3] = 0.134101;
803 fNonLinearityParams[4] = 163.282;
804 fNonLinearityParams[5] = 23.6904;
805 fNonLinearityParams[6] = 0.978;
806 }
a31af82c 807
8517fbd8 808 if (fNonLinearityFunction == kBeamTestCorrectedv2) {
a31af82c 809 fNonLinearityParams[0] = 0.983504;
810 fNonLinearityParams[1] = 0.210106;
811 fNonLinearityParams[2] = 0.897274;
812 fNonLinearityParams[3] = 0.0829064;
813 fNonLinearityParams[4] = 152.299;
814 fNonLinearityParams[5] = 31.5028;
815 fNonLinearityParams[6] = 0.968;
816 }
d24a0a76 817
818 if (fNonLinearityFunction == kSDMv5) {
819 fNonLinearityParams[0] = 1.0;
820 fNonLinearityParams[1] = 6.64778e-02;
821 fNonLinearityParams[2] = 1.570;
822 fNonLinearityParams[3] = 9.67998e-02;
823 fNonLinearityParams[4] = 2.19381e+02;
824 fNonLinearityParams[5] = 6.31604e+01;
825 fNonLinearityParams[6] = 1.01286;
826 }
827
828 if (fNonLinearityFunction == kPi0MCv5) {
829 fNonLinearityParams[0] = 1.0;
830 fNonLinearityParams[1] = 6.64778e-02;
831 fNonLinearityParams[2] = 1.570;
832 fNonLinearityParams[3] = 9.67998e-02;
833 fNonLinearityParams[4] = 2.19381e+02;
834 fNonLinearityParams[5] = 6.31604e+01;
835 fNonLinearityParams[6] = 1.01286;
836 }
837
7e0ecb89 838}
839
a520bcd0 840//_________________________________________________________
fb8eab96 841Float_t AliEMCALRecoUtils::GetDepth(Float_t energy,
842 Int_t iParticle,
843 Int_t iSM) const
094786cc 844{
845 //Calculate shower depth for a given cluster energy and particle type
846
847 // parameters
cb231979 848 Float_t x0 = 1.31;
094786cc 849 Float_t ecr = 8;
850 Float_t depth = 0;
a840d589 851 Float_t arg = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV
094786cc 852
853 switch ( iParticle )
854 {
855 case kPhoton:
a840d589 856 if (arg < 1)
857 depth = 0;
858 else
859 depth = x0 * (TMath::Log(arg) + 0.5);
094786cc 860 break;
861
862 case kElectron:
a840d589 863 if (arg < 1)
864 depth = 0;
865 else
866 depth = x0 * (TMath::Log(arg) - 0.5);
094786cc 867 break;
868
869 case kHadron:
870 // hadron
871 // boxes anc. here
8517fbd8 872 if (gGeoManager) {
094786cc 873 gGeoManager->cd("ALIC_1/XEN1_1");
874 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
875 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
8517fbd8 876 if (geoSM) {
fd6df01c 877 TGeoVolume *geoSMVol = geoSM->GetVolume();
878 TGeoShape *geoSMShape = geoSMVol->GetShape();
879 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
8517fbd8 880 if (geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
fd6df01c 881 else AliFatal("Null GEANT box");
7f5392da 882 }
883 else AliFatal("NULL GEANT node matrix");
094786cc 884 }
7f5392da 885 else
886 {//electron
a840d589 887 if (arg < 1)
888 depth = 0;
889 else
890 depth = x0 * (TMath::Log(arg) - 0.5);
094786cc 891 }
892
893 break;
894
895 default://photon
a840d589 896 if (arg < 1)
897 depth = 0;
898 else
899 depth = x0 * (TMath::Log(arg) + 0.5);
094786cc 900 }
901
902 return depth;
094786cc 903}
904
88b96ad8 905//____________________________________________________________________
906void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
907 AliVCaloCells* cells,
908 const AliVCluster* clu,
a520bcd0 909 Int_t & absId,
910 Int_t & iSupMod,
911 Int_t & ieta,
912 Int_t & iphi,
913 Bool_t & shared)
d9b3567c 914{
915 //For a given CaloCluster gets the absId of the cell
916 //with maximum energy deposit.
917
918 Double_t eMax = -1.;
919 Double_t eCell = -1.;
094786cc 920 Float_t fraction = 1.;
921 Float_t recalFactor = 1.;
d9b3567c 922 Int_t cellAbsId = -1 ;
094786cc 923
d9b3567c 924 Int_t iTower = -1;
925 Int_t iIphi = -1;
926 Int_t iIeta = -1;
cb231979 927 Int_t iSupMod0= -1;
2aeb4226 928
8517fbd8 929 if (!clu) {
2aeb4226 930 AliInfo("Cluster pointer null!");
931 absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
932 return;
933 }
934
8517fbd8 935 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 936 cellAbsId = clu->GetCellAbsId(iDig);
937 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 938 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
8517fbd8 939 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 940 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
941 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
8517fbd8 942 if (iDig==0) {
7f5392da 943 iSupMod0=iSupMod;
8517fbd8 944 } else if (iSupMod0!=iSupMod) {
cb231979 945 shared = kTRUE;
946 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
947 }
8517fbd8 948 if (!fCellsRecalibrated && IsRecalibrationOn()) {
094786cc 949 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
950 }
951 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 952 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
8517fbd8 953 if (eCell > eMax) {
d9b3567c 954 eMax = eCell;
955 absId = cellAbsId;
956 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
957 }
958 }// cell loop
959
960 //Get from the absid the supermodule, tower and eta/phi numbers
961 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
962 //Gives SuperModule and Tower numbers
963 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 964 iIphi, iIeta,iphi,ieta);
965 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
966 //printf("Max end---\n");
d9b3567c 967}
968
88b96ad8 969//______________________________________
970void AliEMCALRecoUtils::InitParameters()
971{
972 // Initialize data members with default values
973
974 fParticleType = kPhoton;
975 fPosAlgo = kUnchanged;
976 fW0 = 4.5;
977
978 fNonLinearityFunction = kNoCorrection;
979 fNonLinearThreshold = 30;
980
981 fExoticCellFraction = 0.97;
982 fExoticCellDiffTime = 1e6;
31ba2b24 983 fExoticCellMinAmplitude = 4.0;
88b96ad8 984
a6a1e3ab 985 fAODFilterMask = 128;
986 fAODHybridTracks = kFALSE;
987 fAODTPCOnlyTracks = kTRUE;
88b96ad8 988
989 fCutEtaPhiSum = kTRUE;
990 fCutEtaPhiSeparate = kFALSE;
991
992 fCutR = 0.05;
993 fCutEta = 0.025;
994 fCutPhi = 0.05;
995
996 fClusterWindow = 100;
997 fMass = 0.139;
998
999 fStepSurface = 20.;
1000 fStepCluster = 5.;
1001 fTrackCutsType = kLooseCut;
1002
1003 fCutMinTrackPt = 0;
1004 fCutMinNClusterTPC = -1;
1005 fCutMinNClusterITS = -1;
1006
1007 fCutMaxChi2PerClusterTPC = 1e10;
1008 fCutMaxChi2PerClusterITS = 1e10;
1009
1010 fCutRequireTPCRefit = kFALSE;
1011 fCutRequireITSRefit = kFALSE;
1012 fCutAcceptKinkDaughters = kFALSE;
1013
1014 fCutMaxDCAToVertexXY = 1e10;
1015 fCutMaxDCAToVertexZ = 1e10;
1016 fCutDCAToVertex2D = kFALSE;
1017
42ceff04 1018 fCutRequireITSStandAlone = kFALSE; //MARCEL
1019 fCutRequireITSpureSA = kFALSE; //Marcel
88b96ad8 1020
1021 //Misalignment matrices
8517fbd8 1022 for (Int_t i = 0; i < 15 ; i++)
7f5392da 1023 {
88b96ad8 1024 fMisalTransShift[i] = 0.;
1025 fMisalRotShift[i] = 0.;
1026 }
1027
1028 //Non linearity
8517fbd8 1029 for (Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
88b96ad8 1030
a31af82c 1031 //For kBeamTestCorrectedv2 case, but default is no correction
1032 fNonLinearityParams[0] = 0.983504;
1033 fNonLinearityParams[1] = 0.210106;
1034 fNonLinearityParams[2] = 0.897274;
1035 fNonLinearityParams[3] = 0.0829064;
1036 fNonLinearityParams[4] = 152.299;
1037 fNonLinearityParams[5] = 31.5028;
1038 fNonLinearityParams[6] = 0.968;
88b96ad8 1039
1040 //Cluster energy smearing
1041 fSmearClusterEnergy = kFALSE;
1042 fSmearClusterParam[0] = 0.07; // * sqrt E term
1043 fSmearClusterParam[1] = 0.00; // * E term
1044 fSmearClusterParam[2] = 0.00; // constant
88b96ad8 1045}
1046
1047//_____________________________________________________
1048void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
1049{
841dbf60 1050 //Init EMCAL recalibration factors
1051 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
fbddd006 1052 //In order to avoid rewriting the same histograms
841dbf60 1053 Bool_t oldStatus = TH1::AddDirectoryStatus();
1054 TH1::AddDirectory(kFALSE);
1055
c3d9f926 1056 fEMCALRecalibrationFactors = new TObjArray(12);
1057 for (int i = 0; i < 12; i++)
841dbf60 1058 fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1059 Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1060 //Init the histograms with 1
7f5392da 1061 for (Int_t sm = 0; sm < 12; sm++)
1062 {
1063 for (Int_t i = 0; i < 48; i++)
1064 {
1065 for (Int_t j = 0; j < 24; j++)
1066 {
841dbf60 1067 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1068 }
1069 }
1070 }
7f5392da 1071
841dbf60 1072 fEMCALRecalibrationFactors->SetOwner(kTRUE);
1073 fEMCALRecalibrationFactors->Compress();
fbddd006 1074
841dbf60 1075 //In order to avoid rewriting the same histograms
fbddd006 1076 TH1::AddDirectory(oldStatus);
094786cc 1077}
1078
a520bcd0 1079//_________________________________________________________
841dbf60 1080void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1081{
1082 //Init EMCAL recalibration factors
1083 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1084 //In order to avoid rewriting the same histograms
1085 Bool_t oldStatus = TH1::AddDirectoryStatus();
1086 TH1::AddDirectory(kFALSE);
1087
1088 fEMCALTimeRecalibrationFactors = new TObjArray(4);
1089 for (int i = 0; i < 4; i++)
3bfc4732 1090 fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1091 Form("hAllTimeAvBC%d",i),
c3d9f926 1092 48*24*12,0.,48*24*12) );
841dbf60 1093 //Init the histograms with 1
7f5392da 1094 for (Int_t bc = 0; bc < 4; bc++)
1095 {
c3d9f926 1096 for (Int_t i = 0; i < 48*24*12; i++)
841dbf60 1097 SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
3bfc4732 1098 }
841dbf60 1099
1100 fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1101 fEMCALTimeRecalibrationFactors->Compress();
1102
1103 //In order to avoid rewriting the same histograms
fbddd006 1104 TH1::AddDirectory(oldStatus);
3bfc4732 1105}
094786cc 1106
a520bcd0 1107//____________________________________________________
841dbf60 1108void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1109{
1110 //Init EMCAL bad channels map
1111 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1112 //In order to avoid rewriting the same histograms
1113 Bool_t oldStatus = TH1::AddDirectoryStatus();
608c80a3 1114 TH1::AddDirectory(kFALSE);
1115
c3d9f926 1116 fEMCALBadChannelMap = new TObjArray(12);
841dbf60 1117 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
7f5392da 1118 for (int i = 0; i < 12; i++)
1119 {
841dbf60 1120 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1121 }
1122
1123 fEMCALBadChannelMap->SetOwner(kTRUE);
1124 fEMCALBadChannelMap->Compress();
1125
1126 //In order to avoid rewriting the same histograms
fbddd006 1127 TH1::AddDirectory(oldStatus);
fd6df01c 1128}
1129
88b96ad8 1130//____________________________________________________________________________
1131void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1132 AliVCluster * cluster,
1133 AliVCaloCells * cells,
fb8eab96 1134 Int_t bc)
88b96ad8 1135{
841dbf60 1136 // Recalibrate the cluster energy and Time, considering the recalibration map
3bfc4732 1137 // and the energy of the cells and time that compose the cluster.
1138 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
fbddd006 1139
8517fbd8 1140 if (!cluster) {
2aeb4226 1141 AliInfo("Cluster pointer null!");
1142 return;
1143 }
1144
841dbf60 1145 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1146 UShort_t * index = cluster->GetCellsAbsId() ;
1147 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1148 Int_t ncells = cluster->GetNCells();
fbddd006 1149
841dbf60 1150 //Initialize some used variables
1151 Float_t energy = 0;
1152 Int_t absId =-1;
3bfc4732 1153 Int_t icol =-1, irow =-1, imod=1;
841dbf60 1154 Float_t factor = 1, frac = 0;
3bfc4732 1155 Int_t absIdMax = -1;
1156 Float_t emax = 0;
1157
841dbf60 1158 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
8517fbd8 1159 for (Int_t icell = 0; icell < ncells; icell++)
7f5392da 1160 {
841dbf60 1161 absId = index[icell];
1162 frac = fraction[icell];
8517fbd8 1163 if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
3bfc4732 1164
8517fbd8 1165 if (!fCellsRecalibrated && IsRecalibrationOn()) {
3bfc4732 1166 // Energy
1167 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1168 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
8517fbd8 1169 if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1170 continue;
fbddd006 1171 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
3bfc4732 1172 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1173
1174 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1175 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1176
1177 }
1178
841dbf60 1179 energy += cells->GetCellAmplitude(absId)*factor*frac;
3bfc4732 1180
8517fbd8 1181 if (emax < cells->GetCellAmplitude(absId)*factor*frac) {
3bfc4732 1182 emax = cells->GetCellAmplitude(absId)*factor*frac;
1183 absIdMax = absId;
1184 }
841dbf60 1185 }
fbddd006 1186
3a5708cd 1187 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
3bfc4732 1188
3a5708cd 1189 cluster->SetE(energy);
3bfc4732 1190
bc8330e9 1191 // Recalculate time of cluster
3a5708cd 1192 Double_t timeorg = cluster->GetTOF();
bc8330e9 1193
1194 Double_t time = cells->GetCellTime(absIdMax);
8517fbd8 1195 if (!fCellsRecalibrated && IsTimeRecalibrationOn())
3a5708cd 1196 RecalibrateCellTime(absIdMax,bc,time);
3bfc4732 1197
bc8330e9 1198 cluster->SetTOF(time);
3bfc4732 1199
bc8330e9 1200 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
3bfc4732 1201}
1202
a520bcd0 1203//_____________________________________________________________
1204void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells,
fb8eab96 1205 Int_t bc)
841dbf60 1206{
1207 // Recalibrate the cells time and energy, considering the recalibration map and the energy
3bfc4732 1208 // of the cells that compose the cluster.
1209 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1210
8517fbd8 1211 if (!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn())
1212 return;
3bfc4732 1213
8517fbd8 1214 if (!cells) {
3bfc4732 1215 AliInfo("Cells pointer null!");
1216 return;
1217 }
1218
fbddd006 1219 Short_t absId =-1;
a7e5a381 1220 Bool_t accept = kFALSE;
1221 Float_t ecell = 0;
1222 Double_t tcell = 0;
fbddd006 1223 Double_t ecellin = 0;
1224 Double_t tcellin = 0;
60d77596 1225 Int_t mclabel = -1;
fbddd006 1226 Double_t efrac = 0;
3bfc4732 1227
841dbf60 1228 Int_t nEMcell = cells->GetNumberOfCells() ;
7f5392da 1229 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1230 {
fbddd006 1231 cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
3bfc4732 1232
a7e5a381 1233 accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
8517fbd8 1234 if (!accept) {
a7e5a381 1235 ecell = 0;
fbddd006 1236 tcell = -1;
a7e5a381 1237 }
3bfc4732 1238
1239 //Set new values
fbddd006 1240 cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
3bfc4732 1241 }
841dbf60 1242
1243 fCellsRecalibrated = kTRUE;
094786cc 1244}
1245
88b96ad8 1246//_______________________________________________________________________________________________________
fb8eab96 1247void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime) const
7d692da6 1248{
841dbf60 1249 // Recalibrate time of cell with absID considering the recalibration map
3bfc4732 1250 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
7d692da6 1251
8517fbd8 1252 if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
7d692da6 1253 celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; ;
3bfc4732 1254 }
3bfc4732 1255}
1256
a520bcd0 1257//______________________________________________________________________________
1258void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1259 AliVCaloCells* cells,
1260 AliVCluster* clu)
d9b3567c 1261{
1262 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1263
8517fbd8 1264 if (!clu) {
2aeb4226 1265 AliInfo("Cluster pointer null!");
1266 return;
1267 }
1268
8517fbd8 1269 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1270 else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1271 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 1272}
1273
a520bcd0 1274//_____________________________________________________________________________________________
1275void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom,
1276 AliVCaloCells* cells,
1277 AliVCluster* clu)
094786cc 1278{
1279 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1280 // The algorithm is a copy of what is done in AliEMCALRecPoint
1281
1282 Double_t eCell = 0.;
1283 Float_t fraction = 1.;
1284 Float_t recalFactor = 1.;
1285
1286 Int_t absId = -1;
1287 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1288 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1289 Float_t weight = 0., totalWeight=0.;
1290 Float_t newPos[3] = {0,0,0};
1291 Double_t pLocal[3], pGlobal[3];
cb231979 1292 Bool_t shared = kFALSE;
1293
094786cc 1294 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
a840d589 1295 if (clEnergy <= 0)
1296 return;
cb231979 1297 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1298 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1299
83bfd77a 1300 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1301
7f5392da 1302 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1303 {
841dbf60 1304 absId = clu->GetCellAbsId(iDig);
1305 fraction = clu->GetCellAmplitudeFraction(iDig);
8517fbd8 1306 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1307
8517fbd8 1308 if (!fCellsRecalibrated) {
3bfc4732 1309 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
fbddd006 1310 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
8517fbd8 1311 if (IsRecalibrationOn()) {
3bfc4732 1312 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1313 }
094786cc 1314 }
3bfc4732 1315
094786cc 1316 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1317
1318 weight = GetCellWeight(eCell,clEnergy);
1319 totalWeight += weight;
3bfc4732 1320
094786cc 1321 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 1322 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 1323 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 1324 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1325
8517fbd8 1326 for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
094786cc 1327 }// cell loop
1328
8517fbd8 1329 if (totalWeight>0) {
1330 for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
094786cc 1331 }
1332
094786cc 1333 //Float_t pos[]={0,0,0};
1334 //clu->GetPosition(pos);
1335 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 1336 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 1337
8517fbd8 1338 if (iSupModMax > 1) { //sector 1
841dbf60 1339 newPos[0] +=fMisalTransShift[3];//-=3.093;
1340 newPos[1] +=fMisalTransShift[4];//+=6.82;
1341 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 1342 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
8517fbd8 1343 } else { //sector 0
841dbf60 1344 newPos[0] +=fMisalTransShift[0];//+=1.134;
1345 newPos[1] +=fMisalTransShift[1];//+=8.2;
1346 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 1347 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
841dbf60 1348 }
83bfd77a 1349 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1350
094786cc 1351 clu->SetPosition(newPos);
094786cc 1352}
1353
a520bcd0 1354//____________________________________________________________________________________________
1355void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom,
1356 AliVCaloCells* cells,
1357 AliVCluster* clu)
094786cc 1358{
1359 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1360 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1361
1362 Double_t eCell = 1.;
1363 Float_t fraction = 1.;
1364 Float_t recalFactor = 1.;
1365
1366 Int_t absId = -1;
d9b3567c 1367 Int_t iTower = -1;
094786cc 1368 Int_t iIphi = -1, iIeta = -1;
841dbf60 1369 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 1370 Int_t iphi = -1, ieta =-1;
cb231979 1371 Bool_t shared = kFALSE;
1372
d9b3567c 1373 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
a1c0f640 1374
a840d589 1375 if (clEnergy <= 0)
1376 return;
cb231979 1377 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1378 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1379
d9b3567c 1380 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 1381 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1382 Int_t startingSM = -1;
d9b3567c 1383
7f5392da 1384 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1385 {
094786cc 1386 absId = clu->GetCellAbsId(iDig);
1387 fraction = clu->GetCellAmplitudeFraction(iDig);
8517fbd8 1388 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1389
8517fbd8 1390 if (iDig==0) startingSM = iSupMod;
1391 else if (iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 1392
1393 eCell = cells->GetCellAmplitude(absId);
d9b3567c 1394
3bfc4732 1395 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
fbddd006 1396 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
3bfc4732 1397
7f5392da 1398 if (!fCellsRecalibrated)
1399 {
8517fbd8 1400 if (IsRecalibrationOn()) {
3bfc4732 1401 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3bfc4732 1402 }
094786cc 1403 }
3bfc4732 1404
094786cc 1405 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 1406
094786cc 1407 weight = GetCellWeight(eCell,clEnergy);
8517fbd8 1408 if (weight < 0) weight = 0;
d9b3567c 1409 totalWeight += weight;
1410 weightedCol += ieta*weight;
1411 weightedRow += iphi*weight;
1412
1413 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
841dbf60 1414 }// cell loop
094786cc 1415
d9b3567c 1416 Float_t xyzNew[]={0.,0.,0.};
8517fbd8 1417 if (areInSameSM == kTRUE) {
d9b3567c 1418 //printf("In Same SM\n");
1419 weightedCol = weightedCol/totalWeight;
1420 weightedRow = weightedRow/totalWeight;
094786cc 1421 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
7f5392da 1422 }
1423 else
1424 {
d9b3567c 1425 //printf("In Different SM\n");
094786cc 1426 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 1427 }
d9b3567c 1428
094786cc 1429 clu->SetPosition(xyzNew);
d9b3567c 1430}
1431
a520bcd0 1432//___________________________________________________________________________________________
1433void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
1434 AliVCaloCells* cells,
1435 AliVCluster * cluster)
841dbf60 1436{
cb231979 1437 //re-evaluate distance to bad channel with updated bad map
1438
8517fbd8 1439 if (!fRecalDistToBadChannels) return;
cb231979 1440
8517fbd8 1441 if (!cluster)
7f5392da 1442 {
2aeb4226 1443 AliInfo("Cluster pointer null!");
1444 return;
1445 }
1446
841dbf60 1447 //Get channels map of the supermodule where the cluster is.
fbddd006 1448 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
cb231979 1449 Bool_t shared = kFALSE;
1450 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1451 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1452
fbddd006 1453 Int_t dRrow, dRcol;
841dbf60 1454 Float_t minDist = 10000.;
1455 Float_t dist = 0.;
cb231979 1456
1457 //Loop on tower status map
8517fbd8 1458 for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
7f5392da 1459 {
8517fbd8 1460 for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
7f5392da 1461 {
841dbf60 1462 //Check if tower is bad.
8517fbd8 1463 if (hMap->GetBinContent(icol,irow)==0) continue;
cb231979 1464 //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 1465 // iSupMod,icol, irow, icolM,irowM);
cb231979 1466
1467 dRrow=TMath::Abs(irowM-irow);
1468 dRcol=TMath::Abs(icolM-icol);
1469 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
8517fbd8 1470 if (dist < minDist)
7f5392da 1471 {
cb231979 1472 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1473 minDist = dist;
1474 }
841dbf60 1475 }
1476 }
cb231979 1477
841dbf60 1478 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
7f5392da 1479 if (shared)
1480 {
841dbf60 1481 TH2D* hMap2 = 0;
1482 Int_t iSupMod2 = -1;
cb231979 1483
841dbf60 1484 //The only possible combinations are (0,1), (2,3) ... (8,9)
8517fbd8 1485 if (iSupMod%2) iSupMod2 = iSupMod-1;
1486 else iSupMod2 = iSupMod+1;
841dbf60 1487 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
cb231979 1488
841dbf60 1489 //Loop on tower status map of second super module
8517fbd8 1490 for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
7f5392da 1491 {
8517fbd8 1492 for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
7f5392da 1493 {
841dbf60 1494 //Check if tower is bad.
8517fbd8 1495 if (hMap2->GetBinContent(icol,irow)==0)
1496 continue;
841dbf60 1497 //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",
1498 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
cb231979 1499 dRrow=TMath::Abs(irow-irowM);
1500
8517fbd8 1501 if (iSupMod%2) {
841dbf60 1502 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
8517fbd8 1503 } else {
cb231979 1504 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
841dbf60 1505 }
cb231979 1506
841dbf60 1507 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
8517fbd8 1508 if (dist < minDist) minDist = dist;
841dbf60 1509 }
1510 }
1511 }// shared cluster in 2 SuperModules
78467229 1512
6fe0e6d0 1513 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1514 cluster->SetDistanceToBadChannel(minDist);
cb231979 1515}
1516
a520bcd0 1517//__________________________________________________________________
841dbf60 1518void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1519{
83bfd77a 1520 //re-evaluate identification parameters with bayesian
2aeb4226 1521
8517fbd8 1522 if (!cluster) {
2aeb4226 1523 AliInfo("Cluster pointer null!");
1524 return;
1525 }
1526
8517fbd8 1527 if (cluster->GetM02() != 0)
83bfd77a 1528 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1529
00a38d07 1530 Float_t pidlist[AliPID::kSPECIESCN+1];
8517fbd8 1531 for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
841dbf60 1532
83bfd77a 1533 cluster->SetPID(pidlist);
83bfd77a 1534}
1535
f0e9e976 1536//___________________________________________________________________________________________________________________
a520bcd0 1537void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1538 AliVCaloCells* cells,
f0e9e976 1539 AliVCluster * cluster,
1540 Float_t & l0, Float_t & l1,
1541 Float_t & disp, Float_t & dEta, Float_t & dPhi,
1542 Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
83bfd77a 1543{
1544 // Calculates new center of gravity in the local EMCAL-module coordinates
1545 // and tranfers into global ALICE coordinates
1546 // Calculates Dispersion and main axis
1547
8517fbd8 1548 if (!cluster) {
2aeb4226 1549 AliInfo("Cluster pointer null!");
1550 return;
1551 }
1552
83bfd77a 1553 Double_t eCell = 0.;
1554 Float_t fraction = 1.;
1555 Float_t recalFactor = 1.;
1556
f0e9e976 1557 Int_t iSupMod = -1;
1558 Int_t iTower = -1;
1559 Int_t iIphi = -1;
1560 Int_t iIeta = -1;
1561 Int_t iphi = -1;
1562 Int_t ieta = -1;
1563 Double_t etai = -1.;
1564 Double_t phii = -1.;
1565
1566 Int_t nstat = 0 ;
1567 Float_t wtot = 0.;
1568 Double_t w = 0.;
1569 Double_t etaMean = 0.;
1570 Double_t phiMean = 0.;
a1c0f640 1571
1572 //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added
1573 // and to check if the cluster is between 2 SM in eta
1574 Int_t iSM0 = -1;
1575 Bool_t shared = kFALSE;
1576 Float_t energy = 0;
1577
8517fbd8 1578 for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
a1c0f640 1579 {
1580 //Get from the absid the supermodule, tower and eta/phi numbers
1581 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1582 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1583
1584 //Check if there are cells of different SM
8517fbd8 1585 if (iDigit == 0 ) iSM0 = iSupMod;
1586 else if (iSupMod!= iSM0) shared = kTRUE;
a1c0f640 1587
1588 //Get the cell energy, if recalibration is on, apply factors
1589 fraction = cluster->GetCellAmplitudeFraction(iDigit);
8517fbd8 1590 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
83bfd77a 1591
8517fbd8 1592 if (IsRecalibrationOn()) {
a1c0f640 1593 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1594 }
1595
1596 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1597
1598 energy += eCell;
1599
1600 }//cell loop
1601
83bfd77a 1602 //Loop on cells
8517fbd8 1603 for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
7f5392da 1604 {
83bfd77a 1605 //Get from the absid the supermodule, tower and eta/phi numbers
1606 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1607 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1608
1609 //Get the cell energy, if recalibration is on, apply factors
1610 fraction = cluster->GetCellAmplitudeFraction(iDigit);
8517fbd8 1611 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1612
8517fbd8 1613 if (!fCellsRecalibrated) {
1614 if (IsRecalibrationOn()) {
3bfc4732 1615 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1616 }
83bfd77a 1617 }
3bfc4732 1618
83bfd77a 1619 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1620
a1c0f640 1621 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1622 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
8517fbd8 1623 if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
a1c0f640 1624
8517fbd8 1625 if (cluster->E() > 0 && eCell > 0) {
83bfd77a 1626 w = GetCellWeight(eCell,cluster->E());
1627
1628 etai=(Double_t)ieta;
fbddd006 1629 phii=(Double_t)iphi;
f0e9e976 1630
8517fbd8 1631 if (w > 0.0) {
83bfd77a 1632 wtot += w ;
fbddd006 1633 nstat++;
83bfd77a 1634 //Shower shape
f0e9e976 1635 sEta += w * etai * etai ;
1636 etaMean += w * etai ;
1637 sPhi += w * phii * phii ;
1638 phiMean += w * phii ;
1639 sEtaPhi += w * etai * phii ;
83bfd77a 1640 }
8517fbd8 1641 } else
83bfd77a 1642 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1643 }//cell loop
1644
fbddd006 1645 //Normalize to the weight
8517fbd8 1646 if (wtot > 0) {
f0e9e976 1647 etaMean /= wtot ;
1648 phiMean /= wtot ;
8517fbd8 1649 } else
83bfd77a 1650 AliError(Form("Wrong weight %f\n", wtot));
1651
fbddd006 1652 //Calculate dispersion
8517fbd8 1653 for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
7f5392da 1654 {
83bfd77a 1655 //Get from the absid the supermodule, tower and eta/phi numbers
1656 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1657 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1658
1659 //Get the cell energy, if recalibration is on, apply factors
1660 fraction = cluster->GetCellAmplitudeFraction(iDigit);
8517fbd8 1661 if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1662 if (IsRecalibrationOn()) {
83bfd77a 1663 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1664 }
1665 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1666
a1c0f640 1667 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1668 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
8517fbd8 1669 if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
a1c0f640 1670
8517fbd8 1671 if (cluster->E() > 0 && eCell > 0) {
83bfd77a 1672 w = GetCellWeight(eCell,cluster->E());
1673
1674 etai=(Double_t)ieta;
fbddd006 1675 phii=(Double_t)iphi;
8517fbd8 1676 if (w > 0.0) {
f0e9e976 1677 disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
1678 dEta += w * (etai-etaMean)*(etai-etaMean) ;
1679 dPhi += w * (phii-phiMean)*(phii-phiMean) ;
1680 }
8517fbd8 1681 } else
83bfd77a 1682 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1683 }// cell loop
1684
1685 //Normalize to the weigth and set shower shape parameters
8517fbd8 1686 if (wtot > 0 && nstat > 1) {
f0e9e976 1687 disp /= wtot ;
1688 dEta /= wtot ;
1689 dPhi /= wtot ;
1690 sEta /= wtot ;
1691 sPhi /= wtot ;
1692 sEtaPhi /= wtot ;
1693
1694 sEta -= etaMean * etaMean ;
1695 sPhi -= phiMean * phiMean ;
1696 sEtaPhi -= etaMean * phiMean ;
1697
1698 l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1699 l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
8517fbd8 1700 } else {
f0e9e976 1701 l0 = 0. ;
1702 l1 = 0. ;
1703 dEta = 0. ; dPhi = 0. ; disp = 0. ;
1704 sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
fbddd006 1705 }
83bfd77a 1706}
1707
f0e9e976 1708//____________________________________________________________________________________________
1709void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1710 AliVCaloCells* cells,
1711 AliVCluster * cluster)
1712{
1713 // Calculates new center of gravity in the local EMCAL-module coordinates
1714 // and tranfers into global ALICE coordinates
1715 // Calculates Dispersion and main axis and puts them into the cluster
1716
1717 Float_t l0 = 0., l1 = 0.;
1718 Float_t disp = 0., dEta = 0., dPhi = 0.;
1719 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1720
1721 AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
1722 dEta, dPhi, sEta, sPhi, sEtaPhi);
1723
1724 cluster->SetM02(l0);
1725 cluster->SetM20(l1);
8517fbd8 1726 if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
f0e9e976 1727
1728}
1729
b540d03f 1730//____________________________________________________________________________
a520bcd0 1731void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
1732 TObjArray * clusterArr,
1733 const AliEMCALGeometry *geom)
bd8c7aef 1734{
1735 //This function should be called before the cluster loop
1736 //Before call this function, please recalculate the cluster positions
1737 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1738 //Store matched cluster indexes and residuals
61160f1f 1739
88b96ad8 1740 fMatchedTrackIndex ->Reset();
bd8c7aef 1741 fMatchedClusterIndex->Reset();
fa4287a2 1742 fResidualPhi->Reset();
1743 fResidualEta->Reset();
bd8c7aef 1744
841dbf60 1745 fMatchedTrackIndex ->Set(1000);
1746 fMatchedClusterIndex->Set(1000);
1747 fResidualPhi->Set(1000);
1748 fResidualEta->Set(1000);
bd8c7aef 1749
1c7a2bf4 1750 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1751 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
57131575 1752
88b96ad8 1753 // init the magnetic field if not already on
8517fbd8 1754 if (!TGeoGlobalMagField::Instance()->GetField()) {
092f179b 1755 if (!event->InitMagneticField())
88b96ad8 1756 {
1757 AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1758 }
88b96ad8 1759 } // Init mag field
1760
42ceff04 1761 if (esdevent) {
1762 UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
1763 UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
1764 Bool_t desc1 = (mask1 >> 3) & 0x1;
1765 Bool_t desc2 = (mask2 >> 3) & 0x1;
1766 if (desc1==0 || desc2==0) {
30e29b2a 1767// AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
1768// mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
1769// mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
42ceff04 1770 fITSTrackSA=kTRUE;
1771 }
1772 }
1773
8fc351e3 1774 TObjArray *clusterArray = 0x0;
8517fbd8 1775 if (!clusterArr) {
092f179b 1776 clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
8517fbd8 1777 for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
092f179b
CL		 1778 {
1779 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8517fbd8 1780 if (geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
1781 continue;
092f179b 1782 clusterArray->AddAt(cluster,icl);
8fc351e3 1783 }
092f179b 1784 }
61160f1f 1785
bd8c7aef 1786 Int_t matched=0;
bb6f5f0b 1787 Double_t cv[21];
1788 for (Int_t i=0; i<21;i++) cv[i]=0;
8517fbd8 1789 for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
bd8c7aef 1790 {
456126ad 1791 AliExternalTrackParam *trackParam = 0;
61160f1f 1792
bb6f5f0b 1793 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
8fc351e3 1794 AliESDtrack *esdTrack = 0;
1795 AliAODTrack *aodTrack = 0;
8517fbd8 1796 if (esdevent) {
8fc351e3 1797 esdTrack = esdevent->GetTrack(itr);
8517fbd8 1798 if (!esdTrack) continue;
1799 if (!IsAccepted(esdTrack)) continue;
1800 if (esdTrack->Pt()<fCutMinTrackPt) continue;
8fc351e3 1801 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
8517fbd8 1802 if (TMath::Abs(esdTrack->Eta())>0.9 || phi <= 10 || phi >= 250 ) continue;
1803 if (!fITSTrackSA)
42ceff04 1804 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
1805 else
1806 trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
61160f1f 1807 }
bb6f5f0b 1808
1809 //If the input event is AOD, the starting point for extrapolation is at vertex
8fc351e3 1810 //AOD tracks are selected according to its filterbit.
8517fbd8 1811 else if (aodevent) {
f15c1f69 1812 aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
1813 if(!aodTrack) AliFatal("Not a standard AOD");
8517fbd8 1814 if (!aodTrack) continue;
a6a1e3ab 1815
8517fbd8 1816 if (fAODTPCOnlyTracks) { // Match with TPC only tracks, default from May 2013, before filter bit 32
a6a1e3ab 1817 //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
c860bb71 1818 if (!aodTrack->IsTPCConstrained()) continue ;
8517fbd8 1819 } else if (fAODHybridTracks) { // Match with hybrid tracks
a6a1e3ab 1820 //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
8517fbd8 1821 if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
1822 } else { // Match with tracks on a mask
a6a1e3ab 1823 //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
8517fbd8 1824 if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
a6a1e3ab 1825 }
1af378e6 1826
8517fbd8 1827 if (aodTrack->Pt()<fCutMinTrackPt) continue;
1af378e6 1828
8fc351e3 1829 Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
8517fbd8 1830 if (TMath::Abs(aodTrack->Eta())>0.9 || phi <= 10 || phi >= 250 )
1831 continue;
61160f1f 1832 Double_t pos[3],mom[3];
1833 aodTrack->GetXYZ(pos);
1834 aodTrack->GetPxPyPz(mom);
1835 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 1836
61160f1f 1837 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1838 }
bd8c7aef 1839
bb6f5f0b 1840 //Return if the input data is not "AOD" or "ESD"
8517fbd8 1841 else {
61160f1f 1842 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
8517fbd8 1843 if (clusterArray) {
092f179b
CL		 1844 clusterArray->Clear();
1845 delete clusterArray;
1846 }
61160f1f 1847 return;
1848 }
1849
8517fbd8 1850 if (!trackParam) continue;
8fc351e3 1851
1852 //Extrapolate the track to EMCal surface
1853 AliExternalTrackParam emcalParam(*trackParam);
a29b2a8a 1854 Float_t eta, phi, pt;
8517fbd8 1855 if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) {
1856 if (aodevent && trackParam) delete trackParam;
1857 if (fITSTrackSA && trackParam) delete trackParam;
092f179b
CL		 1858 continue;
1859 }
8fc351e3 1860
8517fbd8 1861 if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) {
1862 if (aodevent && trackParam) delete trackParam;
1863 if (fITSTrackSA && trackParam) delete trackParam;
092f179b
CL		 1864 continue;
1865 }
8fc351e3 1866
1867 //Find matched clusters
bd8c7aef 1868 Int_t index = -1;
8fc351e3 1869 Float_t dEta = -999, dPhi = -999;
8517fbd8 1870 if (!clusterArr) {
092f179b 1871 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
8517fbd8 1872 } else {
092f179b
CL		 1873 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1874 }
61160f1f 1875
8517fbd8 1876 if (index>-1) {
8fc351e3 1877 fMatchedTrackIndex ->AddAt(itr,matched);
1878 fMatchedClusterIndex ->AddAt(index,matched);
1879 fResidualEta ->AddAt(dEta,matched);
1880 fResidualPhi ->AddAt(dPhi,matched);
bd8c7aef 1881 matched++;
1882 }
8517fbd8 1883 if (aodevent && trackParam) delete trackParam;
1884 if (fITSTrackSA && trackParam) delete trackParam;
bd8c7aef 1885 }//track loop
8fc351e3 1886
8517fbd8 1887 if (clusterArray) {
092f179b
CL		 1888 clusterArray->Clear();
1889 delete clusterArray;
1890 }
b540d03f 1891
1892 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1893
8fc351e3 1894 fMatchedTrackIndex ->Set(matched);
1895 fMatchedClusterIndex ->Set(matched);
1896 fResidualPhi ->Set(matched);
1897 fResidualEta ->Set(matched);
bd8c7aef 1898}
1899
b540d03f 1900//________________________________________________________________________________
a520bcd0 1901Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
1902 const AliVEvent *event,
1903 const AliEMCALGeometry *geom,
1904 Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1905{
1906 //
1907 // This function returns the index of matched cluster to input track
fa4287a2 1908 // Returns -1 if no match is found
bb6f5f0b 1909 Int_t index = -1;
8fc351e3 1910 Double_t phiV = track->Phi()*TMath::RadToDeg();
8517fbd8 1911 if (TMath::Abs(track->Eta())>0.9 || phiV <= 10 || phiV >= 250 ) return index;
42ceff04 1912 AliExternalTrackParam *trackParam = 0;
8517fbd8 1913 if (!fITSTrackSA)
42ceff04 1914 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
1915 else
1916 trackParam = new AliExternalTrackParam(*track);
1917
8517fbd8 1918 if (!trackParam) return index;
8fc351e3 1919 AliExternalTrackParam emcalParam(*trackParam);
a29b2a8a 1920 Float_t eta, phi, pt;
8fc351e3 1921
8517fbd8 1922 if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) {
1923 if (fITSTrackSA) delete trackParam;
1924 return index;
30e29b2a 1925 }
8517fbd8 1926 if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) {
1927 if (fITSTrackSA) delete trackParam;
1928 return index;
30e29b2a 1929 }
1930
8fc351e3 1931 TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
1932
8517fbd8 1933 for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 1934 {
1935 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8517fbd8 1936 if (geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
8fc351e3 1937 clusterArr->AddAt(cluster,icl);
1938 }
1939
1940 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1941 clusterArr->Clear();
1942 delete clusterArr;
8517fbd8 1943 if (fITSTrackSA) delete trackParam;
30e29b2a 1944
8fc351e3 1945 return index;
1946}
1947
7f5392da 1948//_______________________________________________________________________________________________
1949Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
a520bcd0 1950 AliExternalTrackParam *trkParam,
7f5392da 1951 const TObjArray * clusterArr,
a520bcd0 1952 Float_t &dEta, Float_t &dPhi)
8fc351e3 1953{
7f5392da 1954 // Find matched cluster in array
1955
8fc351e3 1956 dEta=-999, dPhi=-999;
1957 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1958 Int_t index = -1;
ee602376 1959 Float_t tmpEta=-999, tmpPhi=-999;
8fc351e3 1960
1961 Double_t exPos[3] = {0.,0.,0.};
8517fbd8 1962 if (!emcalParam->GetXYZ(exPos)) return index;
8fc351e3 1963
1964 Float_t clsPos[3] = {0.,0.,0.};
8517fbd8 1965 for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
092f179b
CL		 1966 {
1967 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
8517fbd8 1968 if (!cluster || !cluster->IsEMCAL()) continue;
092f179b
CL		 1969 cluster->GetPosition(clsPos);
1970 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 1971 if (dR > fClusterWindow) continue;
092f179b
CL		 1972
1973 AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
8517fbd8 1974 if (!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
1975 if (fCutEtaPhiSum) {
092f179b 1976 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
8517fbd8 1977 if (tmpR<dRMax) {
fbddd006 1978 dRMax=tmpR;
1979 dEtaMax=tmpEta;
1980 dPhiMax=tmpPhi;
1981 index=icl;
1982 }
8517fbd8 1983 } else if (fCutEtaPhiSeparate) {
1984 if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax)) {
fbddd006 1985 dEtaMax = tmpEta;
1986 dPhiMax = tmpPhi;
1987 index=icl;
1988 }
8517fbd8 1989 } else {
092f179b
CL		 1990 printf("Error: please specify your cut criteria\n");
1991 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1992 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1993 return index;
1994 }
1995 }
8fc351e3 1996
1997 dEta=dEtaMax;
1998 dPhi=dPhiMax;
1999
bb6f5f0b 2000 return index;
2001}
2002
40891fa2
CL		 2003//------------------------------------------------------------------------------------
2004Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliVTrack *track,
18a566cc 2005 Double_t emcalR, Double_t mass,
2006 Double_t step, Double_t minpt,
2007 Bool_t useMassForTracking)
40891fa2 2008{
8517fbd8 2009 // Extrapolate track to EMCAL surface
40891fa2
CL		 2010
2011 track->SetTrackPhiEtaPtOnEMCal(-999, -999, -999);
2012
6e0ae02d 2013 if (track->Pt()<minpt)
40891fa2
CL		 2014 return kFALSE;
2015
2016 Double_t phi = track->Phi()*TMath::RadToDeg();
2017 if (TMath::Abs(track->Eta())>0.9 || phi <= 10 || phi >= 250)
2018 return kFALSE;
2019
2020 AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(track);
2021 AliAODTrack *aodt = 0;
2022 if (!esdt) {
2023 aodt = dynamic_cast<AliAODTrack*>(track);
2024 if (!aodt)
2025 return kFALSE;
2026 }
2027
26664f75 2028 // Select the mass hypothesis
2029 if ( mass < 0 )
2030 {
40891fa2 2031 Bool_t onlyTPC = kFALSE;
26664f75 2032 if ( mass == -99 ) onlyTPC=kTRUE;
2033
2034 if (esdt)
2035 {
18a566cc 2036 if ( useMassForTracking ) mass = esdt->GetMassForTracking();
2037 else mass = esdt->GetMass(onlyTPC);
2038 }
2039 else
26664f75 2040 {
2041 if ( useMassForTracking ) mass = aodt->GetMassForTracking();
2042 else mass = aodt->M();
2043 }
40891fa2
CL		 2044 }
2045
2046 AliExternalTrackParam *trackParam = 0;
2047 if (esdt) {
0e1c2338 2048 const AliExternalTrackParam *in = esdt->GetInnerParam();
7af66d73
MP		 2049 if (!in)
2050 return kFALSE;
2051 trackParam = new AliExternalTrackParam(*in);
40891fa2
CL		 2052 } else {
2053 Double_t xyz[3] = {0}, pxpypz[3] = {0}, cv[21] = {0};
2054 aodt->PxPyPz(pxpypz);
2055 aodt->XvYvZv(xyz);
2056 aodt->GetCovarianceXYZPxPyPz(cv);
2057 trackParam = new AliExternalTrackParam(xyz,pxpypz,cv,aodt->Charge());
2058 }
2059 if (!trackParam)
2060 return kFALSE;
2061
2062 Float_t etaout=-999, phiout=-999, ptout=-999;
2063 Bool_t ret = ExtrapolateTrackToEMCalSurface(trackParam,
2064 emcalR,
2065 mass,
2066 step,
2067 etaout,
2068 phiout,
2069 ptout);
3ad1e868 2070 delete trackParam;
40891fa2
CL		 2071 if (!ret)
2072 return kFALSE;
40891fa2
CL		 2073 if (TMath::Abs(etaout)>0.75 || (phiout<70*TMath::DegToRad()) || (phiout>190*TMath::DegToRad()))
2074 return kFALSE;
40891fa2
CL		 2075 track->SetTrackPhiEtaPtOnEMCal(phiout, etaout, ptout);
2076 return kTRUE;
2077}
2078
2079
88b96ad8 2080//------------------------------------------------------------------------------------
2081Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
fb8eab96 2082 Double_t emcalR,
2083 Double_t mass,
2084 Double_t step,
88b96ad8 2085 Float_t &eta,
a29b2a8a 2086 Float_t &phi,
2087 Float_t &pt)
ee602376 2088{
88b96ad8 2089 //Extrapolate track to EMCAL surface
2090
a29b2a8a 2091 eta = -999, phi = -999, pt = -999;
8517fbd8 2092 if (!trkParam) return kFALSE;
2093 if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
ee602376 2094 Double_t trkPos[3] = {0.,0.,0.};
8517fbd8 2095 if (!trkParam->GetXYZ(trkPos)) return kFALSE;
ee602376 2096 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2097 eta = trkPosVec.Eta();
2098 phi = trkPosVec.Phi();
a29b2a8a 2099 pt = trkParam->Pt();
8517fbd8 2100 if (phi<0)
ee602376 2101 phi += 2*TMath::Pi();
2102
2103 return kTRUE;
2104}
2105
88b96ad8 2106//-----------------------------------------------------------------------------------
2107Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2108 const Float_t *clsPos,
2109 Double_t mass,
2110 Double_t step,
2111 Float_t &tmpEta,
2112 Float_t &tmpPhi)
ee602376 2113{
2114 //
2115 //Return the residual by extrapolating a track param to a global position
2116 //
2117 tmpEta = -999;
2118 tmpPhi = -999;
8517fbd8 2119 if (!trkParam) return kFALSE;
ee602376 2120 Double_t trkPos[3] = {0.,0.,0.};
2121 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2122 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2123 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
8517fbd8 2124 if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
2125 if (!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
ee602376 2126
2127 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2128 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2129
2130 // track cluster matching
2131 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
2132 tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2133
2134 return kTRUE;
2135}
2136
88b96ad8 2137//----------------------------------------------------------------------------------
2138Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2139 const AliVCluster *cluster,
fb8eab96 2140 Double_t mass,
2141 Double_t step,
88b96ad8 2142 Float_t &tmpEta,
2143 Float_t &tmpPhi)
ee602376 2144{
2145 //
2146 //Return the residual by extrapolating a track param to a cluster
2147 //
2148 tmpEta = -999;
2149 tmpPhi = -999;
8517fbd8 2150 if (!cluster || !trkParam)
2151 return kFALSE;
ee602376 2152
2153 Float_t clsPos[3] = {0.,0.,0.};
2154 cluster->GetPosition(clsPos);
2155
2156 return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2157}
2158
88b96ad8 2159//---------------------------------------------------------------------------------
2160Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2161 const AliVCluster *cluster,
88b96ad8 2162 Float_t &tmpEta,
2163 Float_t &tmpPhi)
bb6f5f0b 2164{
2165 //
ee602376 2166 //Return the residual by extrapolating a track param to a clusterfStepCluster
bb6f5f0b 2167 //
8fc351e3 2168
ee602376 2169 return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
bb6f5f0b 2170}
2171
a520bcd0 2172//_______________________________________________________________________
fb8eab96 2173void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex,
a520bcd0 2174 Float_t &dEta, Float_t &dPhi)
bd8c7aef 2175{
bb6f5f0b 2176 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 2177 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 2178 //Works with ESDs and AODs
bd8c7aef 2179
8517fbd8 2180 if (FindMatchedPosForCluster(clsIndex) >= 999) {
bd8c7aef 2181 AliDebug(2,"No matched tracks found!\n");
fa4287a2 2182 dEta=999.;
2183 dPhi=999.;
bd8c7aef 2184 return;
2185 }
fa4287a2 2186 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2187 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 2188}
841dbf60 2189
88b96ad8 2190//______________________________________________________________________________________________
fa4287a2 2191void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 2192{
2193 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 2194 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 2195 //Works with ESDs and AODs
2196
8517fbd8 2197 if (FindMatchedPosForTrack(trkIndex) >= 999) {
bb6f5f0b 2198 AliDebug(2,"No matched cluster found!\n");
fa4287a2 2199 dEta=999.;
2200 dPhi=999.;
bb6f5f0b 2201 return;
2202 }
fa4287a2 2203 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2204 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 2205}
2206
2207//__________________________________________________________
2208Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2209{
2210 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2211 //Get the index of matched track to this cluster
2212 //Works with ESDs and AODs
2213
8517fbd8 2214 if (IsClusterMatched(clsIndex))
bb6f5f0b 2215 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2216 else
2217 return -1;
bd8c7aef 2218}
2219
b540d03f 2220//__________________________________________________________
bb6f5f0b 2221Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 2222{
bb6f5f0b 2223 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2224 //Get the index of matched cluster to this track
2225 //Works with ESDs and AODs
b540d03f 2226
8517fbd8 2227 if (IsTrackMatched(trkIndex))
bb6f5f0b 2228 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 2229 else
2230 return -1;
2231}
2232
7f5392da 2233//______________________________________________________________
7cdec71f 2234Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 2235{
2236 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2237 //Returns if the cluster has a match
8517fbd8 2238 if (FindMatchedPosForCluster(clsIndex) < 999)
bb6f5f0b 2239 return kTRUE;
2240 else
2241 return kFALSE;
2242}
b540d03f 2243
7f5392da 2244//____________________________________________________________
7cdec71f 2245Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 2246{
bb6f5f0b 2247 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2248 //Returns if the track has a match
8517fbd8 2249 if (FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 2250 return kTRUE;
bd8c7aef 2251 else
2252 return kFALSE;
2253}
bb6f5f0b 2254
7f5392da 2255//______________________________________________________________________
bb6f5f0b 2256UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 2257{
bb6f5f0b 2258 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 2259 //Returns the position of the match in the fMatchedClusterIndex array
2260 Float_t tmpR = fCutR;
81efb149 2261 UInt_t pos = 999;
b540d03f 2262
8517fbd8 2263 for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
7f5392da 2264 {
8517fbd8 2265 if (fMatchedClusterIndex->At(i)==clsIndex) {
841dbf60 2266 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
8517fbd8 2267 if (r<tmpR) {
841dbf60 2268 pos=i;
2269 tmpR=r;
7f5392da 2270 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2271 fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2272 }
841dbf60 2273 }
bb6f5f0b 2274 }
2275 return pos;
2276}
2277
7f5392da 2278//____________________________________________________________________
bb6f5f0b 2279UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2280{
2281 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2282 //Returns the position of the match in the fMatchedTrackIndex array
2283 Float_t tmpR = fCutR;
2284 UInt_t pos = 999;
2285
8517fbd8 2286 for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
7f5392da 2287 {
8517fbd8 2288 if (fMatchedTrackIndex->At(i)==trkIndex) {
841dbf60 2289 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
8517fbd8 2290 if (r<tmpR) {
841dbf60 2291 pos=i;
2292 tmpR=r;
7f5392da 2293 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2294 fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2295 }
841dbf60 2296 }
b540d03f 2297 }
bd8c7aef 2298 return pos;
2299}
2300
a520bcd0 2301//__________________________________________________________________________
2302Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
2303 const AliEMCALGeometry *geom,
fb8eab96 2304 AliVCaloCells* cells, Int_t bc)
b5078f5d 2305{
2306 // check if the cluster survives some quality cut
2307 //
2308 //
2309 Bool_t isGood=kTRUE;
7f5392da 2310
8517fbd8 2311 if (!cluster || !cluster->IsEMCAL()) return kFALSE;
2312 if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
2313 if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
2314 if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
b5078f5d 2315
2316 return isGood;
2317}
2318
b540d03f 2319//__________________________________________________________
bd8c7aef 2320Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2321{
2322 // Given a esd track, return whether the track survive all the cuts
2323
2324 // The different quality parameter are first
2325 // retrieved from the track. then it is found out what cuts the
2326 // track did not survive and finally the cuts are imposed.
2327
2328 UInt_t status = esdTrack->GetStatus();
2329
2330 Int_t nClustersITS = esdTrack->GetITSclusters(0);
2331 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2332
2333 Float_t chi2PerClusterITS = -1;
2334 Float_t chi2PerClusterTPC = -1;
2335 if (nClustersITS!=0)
2336 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2337 if (nClustersTPC!=0)
2338 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 2339
2340
2341 //DCA cuts
8517fbd8 2342 if (fTrackCutsType==kGlobalCut) {
2343 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2344 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2345 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2346 }
82d09e74 2347
bd8c7aef 2348 Float_t b[2];
2349 Float_t bCov[3];
2350 esdTrack->GetImpactParameters(b,bCov);
8517fbd8 2351 if (bCov[0]<=0 || bCov[2]<=0) {
bd8c7aef 2352 AliDebug(1, "Estimated b resolution lower or equal zero!");
2353 bCov[0]=0; bCov[2]=0;
2354 }
2355
2356 Float_t dcaToVertexXY = b[0];
2357 Float_t dcaToVertexZ = b[1];
2358 Float_t dcaToVertex = -1;
2359
2360 if (fCutDCAToVertex2D)
2361 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2362 else
2363 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2364
2365 // cut the track?
2366
2367 Bool_t cuts[kNCuts];
2368 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2369
2370 // track quality cuts
2371 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2372 cuts[0]=kTRUE;
2373 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2374 cuts[1]=kTRUE;
2375 if (nClustersTPC<fCutMinNClusterTPC)
2376 cuts[2]=kTRUE;
2377 if (nClustersITS<fCutMinNClusterITS)
2378 cuts[3]=kTRUE;
2379 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
2380 cuts[4]=kTRUE;
2381 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
2382 cuts[5]=kTRUE;
2383 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2384 cuts[6]=kTRUE;
2385 if (fCutDCAToVertex2D && dcaToVertex > 1)
2386 cuts[7] = kTRUE;
2387 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2388 cuts[8] = kTRUE;
2389 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2390 cuts[9] = kTRUE;
2391
8517fbd8 2392 if (fTrackCutsType==kGlobalCut) {
2393 //Require at least one SPD point + anything else in ITS
2394 if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
2395 cuts[10] = kTRUE;
2396 }
82d09e74 2397
8517fbd8 2398 // ITS
2399 if (fCutRequireITSStandAlone || fCutRequireITSpureSA) {
2400 if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)) {
42ceff04 2401 // TPC tracks
2402 cuts[11] = kTRUE;
8517fbd8 2403 } else {
42ceff04 2404 // ITS standalone tracks
8517fbd8 2405 if (fCutRequireITSStandAlone && !fCutRequireITSpureSA) {
2406 if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
2407 } else if (fCutRequireITSpureSA) {
2408 if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
42ceff04 2409 }
2410 }
2411 }
2412
bd8c7aef 2413 Bool_t cut=kFALSE;
827f9f23 2414 for (Int_t i=0; i<kNCuts; i++)
7f5392da 2415 if (cuts[i]) { cut = kTRUE ; }
bd8c7aef 2416
2417 // cut the track
2418 if (cut)
2419 return kFALSE;
2420 else
2421 return kTRUE;
2422}
827f9f23 2423
88b96ad8 2424//_____________________________________
bd8c7aef 2425void AliEMCALRecoUtils::InitTrackCuts()
2426{
2427 //Intilize the track cut criteria
5f7714ad 2428 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
bd8c7aef 2429 //Also you can customize the cuts using the setters
82d09e74 2430
5f7714ad 2431 switch (fTrackCutsType)
88b96ad8 2432 {
5f7714ad 2433 case kTPCOnlyCut:
88b96ad8 2434 {
2435 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2436 //TPC
2437 SetMinNClustersTPC(70);
2438 SetMaxChi2PerClusterTPC(4);
2439 SetAcceptKinkDaughters(kFALSE);
2440 SetRequireTPCRefit(kFALSE);
2441
2442 //ITS
2443 SetRequireITSRefit(kFALSE);
2444 SetMaxDCAToVertexZ(3.2);
2445 SetMaxDCAToVertexXY(2.4);
2446 SetDCAToVertex2D(kTRUE);
2447
2448 break;
2449 }
2450
5f7714ad 2451 case kGlobalCut:
88b96ad8 2452 {
2453 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2454 //TPC
2455 SetMinNClustersTPC(70);
2456 SetMaxChi2PerClusterTPC(4);
2457 SetAcceptKinkDaughters(kFALSE);
2458 SetRequireTPCRefit(kTRUE);
2459
2460 //ITS
2461 SetRequireITSRefit(kTRUE);
2462 SetMaxDCAToVertexZ(2);
2463 SetMaxDCAToVertexXY();
2464 SetDCAToVertex2D(kFALSE);
2465
2466 break;
2467 }
2468
0e7de35b 2469 case kLooseCut:
88b96ad8 2470 {
2471 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2472 SetMinNClustersTPC(50);
2473 SetAcceptKinkDaughters(kTRUE);
2474
2475 break;
5f7714ad 2476 }
42ceff04 2477
2478 case kITSStandAlone:
2479 {
2480 AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
2481 SetRequireITSRefit(kTRUE);
2482 SetRequireITSStandAlone(kTRUE);
2483 SetITSTrackSA(kTRUE);
2484 break;
2485 }
2486
88b96ad8 2487 }
bd8c7aef 2488}
83bfd77a 2489
57131575 2490
88b96ad8 2491//________________________________________________________________________
dda65b42 2492void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
57131575 2493{
2494 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
dda65b42 2495
57131575 2496 Int_t nTracks = event->GetNumberOfTracks();
7f5392da 2497 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
2498 {
dda65b42 2499 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
7f5392da 2500 if (!track)
2501 {
57131575 2502 AliWarning(Form("Could not receive track %d", iTrack));
2503 continue;
2504 }
7f5392da 2505
fbddd006 2506 Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
57131575 2507 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
dda65b42 2508 /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2509 AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2510 if (esdtrack) {
8517fbd8 2511 if (matchClusIndex != -1)
dda65b42 2512 esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2513 else
2514 esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2515 } else {
2516 AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
8517fbd8 2517 if (matchClusIndex != -1)
dda65b42 2518 aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2519 else
2520 aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2521 }
57131575 2522 }
fbddd006 2523 AliDebug(2,"Track matched to closest cluster");
57131575 2524}
2525
88b96ad8 2526//_________________________________________________________________________
dda65b42 2527void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
57131575 2528{
2529 // Checks if EMC clusters are matched to ESD track.
2530 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2531
7f5392da 2532 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
2533 {
dda65b42 2534 AliVCluster *cluster = event->GetCaloCluster(iClus);
57131575 2535 if (!cluster->IsEMCAL())
2536 continue;
2537
2538 Int_t nTracks = event->GetNumberOfTracks();
2539 TArrayI arrayTrackMatched(nTracks);
2540
2541 // Get the closest track matched to the cluster
2542 Int_t nMatched = 0;
2543 Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
7f5392da 2544 if (matchTrackIndex != -1)
2545 {
57131575 2546 arrayTrackMatched[nMatched] = matchTrackIndex;
2547 nMatched++;
2548 }
2549
2550 // Get all other tracks matched to the cluster
8517fbd8 2551 for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
7f5392da 2552 {
dda65b42 2553 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
0ddff478 2554
2555 if( !track ) continue;
2556
2557 if ( iTrk == matchTrackIndex ) continue;
2558
2559 if ( track->GetEMCALcluster() == iClus )
2560 {
57131575 2561 arrayTrackMatched[nMatched] = iTrk;
2562 ++nMatched;
2563 }
2564 }
2565
2566 //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2567
2568 arrayTrackMatched.Set(nMatched);
dda65b42 2569 AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2570 if (esdcluster)
2571 esdcluster->AddTracksMatched(arrayTrackMatched);
2572 else if (nMatched>0) {
2573 AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2574 if (aodcluster)
2575 aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2576 }
57131575 2577
2578 Float_t eta= -999, phi = -999;
2579 if (matchTrackIndex != -1)
2580 GetMatchedResiduals(iClus, eta, phi);
2581 cluster->SetTrackDistance(phi, eta);
2582 }
2583
fbddd006 2584 AliDebug(2,"Cluster matched to tracks");
57131575 2585}
2586
b540d03f 2587//___________________________________________________
d9b3567c 2588void AliEMCALRecoUtils::Print(const Option_t *) const
2589{
2590 // Print Parameters
2591
2592 printf("AliEMCALRecoUtils Settings: \n");
2593 printf("Misalignment shifts\n");
8517fbd8 2594 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 2595 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2596 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 2597 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
8517fbd8 2598 for (Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 2599
2600 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 2601
fa4287a2 2602 printf("Matching criteria: ");
8517fbd8 2603 if (fCutEtaPhiSum) {
2604 printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
2605 } else if (fCutEtaPhiSeparate) {
2606 printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
2607 } else {
2608 printf("Error\n");
2609 printf("please specify your cut criteria\n");
2610 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2611 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2612 }
fa4287a2 2613
8fc351e3 2614 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);
2615 printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
bd8c7aef 2616
2617 printf("Track cuts: \n");
fa4287a2 2618 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
a6a1e3ab 2619 printf("AOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
bd8c7aef 2620 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2621 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2622 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2623 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2624 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
d9b3567c 2625}