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