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