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