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