]> 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
For AliAnalysisTaskSingleMu: remove check for IsPileupSPD; possibility to add the...
[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
16/* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
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//_______________________________________________________________________________
335Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
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,
464 const 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//___________________________________________________________________________
495Float_t AliEMCALRecoUtils::GetECross(const Int_t absID, const Double_t tcell,
496 AliVCaloCells* cells, const Int_t bc)
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//_____________________________________________________________________________________________
571Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
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,
601 const 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//_________________________________________________________
876Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy,
877 const Int_t iParticle,
878 const 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,
1178 const Int_t bc)
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,
1251 const 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//_______________________________________________________________________________________________________
1294void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const 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 {
88b96ad8 1832 AliInfo("Init the magnetic field\n");
1833 if (esdevent)
1834 {
1835 esdevent->InitMagneticField();
1836 }
1837 else if(aodevent)
1838 {
1839 Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
1840 Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
1841 AliMagF *field = AliMagF::CreateFieldMap(curSol,curDip);
1842 TGeoGlobalMagField::Instance()->SetField(field);
1843 }
1844 else
1845 {
1846 AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1847 }
1848
1849 } // Init mag field
1850
42ceff04 1851 if (esdevent) {
1852 UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
1853 UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
1854 Bool_t desc1 = (mask1 >> 3) & 0x1;
1855 Bool_t desc2 = (mask2 >> 3) & 0x1;
1856 if (desc1==0 || desc2==0) {
30e29b2a 1857// AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
1858// mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
1859// mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
42ceff04 1860 fITSTrackSA=kTRUE;
1861 }
1862 }
1863
8fc351e3 1864 TObjArray *clusterArray = 0x0;
1865 if(!clusterArr)
1866 {
1867 clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
1868 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
fbddd006 1869 {
1870 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1871 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1872 clusterArray->AddAt(cluster,icl);
1873 }
8fc351e3 1874 }
61160f1f 1875
bd8c7aef 1876 Int_t matched=0;
bb6f5f0b 1877 Double_t cv[21];
1878 for (Int_t i=0; i<21;i++) cv[i]=0;
bd8c7aef 1879 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1880 {
456126ad 1881 AliExternalTrackParam *trackParam = 0;
61160f1f 1882
bb6f5f0b 1883 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
8fc351e3 1884 AliESDtrack *esdTrack = 0;
1885 AliAODTrack *aodTrack = 0;
1c7a2bf4 1886 if(esdevent)
61160f1f 1887 {
8fc351e3 1888 esdTrack = esdevent->GetTrack(itr);
1889 if(!esdTrack) continue;
1890 if(!IsAccepted(esdTrack)) continue;
61160f1f 1891 if(esdTrack->Pt()<fCutMinTrackPt) continue;
8fc351e3 1892 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
1893 if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
42ceff04 1894 if(!fITSTrackSA)
1895 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
1896 else
1897 trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
61160f1f 1898 }
bb6f5f0b 1899
1900 //If the input event is AOD, the starting point for extrapolation is at vertex
8fc351e3 1901 //AOD tracks are selected according to its filterbit.
1c7a2bf4 1902 else if(aodevent)
61160f1f 1903 {
8fc351e3 1904 aodTrack = aodevent->GetTrack(itr);
61160f1f 1905 if(!aodTrack) continue;
a6a1e3ab 1906
1907 if(fAODTPCOnlyTracks) // Match with TPC only tracks, default from May 2013, before filter bit 32
1908 {
1909 //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
1910 if(!aodTrack->IsTPCOnly()) continue ;
1911 }
1912 else if(fAODHybridTracks) // Match with hybrid tracks
1913 {
1914 //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
1915 if(!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
1916 }
1917 else // Match with tracks on a mask
1918 {
1919 //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
1920 if(!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
1921 }
1af378e6 1922
61160f1f 1923 if(aodTrack->Pt()<fCutMinTrackPt) continue;
1af378e6 1924
8fc351e3 1925 Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
1926 if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
61160f1f 1927 Double_t pos[3],mom[3];
1928 aodTrack->GetXYZ(pos);
1929 aodTrack->GetPxPyPz(mom);
1930 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 1931
61160f1f 1932 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1933 }
bd8c7aef 1934
bb6f5f0b 1935 //Return if the input data is not "AOD" or "ESD"
1936 else
61160f1f 1937 {
1938 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
8fc351e3 1939 if(clusterArray)
fbddd006 1940 {
1941 clusterArray->Clear();
1942 delete clusterArray;
1943 }
61160f1f 1944 return;
1945 }
1946
bb6f5f0b 1947 if(!trackParam) continue;
8fc351e3 1948
1949 //Extrapolate the track to EMCal surface
1950 AliExternalTrackParam emcalParam(*trackParam);
a29b2a8a 1951 Float_t eta, phi, pt;
1952 if(!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
8fc351e3 1953 {
fbddd006 1954 if(aodevent && trackParam) delete trackParam;
30e29b2a 1955 if(fITSTrackSA && trackParam) delete trackParam;
fbddd006 1956 continue;
8fc351e3 1957 }
1958
150f4870 1959// if(esdevent)
1960// {
fbddd006 1961// esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
150f4870 1962// }
8fc351e3 1963
1964 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
1965 {
fbddd006 1966 if(aodevent && trackParam) delete trackParam;
30e29b2a 1967 if(fITSTrackSA && trackParam) delete trackParam;
fbddd006 1968 continue;
8fc351e3 1969 }
1970
1971
1972 //Find matched clusters
bd8c7aef 1973 Int_t index = -1;
8fc351e3 1974 Float_t dEta = -999, dPhi = -999;
1975 if(!clusterArr)
61160f1f 1976 {
fbddd006 1977 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
8fc351e3 1978 }
1979 else
61160f1f 1980 {
fbddd006 1981 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
8fc351e3 1982 }
61160f1f 1983
bd8c7aef 1984 if(index>-1)
1985 {
8fc351e3 1986 fMatchedTrackIndex ->AddAt(itr,matched);
1987 fMatchedClusterIndex ->AddAt(index,matched);
1988 fResidualEta ->AddAt(dEta,matched);
1989 fResidualPhi ->AddAt(dPhi,matched);
bd8c7aef 1990 matched++;
1991 }
456126ad 1992 if(aodevent && trackParam) delete trackParam;
30e29b2a 1993 if(fITSTrackSA && trackParam) delete trackParam;
bd8c7aef 1994 }//track loop
8fc351e3 1995
1996 if(clusterArray)
1997 {
1998 clusterArray->Clear();
1999 delete clusterArray;
2000 }
b540d03f 2001
2002 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2003
8fc351e3 2004 fMatchedTrackIndex ->Set(matched);
2005 fMatchedClusterIndex ->Set(matched);
2006 fResidualPhi ->Set(matched);
2007 fResidualEta ->Set(matched);
bd8c7aef 2008}
2009
b540d03f 2010//________________________________________________________________________________
a520bcd0 2011Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
2012 const AliVEvent *event,
2013 const AliEMCALGeometry *geom,
2014 Float_t &dEta, Float_t &dPhi)
bb6f5f0b 2015{
2016 //
2017 // This function returns the index of matched cluster to input track
fa4287a2 2018 // Returns -1 if no match is found
bb6f5f0b 2019 Int_t index = -1;
8fc351e3 2020 Double_t phiV = track->Phi()*TMath::RadToDeg();
2021 if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
42ceff04 2022 AliExternalTrackParam *trackParam = 0;
2023 if(!fITSTrackSA)
2024 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
2025 else
2026 trackParam = new AliExternalTrackParam(*track);
2027
8fc351e3 2028 if(!trackParam) return index;
2029 AliExternalTrackParam emcalParam(*trackParam);
a29b2a8a 2030 Float_t eta, phi, pt;
8fc351e3 2031
a29b2a8a 2032 if(!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) {
30e29b2a 2033 if(fITSTrackSA) delete trackParam;
2034 return index;
2035 }
2036 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){
2037 if(fITSTrackSA) delete trackParam;
2038 return index;
2039 }
2040
8fc351e3 2041 TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2042
bb6f5f0b 2043 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 2044 {
2045 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8fc351e3 2046 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2047 clusterArr->AddAt(cluster,icl);
2048 }
2049
2050 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2051 clusterArr->Clear();
2052 delete clusterArr;
30e29b2a 2053 if(fITSTrackSA) delete trackParam;
2054
8fc351e3 2055 return index;
2056}
2057
7f5392da 2058//_______________________________________________________________________________________________
2059Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
a520bcd0 2060 AliExternalTrackParam *trkParam,
7f5392da 2061 const TObjArray * clusterArr,
a520bcd0 2062 Float_t &dEta, Float_t &dPhi)
8fc351e3 2063{
7f5392da 2064 // Find matched cluster in array
2065
8fc351e3 2066 dEta=-999, dPhi=-999;
2067 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2068 Int_t index = -1;
ee602376 2069 Float_t tmpEta=-999, tmpPhi=-999;
8fc351e3 2070
2071 Double_t exPos[3] = {0.,0.,0.};
2072 if(!emcalParam->GetXYZ(exPos)) return index;
2073
2074 Float_t clsPos[3] = {0.,0.,0.};
2075 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
bb6f5f0b 2076 {
8fc351e3 2077 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2078 if(!cluster || !cluster->IsEMCAL()) continue;
2079 cluster->GetPosition(clsPos);
2080 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));
2081 if(dR > fClusterWindow) continue;
2082
2083 AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
ee602376 2084 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
8fc351e3 2085 if(fCutEtaPhiSum)
2086 {
2087 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2088 if(tmpR<dRMax)
fbddd006 2089 {
2090 dRMax=tmpR;
2091 dEtaMax=tmpEta;
2092 dPhiMax=tmpPhi;
2093 index=icl;
2094 }
8fc351e3 2095 }
2096 else if(fCutEtaPhiSeparate)
2097 {
2098 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
fbddd006 2099 {
2100 dEtaMax = tmpEta;
2101 dPhiMax = tmpPhi;
2102 index=icl;
2103 }
8fc351e3 2104 }
2105 else
2106 {
2107 printf("Error: please specify your cut criteria\n");
2108 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2109 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2110 return index;
2111 }
61160f1f 2112 }
8fc351e3 2113
2114 dEta=dEtaMax;
2115 dPhi=dPhiMax;
2116
bb6f5f0b 2117 return index;
2118}
2119
88b96ad8 2120//------------------------------------------------------------------------------------
2121Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
a520bcd0 2122 const Double_t emcalR,
2123 const Double_t mass,
2124 const Double_t step,
88b96ad8 2125 Float_t &eta,
a29b2a8a 2126 Float_t &phi,
2127 Float_t &pt)
ee602376 2128{
88b96ad8 2129 //Extrapolate track to EMCAL surface
2130
a29b2a8a 2131 eta = -999, phi = -999, pt = -999;
ee602376 2132 if(!trkParam) return kFALSE;
2133 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
2134 Double_t trkPos[3] = {0.,0.,0.};
2135 if(!trkParam->GetXYZ(trkPos)) return kFALSE;
2136 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2137 eta = trkPosVec.Eta();
2138 phi = trkPosVec.Phi();
a29b2a8a 2139 pt = trkParam->Pt();
ee602376 2140 if(phi<0)
2141 phi += 2*TMath::Pi();
2142
2143 return kTRUE;
2144}
2145
88b96ad8 2146//-----------------------------------------------------------------------------------
2147Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2148 const Float_t *clsPos,
2149 Double_t mass,
2150 Double_t step,
2151 Float_t &tmpEta,
2152 Float_t &tmpPhi)
ee602376 2153{
2154 //
2155 //Return the residual by extrapolating a track param to a global position
2156 //
2157 tmpEta = -999;
2158 tmpPhi = -999;
2159 if(!trkParam) return kFALSE;
2160 Double_t trkPos[3] = {0.,0.,0.};
2161 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2162 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2163 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
2164 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
2165 if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
2166
2167 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2168 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2169
2170 // track cluster matching
2171 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
2172 tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2173
2174 return kTRUE;
2175}
2176
88b96ad8 2177//----------------------------------------------------------------------------------
2178Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2179 const AliVCluster *cluster,
a520bcd0 2180 const Double_t mass,
2181 const Double_t step,
88b96ad8 2182 Float_t &tmpEta,
2183 Float_t &tmpPhi)
ee602376 2184{
2185 //
2186 //Return the residual by extrapolating a track param to a cluster
2187 //
2188 tmpEta = -999;
2189 tmpPhi = -999;
2190 if(!cluster || !trkParam) return kFALSE;
2191
2192 Float_t clsPos[3] = {0.,0.,0.};
2193 cluster->GetPosition(clsPos);
2194
2195 return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2196}
2197
88b96ad8 2198//---------------------------------------------------------------------------------
2199Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2200 const AliVCluster *cluster,
88b96ad8 2201 Float_t &tmpEta,
2202 Float_t &tmpPhi)
bb6f5f0b 2203{
2204 //
ee602376 2205 //Return the residual by extrapolating a track param to a clusterfStepCluster
bb6f5f0b 2206 //
8fc351e3 2207
ee602376 2208 return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
bb6f5f0b 2209}
2210
a520bcd0 2211//_______________________________________________________________________
2212void AliEMCALRecoUtils::GetMatchedResiduals(const Int_t clsIndex,
2213 Float_t &dEta, Float_t &dPhi)
bd8c7aef 2214{
bb6f5f0b 2215 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 2216 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 2217 //Works with ESDs and AODs
bd8c7aef 2218
bb6f5f0b 2219 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 2220 {
2221 AliDebug(2,"No matched tracks found!\n");
fa4287a2 2222 dEta=999.;
2223 dPhi=999.;
bd8c7aef 2224 return;
2225 }
fa4287a2 2226 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2227 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 2228}
841dbf60 2229
88b96ad8 2230//______________________________________________________________________________________________
fa4287a2 2231void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 2232{
2233 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 2234 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 2235 //Works with ESDs and AODs
2236
2237 if( FindMatchedPosForTrack(trkIndex) >= 999 )
2238 {
2239 AliDebug(2,"No matched cluster found!\n");
fa4287a2 2240 dEta=999.;
2241 dPhi=999.;
bb6f5f0b 2242 return;
2243 }
fa4287a2 2244 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2245 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 2246}
2247
2248//__________________________________________________________
2249Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2250{
2251 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2252 //Get the index of matched track to this cluster
2253 //Works with ESDs and AODs
2254
2255 if(IsClusterMatched(clsIndex))
2256 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2257 else
2258 return -1;
bd8c7aef 2259}
2260
b540d03f 2261//__________________________________________________________
bb6f5f0b 2262Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 2263{
bb6f5f0b 2264 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2265 //Get the index of matched cluster to this track
2266 //Works with ESDs and AODs
b540d03f 2267
bb6f5f0b 2268 if(IsTrackMatched(trkIndex))
2269 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 2270 else
2271 return -1;
2272}
2273
7f5392da 2274//______________________________________________________________
7cdec71f 2275Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 2276{
2277 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2278 //Returns if the cluster has a match
2279 if(FindMatchedPosForCluster(clsIndex) < 999)
2280 return kTRUE;
2281 else
2282 return kFALSE;
2283}
b540d03f 2284
7f5392da 2285//____________________________________________________________
7cdec71f 2286Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 2287{
bb6f5f0b 2288 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2289 //Returns if the track has a match
2290 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 2291 return kTRUE;
bd8c7aef 2292 else
2293 return kFALSE;
2294}
bb6f5f0b 2295
7f5392da 2296//______________________________________________________________________
bb6f5f0b 2297UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 2298{
bb6f5f0b 2299 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 2300 //Returns the position of the match in the fMatchedClusterIndex array
2301 Float_t tmpR = fCutR;
81efb149 2302 UInt_t pos = 999;
b540d03f 2303
7f5392da 2304 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
2305 {
2306 if(fMatchedClusterIndex->At(i)==clsIndex)
2307 {
841dbf60 2308 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
7f5392da 2309 if(r<tmpR)
2310 {
841dbf60 2311 pos=i;
2312 tmpR=r;
7f5392da 2313 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2314 fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2315 }
841dbf60 2316 }
bb6f5f0b 2317 }
2318 return pos;
2319}
2320
7f5392da 2321//____________________________________________________________________
bb6f5f0b 2322UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2323{
2324 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2325 //Returns the position of the match in the fMatchedTrackIndex array
2326 Float_t tmpR = fCutR;
2327 UInt_t pos = 999;
2328
7f5392da 2329 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
2330 {
2331 if(fMatchedTrackIndex->At(i)==trkIndex)
2332 {
841dbf60 2333 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
7f5392da 2334 if(r<tmpR)
2335 {
841dbf60 2336 pos=i;
2337 tmpR=r;
7f5392da 2338 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2339 fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2340 }
841dbf60 2341 }
b540d03f 2342 }
bd8c7aef 2343 return pos;
2344}
2345
a520bcd0 2346//__________________________________________________________________________
2347Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
2348 const AliEMCALGeometry *geom,
a7e5a381 2349 AliVCaloCells* cells,const Int_t bc)
b5078f5d 2350{
2351 // check if the cluster survives some quality cut
2352 //
2353 //
2354 Bool_t isGood=kTRUE;
7f5392da 2355
a7e5a381 2356 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
2357
fa4287a2 2358 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
a7e5a381 2359
fa4287a2 2360 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
a7e5a381 2361
2362 if(IsExoticCluster(cluster, cells,bc)) return kFALSE;
b5078f5d 2363
2364 return isGood;
2365}
2366
b540d03f 2367//__________________________________________________________
bd8c7aef 2368Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2369{
2370 // Given a esd track, return whether the track survive all the cuts
2371
2372 // The different quality parameter are first
2373 // retrieved from the track. then it is found out what cuts the
2374 // track did not survive and finally the cuts are imposed.
2375
2376 UInt_t status = esdTrack->GetStatus();
2377
2378 Int_t nClustersITS = esdTrack->GetITSclusters(0);
2379 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2380
2381 Float_t chi2PerClusterITS = -1;
2382 Float_t chi2PerClusterTPC = -1;
2383 if (nClustersITS!=0)
2384 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2385 if (nClustersTPC!=0)
2386 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 2387
2388
2389 //DCA cuts
827f9f23 2390 if(fTrackCutsType==kGlobalCut)
2391 {
2392 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2393 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2394 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2395 }
82d09e74 2396
2397
bd8c7aef 2398 Float_t b[2];
2399 Float_t bCov[3];
2400 esdTrack->GetImpactParameters(b,bCov);
7f5392da 2401 if (bCov[0]<=0 || bCov[2]<=0)
2402 {
bd8c7aef 2403 AliDebug(1, "Estimated b resolution lower or equal zero!");
2404 bCov[0]=0; bCov[2]=0;
2405 }
2406
2407 Float_t dcaToVertexXY = b[0];
2408 Float_t dcaToVertexZ = b[1];
2409 Float_t dcaToVertex = -1;
2410
2411 if (fCutDCAToVertex2D)
2412 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2413 else
2414 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2415
2416 // cut the track?
2417
2418 Bool_t cuts[kNCuts];
2419 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2420
2421 // track quality cuts
2422 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2423 cuts[0]=kTRUE;
2424 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2425 cuts[1]=kTRUE;
2426 if (nClustersTPC<fCutMinNClusterTPC)
2427 cuts[2]=kTRUE;
2428 if (nClustersITS<fCutMinNClusterITS)
2429 cuts[3]=kTRUE;
2430 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
2431 cuts[4]=kTRUE;
2432 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
2433 cuts[5]=kTRUE;
2434 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2435 cuts[6]=kTRUE;
2436 if (fCutDCAToVertex2D && dcaToVertex > 1)
2437 cuts[7] = kTRUE;
2438 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2439 cuts[8] = kTRUE;
2440 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2441 cuts[9] = kTRUE;
2442
827f9f23 2443 if(fTrackCutsType==kGlobalCut)
2444 {
2445 //Require at least one SPD point + anything else in ITS
2446 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
fbddd006 2447 cuts[10] = kTRUE;
827f9f23 2448 }
82d09e74 2449
42ceff04 2450 // ITS
2451 if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
2452 if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
2453 // TPC tracks
2454 cuts[11] = kTRUE;
2455 }else{
2456 // ITS standalone tracks
2457 if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
2458 if(status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
2459 }else if(fCutRequireITSpureSA){
2460 if(!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
2461 }
2462 }
2463 }
2464
bd8c7aef 2465 Bool_t cut=kFALSE;
827f9f23 2466 for (Int_t i=0; i<kNCuts; i++)
7f5392da 2467 if (cuts[i]) { cut = kTRUE ; }
bd8c7aef 2468
2469 // cut the track
2470 if (cut)
2471 return kFALSE;
2472 else
2473 return kTRUE;
2474}
827f9f23 2475
88b96ad8 2476//_____________________________________
bd8c7aef 2477void AliEMCALRecoUtils::InitTrackCuts()
2478{
2479 //Intilize the track cut criteria
5f7714ad 2480 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
bd8c7aef 2481 //Also you can customize the cuts using the setters
82d09e74 2482
5f7714ad 2483 switch (fTrackCutsType)
88b96ad8 2484 {
5f7714ad 2485 case kTPCOnlyCut:
88b96ad8 2486 {
2487 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2488 //TPC
2489 SetMinNClustersTPC(70);
2490 SetMaxChi2PerClusterTPC(4);
2491 SetAcceptKinkDaughters(kFALSE);
2492 SetRequireTPCRefit(kFALSE);
2493
2494 //ITS
2495 SetRequireITSRefit(kFALSE);
2496 SetMaxDCAToVertexZ(3.2);
2497 SetMaxDCAToVertexXY(2.4);
2498 SetDCAToVertex2D(kTRUE);
2499
2500 break;
2501 }
2502
5f7714ad 2503 case kGlobalCut:
88b96ad8 2504 {
2505 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2506 //TPC
2507 SetMinNClustersTPC(70);
2508 SetMaxChi2PerClusterTPC(4);
2509 SetAcceptKinkDaughters(kFALSE);
2510 SetRequireTPCRefit(kTRUE);
2511
2512 //ITS
2513 SetRequireITSRefit(kTRUE);
2514 SetMaxDCAToVertexZ(2);
2515 SetMaxDCAToVertexXY();
2516 SetDCAToVertex2D(kFALSE);
2517
2518 break;
2519 }
2520
0e7de35b 2521 case kLooseCut:
88b96ad8 2522 {
2523 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2524 SetMinNClustersTPC(50);
2525 SetAcceptKinkDaughters(kTRUE);
2526
2527 break;
5f7714ad 2528 }
42ceff04 2529
2530 case kITSStandAlone:
2531 {
2532 AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
2533 SetRequireITSRefit(kTRUE);
2534 SetRequireITSStandAlone(kTRUE);
2535 SetITSTrackSA(kTRUE);
2536 break;
2537 }
2538
88b96ad8 2539 }
bd8c7aef 2540}
83bfd77a 2541
57131575 2542
88b96ad8 2543//________________________________________________________________________
dda65b42 2544void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
57131575 2545{
2546 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
dda65b42 2547
57131575 2548 Int_t nTracks = event->GetNumberOfTracks();
7f5392da 2549 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
2550 {
dda65b42 2551 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
7f5392da 2552 if (!track)
2553 {
57131575 2554 AliWarning(Form("Could not receive track %d", iTrack));
2555 continue;
2556 }
7f5392da 2557
fbddd006 2558 Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
57131575 2559 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
dda65b42 2560 /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2561 AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2562 if (esdtrack) {
2563 if(matchClusIndex != -1)
2564 esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2565 else
2566 esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2567 } else {
2568 AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
2569 if(matchClusIndex != -1)
2570 aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2571 else
2572 aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2573 }
2574
57131575 2575 }
fbddd006 2576 AliDebug(2,"Track matched to closest cluster");
57131575 2577}
2578
88b96ad8 2579//_________________________________________________________________________
dda65b42 2580void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
57131575 2581{
2582 // Checks if EMC clusters are matched to ESD track.
2583 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2584
7f5392da 2585 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
2586 {
dda65b42 2587 AliVCluster *cluster = event->GetCaloCluster(iClus);
57131575 2588 if (!cluster->IsEMCAL())
2589 continue;
2590
2591 Int_t nTracks = event->GetNumberOfTracks();
2592 TArrayI arrayTrackMatched(nTracks);
2593
2594 // Get the closest track matched to the cluster
2595 Int_t nMatched = 0;
2596 Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
7f5392da 2597 if (matchTrackIndex != -1)
2598 {
57131575 2599 arrayTrackMatched[nMatched] = matchTrackIndex;
2600 nMatched++;
2601 }
2602
2603 // Get all other tracks matched to the cluster
7f5392da 2604 for(Int_t iTrk=0; iTrk<nTracks; ++iTrk)
2605 {
dda65b42 2606 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
57131575 2607 if(iTrk == matchTrackIndex) continue;
7f5392da 2608 if(track->GetEMCALcluster() == iClus)
2609 {
57131575 2610 arrayTrackMatched[nMatched] = iTrk;
2611 ++nMatched;
2612 }
2613 }
2614
2615 //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2616
2617 arrayTrackMatched.Set(nMatched);
dda65b42 2618 AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2619 if (esdcluster)
2620 esdcluster->AddTracksMatched(arrayTrackMatched);
2621 else if (nMatched>0) {
2622 AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2623 if (aodcluster)
2624 aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2625 }
57131575 2626
2627 Float_t eta= -999, phi = -999;
2628 if (matchTrackIndex != -1)
2629 GetMatchedResiduals(iClus, eta, phi);
2630 cluster->SetTrackDistance(phi, eta);
2631 }
2632
fbddd006 2633 AliDebug(2,"Cluster matched to tracks");
57131575 2634}
2635
b540d03f 2636//___________________________________________________
d9b3567c 2637void AliEMCALRecoUtils::Print(const Option_t *) const
2638{
2639 // Print Parameters
2640
2641 printf("AliEMCALRecoUtils Settings: \n");
2642 printf("Misalignment shifts\n");
2a71e873 2643 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,
2644 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2645 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 2646 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
2647 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 2648
2649 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 2650
fa4287a2 2651 printf("Matching criteria: ");
2652 if(fCutEtaPhiSum)
2653 {
8fc351e3 2654 printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
fa4287a2 2655 }
2656 else if(fCutEtaPhiSeparate)
2657 {
8fc351e3 2658 printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
fa4287a2 2659 }
2660 else
2661 {
2662 printf("Error\n");
2663 printf("please specify your cut criteria\n");
2664 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2665 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2666 }
2667
8fc351e3 2668 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);
2669 printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
bd8c7aef 2670
2671 printf("Track cuts: \n");
fa4287a2 2672 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
a6a1e3ab 2673 printf("AOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
bd8c7aef 2674 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2675 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2676 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2677 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2678 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
d9b3567c 2679}
96957075 2680