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