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