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