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