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