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