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