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