]> 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
bug fix - update to track cuts by Rongrong
[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
094786cc 286//__________________________________________________
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
01d44f1f 432//__________________________________________________
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
d9b3567c 455//__________________________________________________
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}
7e0ecb89 540//__________________________________________________
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
d9b3567c 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
fd6df01c 760//________________________________________________________________
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
094786cc 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
888 if(!IsRecalibrationOn()) return;
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
3bfc4732 929//________________________________________________________________
930void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime){
931 // Recalibrate time of cell with absID considering the recalibration map
932 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
933
934 if(!fCellsRecalibrated && IsTimeRecalibrationOn()){
935// printf("cell time org %g, ",celltime);
094786cc 936
3bfc4732 937 Double_t timeBCoffset = 0.;
938 if( bc%4 ==0 || bc%4==1) timeBCoffset = 100.*1.e-9; //in ns
939
940 Double_t celloffset = GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;
941
942// printf("absId %d, time %f bc %d-%d: bc0 %f, bc1 %f, bc2 %f, bc3 %f \n", absId, celltime*1.e9,bc, bc%4,
943// GetEMCALChannelTimeRecalibrationFactor(0,absId),GetEMCALChannelTimeRecalibrationFactor(1,absId),
944// GetEMCALChannelTimeRecalibrationFactor(2,absId),GetEMCALChannelTimeRecalibrationFactor(3,absId));
945
946 celltime -= timeBCoffset ;
947 celltime -= celloffset ;
948// printf("new %g\n",celltime);
949 }
950
951}
952
d9b3567c 953//__________________________________________________
094786cc 954void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 955{
956 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
957
2aeb4226 958 if(!clu){
959 AliInfo("Cluster pointer null!");
960 return;
961 }
962
094786cc 963 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
964 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 965 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 966
967}
968
969//__________________________________________________
970void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
971{
972 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
973 // The algorithm is a copy of what is done in AliEMCALRecPoint
974
975 Double_t eCell = 0.;
976 Float_t fraction = 1.;
977 Float_t recalFactor = 1.;
978
979 Int_t absId = -1;
980 Int_t iTower = -1, iIphi = -1, iIeta = -1;
981 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
982 Float_t weight = 0., totalWeight=0.;
983 Float_t newPos[3] = {0,0,0};
984 Double_t pLocal[3], pGlobal[3];
cb231979 985 Bool_t shared = kFALSE;
986
094786cc 987 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
cb231979 988 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 989 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
990
83bfd77a 991 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
992
094786cc 993 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 994
3bfc4732 995 absId = clu->GetCellAbsId(iDig);
996 fraction = clu->GetCellAmplitudeFraction(iDig);
997 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
998
999 if(!fCellsRecalibrated){
1000
1001 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1002 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1003
1004 if(IsRecalibrationOn()) {
1005 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1006 }
094786cc 1007 }
3bfc4732 1008
094786cc 1009 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1010
1011 weight = GetCellWeight(eCell,clEnergy);
1012 totalWeight += weight;
3bfc4732 1013
094786cc 1014 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 1015 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 1016 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 1017 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1018
094786cc 1019 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1020
1021 }// cell loop
1022
1023 if(totalWeight>0){
1024 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1025 }
1026
094786cc 1027 //Float_t pos[]={0,0,0};
1028 //clu->GetPosition(pos);
1029 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 1030 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 1031
1032 if(iSupModMax > 1) {//sector 1
1033 newPos[0] +=fMisalTransShift[3];//-=3.093;
1034 newPos[1] +=fMisalTransShift[4];//+=6.82;
1035 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 1036 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1037
094786cc 1038 }
1039 else {//sector 0
1040 newPos[0] +=fMisalTransShift[0];//+=1.134;
1041 newPos[1] +=fMisalTransShift[1];//+=8.2;
1042 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 1043 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1044
094786cc 1045 }
83bfd77a 1046 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1047
094786cc 1048 clu->SetPosition(newPos);
1049
094786cc 1050}
1051
1052//__________________________________________________
1053void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
1054{
1055 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1056 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1057
1058 Double_t eCell = 1.;
1059 Float_t fraction = 1.;
1060 Float_t recalFactor = 1.;
1061
1062 Int_t absId = -1;
d9b3567c 1063 Int_t iTower = -1;
094786cc 1064 Int_t iIphi = -1, iIeta = -1;
1065 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 1066 Int_t iphi = -1, ieta =-1;
cb231979 1067 Bool_t shared = kFALSE;
1068
d9b3567c 1069 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
cb231979 1070 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1071 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1072
d9b3567c 1073 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 1074 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1075 Int_t startingSM = -1;
d9b3567c 1076
1077 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 1078 absId = clu->GetCellAbsId(iDig);
1079 fraction = clu->GetCellAmplitudeFraction(iDig);
1080 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1081
d9b3567c 1082 if (iDig==0) startingSM = iSupMod;
1083 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 1084
1085 eCell = cells->GetCellAmplitude(absId);
d9b3567c 1086
3bfc4732 1087 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1088 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1089
1090 if(!fCellsRecalibrated){
1091
1092 if(IsRecalibrationOn()) {
1093
1094 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1095
1096 }
094786cc 1097 }
3bfc4732 1098
094786cc 1099 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 1100
094786cc 1101 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 1102 if(weight < 0) weight = 0;
1103 totalWeight += weight;
1104 weightedCol += ieta*weight;
1105 weightedRow += iphi*weight;
1106
1107 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1108
094786cc 1109 }// cell loop
1110
d9b3567c 1111 Float_t xyzNew[]={0.,0.,0.};
1112 if(areInSameSM == kTRUE) {
1113 //printf("In Same SM\n");
1114 weightedCol = weightedCol/totalWeight;
1115 weightedRow = weightedRow/totalWeight;
094786cc 1116 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 1117 }
1118 else {
1119 //printf("In Different SM\n");
094786cc 1120 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 1121 }
d9b3567c 1122
094786cc 1123 clu->SetPosition(xyzNew);
d9b3567c 1124
1125}
1126
cb231979 1127//____________________________________________________________________________
1128void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
1129
1130 //re-evaluate distance to bad channel with updated bad map
1131
78467229 1132 if(!fRecalDistToBadChannels) return;
cb231979 1133
2aeb4226 1134 if(!cluster){
1135 AliInfo("Cluster pointer null!");
1136 return;
1137 }
1138
cb231979 1139 //Get channels map of the supermodule where the cluster is.
cb231979 1140 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
1141 Bool_t shared = kFALSE;
1142 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1143 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1144
1145 Int_t dRrow, dRcol;
1146 Float_t minDist = 10000.;
1147 Float_t dist = 0.;
1148
1149 //Loop on tower status map
1150 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1151 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1152 //Check if tower is bad.
1153 if(hMap->GetBinContent(icol,irow)==0) continue;
1154 //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 1155 // iSupMod,icol, irow, icolM,irowM);
cb231979 1156
1157 dRrow=TMath::Abs(irowM-irow);
1158 dRcol=TMath::Abs(icolM-icol);
1159 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1160 if(dist < minDist){
1161 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1162 minDist = dist;
1163 }
1164
1165 }
1166 }
1167
1168 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
1169 if (shared) {
1170 TH2D* hMap2 = 0;
1171 Int_t iSupMod2 = -1;
1172
1173 //The only possible combinations are (0,1), (2,3) ... (8,9)
1174 if(iSupMod%2) iSupMod2 = iSupMod-1;
1175 else iSupMod2 = iSupMod+1;
1176 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
1177
1178 //Loop on tower status map of second super module
1179 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1180 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1181 //Check if tower is bad.
1182 if(hMap2->GetBinContent(icol,irow)==0) continue;
1183 //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",
1184 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
1185
1186 dRrow=TMath::Abs(irow-irowM);
1187
1188 if(iSupMod%2) {
1189 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1190 }
1191 else {
1192 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1193 }
1194
1195 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1196 if(dist < minDist) minDist = dist;
1197
1198 }
1199 }
1200
1201 }// shared cluster in 2 SuperModules
78467229 1202
6fe0e6d0 1203 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1204 cluster->SetDistanceToBadChannel(minDist);
cb231979 1205
1206}
1207
83bfd77a 1208//____________________________________________________________________________
1209void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
1210
1211 //re-evaluate identification parameters with bayesian
2aeb4226 1212
1213 if(!cluster){
1214 AliInfo("Cluster pointer null!");
1215 return;
1216 }
1217
83bfd77a 1218 if ( cluster->GetM02() != 0)
1219 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1220
1221 Float_t pidlist[AliPID::kSPECIESN+1];
1222 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1223
1224 cluster->SetPID(pidlist);
1225
1226}
1227
1228//____________________________________________________________________________
1229void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
1230{
1231 // Calculates new center of gravity in the local EMCAL-module coordinates
1232 // and tranfers into global ALICE coordinates
1233 // Calculates Dispersion and main axis
1234
2aeb4226 1235 if(!cluster){
1236 AliInfo("Cluster pointer null!");
1237 return;
1238 }
1239
83bfd77a 1240 Int_t nstat = 0;
1241 Float_t wtot = 0. ;
1242 Double_t eCell = 0.;
1243 Float_t fraction = 1.;
1244 Float_t recalFactor = 1.;
1245
1246 Int_t iSupMod = -1;
1247 Int_t iTower = -1;
1248 Int_t iIphi = -1;
1249 Int_t iIeta = -1;
1250 Int_t iphi = -1;
1251 Int_t ieta = -1;
1252 Double_t etai = -1.;
1253 Double_t phii = -1.;
1254
1255 Double_t w = 0.;
1256 Double_t d = 0.;
1257 Double_t dxx = 0.;
1258 Double_t dzz = 0.;
1259 Double_t dxz = 0.;
1260 Double_t xmean = 0.;
1261 Double_t zmean = 0.;
1262
1263 //Loop on cells
1264 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1265
1266 //Get from the absid the supermodule, tower and eta/phi numbers
1267 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1268 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1269
1270 //Get the cell energy, if recalibration is on, apply factors
1271 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1272 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1273
1274 if(!fCellsRecalibrated){
1275
1276 if(IsRecalibrationOn()) {
1277 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1278 }
1279
83bfd77a 1280 }
3bfc4732 1281
83bfd77a 1282 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1283
1284 if(cluster->E() > 0 && eCell > 0){
1285
1286 w = GetCellWeight(eCell,cluster->E());
1287
1288 etai=(Double_t)ieta;
1289 phii=(Double_t)iphi;
1290 if(w > 0.0) {
1291 wtot += w ;
1292 nstat++;
1293 //Shower shape
1294 dxx += w * etai * etai ;
1295 xmean+= w * etai ;
1296 dzz += w * phii * phii ;
1297 zmean+= w * phii ;
1298 dxz += w * etai * phii ;
1299 }
1300 }
1301 else
1302 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1303 }//cell loop
1304
1305 //Normalize to the weight
1306 if (wtot > 0) {
1307 xmean /= wtot ;
1308 zmean /= wtot ;
1309 }
1310 else
1311 AliError(Form("Wrong weight %f\n", wtot));
1312
1313 //Calculate dispersion
1314 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1315
1316 //Get from the absid the supermodule, tower and eta/phi numbers
1317 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1318 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1319
1320 //Get the cell energy, if recalibration is on, apply factors
1321 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1322 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1323 if(IsRecalibrationOn()) {
1324 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1325 }
1326 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1327
1328 if(cluster->E() > 0 && eCell > 0){
1329
1330 w = GetCellWeight(eCell,cluster->E());
1331
1332 etai=(Double_t)ieta;
1333 phii=(Double_t)iphi;
1334 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1335 }
1336 else
1337 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1338 }// cell loop
1339
1340 //Normalize to the weigth and set shower shape parameters
1341 if (wtot > 0 && nstat > 1) {
1342 d /= wtot ;
1343 dxx /= wtot ;
1344 dzz /= wtot ;
1345 dxz /= wtot ;
1346 dxx -= xmean * xmean ;
1347 dzz -= zmean * zmean ;
1348 dxz -= xmean * zmean ;
1349 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1350 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1351 }
1352 else{
1353 d=0. ;
1354 cluster->SetM20(0.) ;
1355 cluster->SetM02(0.) ;
1356 }
1357
1358 if (d>=0)
1359 cluster->SetDispersion(TMath::Sqrt(d)) ;
1360 else
1361 cluster->SetDispersion(0) ;
83bfd77a 1362}
1363
b540d03f 1364//____________________________________________________________________________
b5078f5d 1365void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
bd8c7aef 1366{
1367 //This function should be called before the cluster loop
1368 //Before call this function, please recalculate the cluster positions
1369 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1370 //Store matched cluster indexes and residuals
61160f1f 1371
fa4287a2 1372 fMatchedTrackIndex->Reset();
bd8c7aef 1373 fMatchedClusterIndex->Reset();
fa4287a2 1374 fResidualPhi->Reset();
1375 fResidualEta->Reset();
bd8c7aef 1376
fa4287a2 1377 fMatchedTrackIndex->Set(500);
67241b7e 1378 fMatchedClusterIndex->Set(500);
fa4287a2 1379 fResidualPhi->Set(500);
1380 fResidualEta->Set(500);
bd8c7aef 1381
1c7a2bf4 1382 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1383 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
61160f1f 1384
bd8c7aef 1385 Int_t matched=0;
bb6f5f0b 1386 Double_t cv[21];
1387 for (Int_t i=0; i<21;i++) cv[i]=0;
bd8c7aef 1388 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1389 {
456126ad 1390 AliExternalTrackParam *trackParam = 0;
61160f1f 1391
bb6f5f0b 1392 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1c7a2bf4 1393 if(esdevent)
61160f1f 1394 {
1395 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1396 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1397 if(esdTrack->Pt()<fCutMinTrackPt) continue;
97c0d532 1398 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
61160f1f 1399 }
bb6f5f0b 1400
1401 //If the input event is AOD, the starting point for extrapolation is at vertex
1402 //AOD tracks are selected according to its bit.
1c7a2bf4 1403 else if(aodevent)
61160f1f 1404 {
1405 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1406 if(!aodTrack) continue;
1407 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1408 if(aodTrack->Pt()<fCutMinTrackPt) continue;
1409 Double_t pos[3],mom[3];
1410 aodTrack->GetXYZ(pos);
1411 aodTrack->GetPxPyPz(mom);
1412 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()));
1413 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1414 }
bd8c7aef 1415
bb6f5f0b 1416 //Return if the input data is not "AOD" or "ESD"
1417 else
61160f1f 1418 {
1419 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1420 return;
1421 }
1422
bb6f5f0b 1423 if(!trackParam) continue;
61160f1f 1424
fa4287a2 1425 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
bd8c7aef 1426 Int_t index = -1;
b540d03f 1427 if(!clusterArr){// get clusters from event
1428 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 1429 {
1430 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1431 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1432 AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation
1433 Float_t tmpEta=-999, tmpPhi=-999;
1434 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1435 if(fCutEtaPhiSum)
1436 {
1437 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1438 if(tmpR<dRMax)
1439 {
1440 dRMax=tmpR;
1441 dEtaMax=tmpEta;
1442 dPhiMax=tmpPhi;
1443 index=icl;
1444 }
1445 }
1446 else if(fCutEtaPhiSeparate)
1447 {
1448 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1449 {
1450 dEtaMax = tmpEta;
1451 dPhiMax = tmpPhi;
1452 index=icl;
1453 }
1454 }
1455 else
1456 {
1457 printf("Error: please specify your cut criteria\n");
1458 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1459 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1460 if(aodevent && trackParam) delete trackParam;
1461 return;
1462 }
1463 }//cluster loop
1464 }
1465 else { // external cluster array, not from ESD event
b540d03f 1466 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
61160f1f 1467 {
1468 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1469 if(!cluster){
1470 AliInfo("Cluster not found!!!");
1471 continue;
1472 }
1473 if(!cluster->IsEMCAL()) continue;
1474 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1475 Float_t tmpEta=-999, tmpPhi=-999;
1476 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1477 if(fCutEtaPhiSum)
1478 {
1479 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1480 if(tmpR<dRMax)
1481 {
1482 dRMax=tmpR;
1483 dEtaMax=tmpEta;
1484 dPhiMax=tmpPhi;
1485 index=icl;
1486 }
1487 }
1488 else if(fCutEtaPhiSeparate)
1489 {
1490 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1491 {
1492 dEtaMax = tmpEta;
1493 dPhiMax = tmpPhi;
1494 index=icl;
1495 }
1496 }
1497 else
1498 {
1499 printf("Error: please specify your cut criteria\n");
1500 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1501 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1502 if(aodevent && trackParam) delete trackParam;
1503 return;
1504 }
b540d03f 1505 }//cluster loop
1506 }// external list of clusters
61160f1f 1507
bd8c7aef 1508 if(index>-1)
1509 {
b540d03f 1510 fMatchedTrackIndex ->AddAt(itr,matched);
bd8c7aef 1511 fMatchedClusterIndex->AddAt(index,matched);
fa4287a2 1512 fResidualEta ->AddAt(dEtaMax,matched);
1513 fResidualPhi ->AddAt(dPhiMax,matched);
bd8c7aef 1514 matched++;
1515 }
456126ad 1516 if(aodevent && trackParam) delete trackParam;
bd8c7aef 1517 }//track loop
b540d03f 1518
1519 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1520
1521 fMatchedTrackIndex ->Set(matched);
bd8c7aef 1522 fMatchedClusterIndex->Set(matched);
fa4287a2 1523 fResidualPhi ->Set(matched);
1524 fResidualEta ->Set(matched);
bd8c7aef 1525}
1526
b540d03f 1527//________________________________________________________________________________
b5078f5d 1528Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
bb6f5f0b 1529{
1530 //
1531 // This function returns the index of matched cluster to input track
fa4287a2 1532 // Returns -1 if no match is found
61160f1f 1533
1534
fa4287a2 1535 Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
bb6f5f0b 1536 Int_t index = -1;
61160f1f 1537
97c0d532 1538 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
61160f1f 1539
bb6f5f0b 1540 if(!trackParam) return index;
1541 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 1542 {
1543 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1544 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1545 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1546 Float_t tmpEta=-999, tmpPhi=-999;
1547 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1548 if(fCutEtaPhiSum)
bb6f5f0b 1549 {
61160f1f 1550 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1551 if(tmpR<dRMax)
fa4287a2 1552 {
1553 dRMax=tmpR;
1554 dEtaMax=tmpEta;
1555 dPhiMax=tmpPhi;
1556 index=icl;
1557 }
61160f1f 1558 }
1559 else if(fCutEtaPhiSeparate)
1560 {
1561 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
fa4287a2 1562 {
1563 dEtaMax = tmpEta;
1564 dPhiMax = tmpPhi;
1565 index=icl;
1566 }
61160f1f 1567 }
1568 else
1569 {
1570 printf("Error: please specify your cut criteria\n");
1571 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1572 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1573 return -1;
1574 }
1575 }//cluster loop
bb6f5f0b 1576 return index;
1577}
1578
1579//________________________________________________________________________________
fa4287a2 1580Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
bb6f5f0b 1581{
1582 //
1583 //Return the residual by extrapolating a track to a cluster
1584 //
1585 if(!trkParam || !cluster) return kFALSE;
1586 Float_t clsPos[3];
1587 Double_t trkPos[3];
1588 cluster->GetPosition(clsPos); //Has been recalculated
1589 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1590 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1591 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
aef4d246 1592 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.8, -1)) return kFALSE;
bb6f5f0b 1593 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
fa4287a2 1594
1595 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1596 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1597
b33358c2 1598 // track cluster matching
1599 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
fa4287a2 1600 tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
1601
bb6f5f0b 1602 return kTRUE;
1603}
1604
1605//________________________________________________________________________________
fa4287a2 1606void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
bd8c7aef 1607{
bb6f5f0b 1608 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 1609 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 1610 //Works with ESDs and AODs
bd8c7aef 1611
bb6f5f0b 1612 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 1613 {
1614 AliDebug(2,"No matched tracks found!\n");
fa4287a2 1615 dEta=999.;
1616 dPhi=999.;
bd8c7aef 1617 return;
1618 }
fa4287a2 1619 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1620 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 1621}
1622//________________________________________________________________________________
fa4287a2 1623void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1624{
1625 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 1626 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 1627 //Works with ESDs and AODs
1628
1629 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1630 {
1631 AliDebug(2,"No matched cluster found!\n");
fa4287a2 1632 dEta=999.;
1633 dPhi=999.;
bb6f5f0b 1634 return;
1635 }
fa4287a2 1636 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1637 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 1638}
1639
1640//__________________________________________________________
1641Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1642{
1643 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1644 //Get the index of matched track to this cluster
1645 //Works with ESDs and AODs
1646
1647 if(IsClusterMatched(clsIndex))
1648 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1649 else
1650 return -1;
bd8c7aef 1651}
1652
b540d03f 1653//__________________________________________________________
bb6f5f0b 1654Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 1655{
bb6f5f0b 1656 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1657 //Get the index of matched cluster to this track
1658 //Works with ESDs and AODs
b540d03f 1659
bb6f5f0b 1660 if(IsTrackMatched(trkIndex))
1661 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 1662 else
1663 return -1;
1664}
1665
bb6f5f0b 1666//__________________________________________________
7cdec71f 1667Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 1668{
1669 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1670 //Returns if the cluster has a match
1671 if(FindMatchedPosForCluster(clsIndex) < 999)
1672 return kTRUE;
1673 else
1674 return kFALSE;
1675}
b540d03f 1676
bd8c7aef 1677//__________________________________________________
7cdec71f 1678Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 1679{
bb6f5f0b 1680 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1681 //Returns if the track has a match
1682 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 1683 return kTRUE;
bd8c7aef 1684 else
1685 return kFALSE;
1686}
bb6f5f0b 1687
b540d03f 1688//__________________________________________________________
bb6f5f0b 1689UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 1690{
bb6f5f0b 1691 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 1692 //Returns the position of the match in the fMatchedClusterIndex array
1693 Float_t tmpR = fCutR;
81efb149 1694 UInt_t pos = 999;
b540d03f 1695
bd8c7aef 1696 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1697 {
fa4287a2 1698 if(fMatchedClusterIndex->At(i)==clsIndex)
1699 {
1700 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1701 if(r<tmpR)
1702 {
1703 pos=i;
1704 tmpR=r;
1705 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1706 }
1707 }
bb6f5f0b 1708 }
1709 return pos;
1710}
1711
1712//__________________________________________________________
1713UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1714{
1715 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1716 //Returns the position of the match in the fMatchedTrackIndex array
1717 Float_t tmpR = fCutR;
1718 UInt_t pos = 999;
1719
1720 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1721 {
fa4287a2 1722 if(fMatchedTrackIndex->At(i)==trkIndex)
1723 {
1724 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1725 if(r<tmpR)
1726 {
1727 pos=i;
1728 tmpR=r;
1729 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1730 }
1731 }
b540d03f 1732 }
bd8c7aef 1733 return pos;
1734}
1735
b5078f5d 1736//__________________________________________________________
1737Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1738{
1739 // check if the cluster survives some quality cut
1740 //
1741 //
1742 Bool_t isGood=kTRUE;
08ff636b 1743 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
fa4287a2 1744 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1745 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1746 if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
b5078f5d 1747
1748 return isGood;
1749}
1750
b540d03f 1751//__________________________________________________________
bd8c7aef 1752Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1753{
1754 // Given a esd track, return whether the track survive all the cuts
1755
1756 // The different quality parameter are first
1757 // retrieved from the track. then it is found out what cuts the
1758 // track did not survive and finally the cuts are imposed.
1759
1760 UInt_t status = esdTrack->GetStatus();
1761
1762 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1763 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1764
1765 Float_t chi2PerClusterITS = -1;
1766 Float_t chi2PerClusterTPC = -1;
1767 if (nClustersITS!=0)
1768 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1769 if (nClustersTPC!=0)
1770 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1771
1772
1773 //DCA cuts
827f9f23 1774 if(fTrackCutsType==kGlobalCut)
1775 {
1776 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1777 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1778 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1779 }
82d09e74 1780
1781
bd8c7aef 1782 Float_t b[2];
1783 Float_t bCov[3];
1784 esdTrack->GetImpactParameters(b,bCov);
1785 if (bCov[0]<=0 || bCov[2]<=0) {
1786 AliDebug(1, "Estimated b resolution lower or equal zero!");
1787 bCov[0]=0; bCov[2]=0;
1788 }
1789
1790 Float_t dcaToVertexXY = b[0];
1791 Float_t dcaToVertexZ = b[1];
1792 Float_t dcaToVertex = -1;
1793
1794 if (fCutDCAToVertex2D)
1795 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1796 else
1797 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1798
1799 // cut the track?
1800
1801 Bool_t cuts[kNCuts];
1802 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1803
1804 // track quality cuts
1805 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1806 cuts[0]=kTRUE;
1807 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1808 cuts[1]=kTRUE;
1809 if (nClustersTPC<fCutMinNClusterTPC)
1810 cuts[2]=kTRUE;
1811 if (nClustersITS<fCutMinNClusterITS)
1812 cuts[3]=kTRUE;
1813 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1814 cuts[4]=kTRUE;
1815 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1816 cuts[5]=kTRUE;
1817 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1818 cuts[6]=kTRUE;
1819 if (fCutDCAToVertex2D && dcaToVertex > 1)
1820 cuts[7] = kTRUE;
1821 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1822 cuts[8] = kTRUE;
1823 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1824 cuts[9] = kTRUE;
1825
827f9f23 1826 if(fTrackCutsType==kGlobalCut)
1827 {
1828 //Require at least one SPD point + anything else in ITS
1829 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1830 cuts[10] = kTRUE;
1831 }
82d09e74 1832
bd8c7aef 1833 Bool_t cut=kFALSE;
827f9f23 1834 for (Int_t i=0; i<kNCuts; i++)
bd8c7aef 1835 if (cuts[i]) {cut = kTRUE;}
1836
1837 // cut the track
1838 if (cut)
1839 return kFALSE;
1840 else
1841 return kTRUE;
1842}
827f9f23 1843
1844
bd8c7aef 1845//__________________________________________________
1846void AliEMCALRecoUtils::InitTrackCuts()
1847{
1848 //Intilize the track cut criteria
5f7714ad 1849 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
bd8c7aef 1850 //Also you can customize the cuts using the setters
82d09e74 1851
5f7714ad 1852 switch (fTrackCutsType)
1853 {
1854 case kTPCOnlyCut:
1855 {
1937171a 1856 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
5f7714ad 1857 //TPC
1858 SetMinNClustersTPC(70);
1859 SetMaxChi2PerClusterTPC(4);
1860 SetAcceptKinkDaughters(kFALSE);
1861 SetRequireTPCRefit(kFALSE);
1862
1863 //ITS
1864 SetRequireITSRefit(kFALSE);
1865 SetMaxDCAToVertexZ(3.2);
1866 SetMaxDCAToVertexXY(2.4);
1867 SetDCAToVertex2D(kTRUE);
1868
1869 break;
1870 }
1871
1872 case kGlobalCut:
1873 {
1937171a 1874 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
5f7714ad 1875 //TPC
1876 SetMinNClustersTPC(70);
1877 SetMaxChi2PerClusterTPC(4);
1878 SetAcceptKinkDaughters(kFALSE);
1879 SetRequireTPCRefit(kTRUE);
1880
1881 //ITS
1882 SetRequireITSRefit(kTRUE);
1883 SetMaxDCAToVertexZ(2);
1884 SetMaxDCAToVertexXY();
1885 SetDCAToVertex2D(kFALSE);
1886
1887 break;
1888 }
0e7de35b 1889
1890 case kLooseCut:
1891 {
1892 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
1893 SetMinNClustersTPC(50);
aef4d246 1894 SetAcceptKinkDaughters(kTRUE);
8e1398ff 1895
1896 break;
0e7de35b 1897 }
5f7714ad 1898 }
bd8c7aef 1899}
83bfd77a 1900
b540d03f 1901//___________________________________________________
d9b3567c 1902void AliEMCALRecoUtils::Print(const Option_t *) const
1903{
1904 // Print Parameters
1905
1906 printf("AliEMCALRecoUtils Settings: \n");
1907 printf("Misalignment shifts\n");
2a71e873 1908 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,
1909 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1910 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1911 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1912 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1913
1914 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1915
fa4287a2 1916 printf("Matching criteria: ");
1917 if(fCutEtaPhiSum)
1918 {
1919 printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1920 }
1921 else if(fCutEtaPhiSeparate)
1922 {
1923 printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1924 }
1925 else
1926 {
1927 printf("Error\n");
1928 printf("please specify your cut criteria\n");
1929 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1930 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1931 }
1932
5f7714ad 1933 printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
bd8c7aef 1934
1935 printf("Track cuts: \n");
fa4287a2 1936 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
bb6f5f0b 1937 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 1938 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1939 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1940 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1941 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1942 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1943
d9b3567c 1944}
96957075 1945
b540d03f 1946//_____________________________________________________________________
3bfc4732 1947void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber){
96957075 1948 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1949 //Do it only once and only if it is requested
1950
3bfc4732 1951 if(!fUseRunCorrectionFactors) return;
1952 if(fRunCorrectionFactorsSet) return;
96957075 1953
3bfc4732 1954 AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber));
96957075 1955
1956 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1957 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1958
1959 SwitchOnRecalibration();
1960 for(Int_t ism = 0; ism < 4; ism++){
1961 for(Int_t icol = 0; icol < 48; icol++){
1962 for(Int_t irow = 0; irow < 24; irow++){
1963 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1964 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1965 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1966 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1967 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1968 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1969 }
1970 }
1971 }
3bfc4732 1972 fRunCorrectionFactorsSet = kTRUE;
96957075 1973}
1974