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