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