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