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