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