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