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