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