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