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