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