]> 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
check if calo ped is zero
[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///////////////////////////////////////////////////////////////////////////////
28
29// --- standard c ---
30
31// standard C++ includes
32//#include <Riostream.h>
33
34// ROOT includes
094786cc 35#include <TGeoManager.h>
36#include <TGeoMatrix.h>
37#include <TGeoBBox.h>
d9b3567c 38
39// STEER includes
d9b3567c 40#include "AliVCluster.h"
41#include "AliVCaloCells.h"
bd8c7aef 42#include "AliVEvent.h"
d9b3567c 43#include "AliLog.h"
83bfd77a 44#include "AliPID.h"
bd8c7aef 45#include "AliESDEvent.h"
bb6f5f0b 46#include "AliAODEvent.h"
bd8c7aef 47#include "AliESDtrack.h"
bb6f5f0b 48#include "AliAODTrack.h"
49#include "AliExternalTrackParam.h"
50#include "AliESDfriendTrack.h"
51#include "AliTrackerBase.h"
b540d03f 52
53// EMCAL includes
54#include "AliEMCALRecoUtils.h"
55#include "AliEMCALGeometry.h"
bd8c7aef 56#include "AliEMCALTrack.h"
96957075 57#include "AliEMCALCalibTimeDepCorrection.h"
b540d03f 58#include "AliEMCALPIDUtils.h"
d9b3567c 59
60ClassImp(AliEMCALRecoUtils)
61
62//______________________________________________
63AliEMCALRecoUtils::AliEMCALRecoUtils():
fd6df01c 64 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
65 fPosAlgo(kUnchanged),fW0(4.),
66 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
78467229 67 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
bd8c7aef 68 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
bb6f5f0b 69 fAODFilterMask(32),
b540d03f 70 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
bb6f5f0b 71 fResidualZ(0x0), fResidualR(0x0), fCutR(10), fCutZ(10), fMass(0.139), fStep(1),
bd8c7aef 72 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
73 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
96957075 74 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
75 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE)
d9b3567c 76{
77//
78 // Constructor.
79 // Initialize all constant values which have to be used
80 // during Reco algorithm execution
81 //
82
b540d03f 83 //Misalignment matrices
fd6df01c 84 for(Int_t i = 0; i < 15 ; i++) {
85 fMisalTransShift[i] = 0.;
b540d03f 86 fMisalRotShift[i] = 0.;
fd6df01c 87 }
b540d03f 88
89 //Non linearity
d9b3567c 90 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
fd6df01c 91 //For kPi0GammaGamma case, but default is no correction
d9b3567c 92 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
93 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
94 fNonLinearityParams[2] = 1.046;
bd8c7aef 95
b540d03f 96 //Track matching
97 fMatchedTrackIndex = new TArrayI();
bd8c7aef 98 fMatchedClusterIndex = new TArrayI();
b540d03f 99 fResidualZ = new TArrayF();
100 fResidualR = new TArrayF();
d9b3567c 101
bd8c7aef 102 InitTrackCuts();
b540d03f 103
104 fPIDUtils = new AliEMCALPIDUtils();
105
bd8c7aef 106
d9b3567c 107}
108
109//______________________________________________________________________
110AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
094786cc 111: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
fd6df01c 112 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
113 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
78467229 114 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
115 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
83bfd77a 116 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
bb6f5f0b 117 fAODFilterMask(reco.fAODFilterMask),
b540d03f 118 fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
bd8c7aef 119 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
120 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
121 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
bb6f5f0b 122 fCutR(reco.fCutR),fCutZ(reco.fCutZ),fMass(reco.fMass), fStep(reco.fStep),
bd8c7aef 123 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
124 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
125 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
126 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
127 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
96957075 128 fPIDUtils(reco.fPIDUtils),
129 fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
d9b3567c 130{
131 //Copy ctor
132
fd6df01c 133 for(Int_t i = 0; i < 15 ; i++) {
134 fMisalRotShift[i] = reco.fMisalRotShift[i];
135 fMisalTransShift[i] = reco.fMisalTransShift[i];
136 }
d9b3567c 137 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
bd8c7aef 138
d9b3567c 139}
140
141
142//______________________________________________________________________
143AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
144{
145 //Assignment operator
146
147 if(this == &reco)return *this;
148 ((TNamed *)this)->operator=(reco);
149
96957075 150 fNonLinearityFunction = reco.fNonLinearityFunction;
151 fParticleType = reco.fParticleType;
152 fPosAlgo = reco.fPosAlgo;
153 fW0 = reco.fW0;
154 fRecalibration = reco.fRecalibration;
094786cc 155 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
96957075 156 fRemoveBadChannels = reco.fRemoveBadChannels;
157 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
158 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
159 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
160 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
bd8c7aef 161
83bfd77a 162
2a71e873 163 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
d9b3567c 164 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
bb6f5f0b 165
166 fAODFilterMask = reco.fAODFilterMask;
d9b3567c 167
96957075 168 fCutR = reco.fCutR;
169 fCutZ = reco.fCutZ;
bb6f5f0b 170 fMass = reco.fMass;
171 fStep = reco.fStep;
bd8c7aef 172
96957075 173 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
174 fCutMinNClusterITS = reco.fCutMinNClusterITS;
175 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
176 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
177 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
178 fCutRequireITSRefit = reco.fCutRequireITSRefit;
179 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
180 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
181 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
182 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
bd8c7aef 183
96957075 184 fPIDUtils = reco.fPIDUtils;
bd8c7aef 185
96957075 186 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
187 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
188
bd8c7aef 189
190 if(reco.fResidualR){
191 // assign or copy construct
192 if(fResidualR){
193 *fResidualR = *reco.fResidualR;
194 }
195 else fResidualR = new TArrayF(*reco.fResidualR);
196 }
197 else{
198 if(fResidualR)delete fResidualR;
199 fResidualR = 0;
200 }
201
202 if(reco.fResidualZ){
203 // assign or copy construct
204 if(fResidualZ){
205 *fResidualZ = *reco.fResidualZ;
206 }
207 else fResidualZ = new TArrayF(*reco.fResidualZ);
208 }
209 else{
210 if(fResidualZ)delete fResidualZ;
211 fResidualZ = 0;
212 }
213
b540d03f 214 if(reco.fMatchedTrackIndex){
215 // assign or copy construct
216 if(fMatchedTrackIndex){
217 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
218 }
219 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
220 }
221 else{
222 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
223 fMatchedTrackIndex = 0;
224 }
bd8c7aef 225
226 if(reco.fMatchedClusterIndex){
227 // assign or copy construct
228 if(fMatchedClusterIndex){
229 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
230 }
231 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
232 }
233 else{
234 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
235 fMatchedClusterIndex = 0;
236 }
237
238
d9b3567c 239 return *this;
240}
241
242
094786cc 243//__________________________________________________
244AliEMCALRecoUtils::~AliEMCALRecoUtils()
245{
246 //Destructor.
247
248 if(fEMCALRecalibrationFactors) {
249 fEMCALRecalibrationFactors->Clear();
250 delete fEMCALRecalibrationFactors;
251 }
fd6df01c 252
253 if(fEMCALBadChannelMap) {
254 fEMCALBadChannelMap->Clear();
255 delete fEMCALBadChannelMap;
256 }
bd8c7aef 257
b540d03f 258 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
bd8c7aef 259 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
260 if(fResidualR) {delete fResidualR; fResidualR=0;}
261 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
262
094786cc 263}
264
fd6df01c 265//_______________________________________________________________
266Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
267{
268 // Given the list of AbsId of the cluster, get the maximum cell and
269 // check if there are fNCellsFromBorder from the calorimeter border
270
271 //If the distance to the border is 0 or negative just exit accept all clusters
272 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
273
cb231979 274 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
275 Bool_t shared = kFALSE;
276 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
fd6df01c 277
83bfd77a 278 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
279 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
fd6df01c 280
281 if(absIdMax==-1) return kFALSE;
282
283 //Check if the cell is close to the borders:
284 Bool_t okrow = kFALSE;
285 Bool_t okcol = kFALSE;
286
287 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
288 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
289 iSM,ieta,iphi));
290 }
291
292 //Check rows/phi
293 if(iSM < 10){
294 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
295 }
296 else{
297 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
298 }
299
300 //Check columns/eta
301 if(!fNoEMCALBorderAtEta0){
302 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
303 }
304 else{
305 if(iSM%2==0){
306 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
307 }
308 else {
309 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
310 }
311 }//eta 0 not checked
312
83bfd77a 313 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
fd6df01c 314 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
315
83bfd77a 316 if (okcol && okrow) {
317 //printf("Accept\n");
318 return kTRUE;
319 }
320 else {
321 //printf("Reject\n");
322 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
323 return kFALSE;
324 }
fd6df01c 325
326}
327
328
329//_________________________________________________________________________________________________________
330Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
331 // Check that in the cluster cells, there is no bad channel of those stored
332 // in fEMCALBadChannelMap or fPHOSBadChannelMap
333
334 if(!fRemoveBadChannels) return kFALSE;
335 if(!fEMCALBadChannelMap) return kFALSE;
336
337 Int_t icol = -1;
338 Int_t irow = -1;
339 Int_t imod = -1;
340 for(Int_t iCell = 0; iCell<nCells; iCell++){
341
342 //Get the column and row
343 Int_t iTower = -1, iIphi = -1, iIeta = -1;
344 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
345 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
346 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
83bfd77a 347 if(GetEMCALChannelStatus(imod, icol, irow)){
348 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
349 return kTRUE;
350 }
fd6df01c 351
352 }// cell cluster loop
353
354 return kFALSE;
355
356}
094786cc 357
d9b3567c 358//__________________________________________________
359Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
360// Correct cluster energy from non linearity functions
361 Float_t energy = cluster->E();
362
363 switch (fNonLinearityFunction) {
364
365 case kPi0MC:
871aee7a 366 {
d9b3567c 367 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
871aee7a 368 //Double_t fNonLinearityParams[0] = 1.001;
369 //Double_t fNonLinearityParams[1] = -0.01264;
370 //Double_t fNonLinearityParams[2] = -0.03632;
371 //Double_t fNonLinearityParams[3] = 0.1798;
372 //Double_t fNonLinearityParams[4] = -0.522;
8cdd1f1f 373 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
d9b3567c 374 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
375 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
376 break;
871aee7a 377 }
d9b3567c 378
379 case kPi0GammaGamma:
871aee7a 380 {
d9b3567c 381 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
96957075 382 //Double_t fNonLinearityParams[0] = 1.04;
383 //Double_t fNonLinearityParams[1] = -0.1445;
871aee7a 384 //Double_t fNonLinearityParams[2] = 1.046;
d9b3567c 385 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
386 break;
871aee7a 387 }
d9b3567c 388
389 case kPi0GammaConversion:
871aee7a 390 {
d9b3567c 391 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
871aee7a 392 //fNonLinearityParams[0] = 0.139393/0.1349766;
393 //fNonLinearityParams[1] = 0.0566186;
394 //fNonLinearityParams[2] = 0.982133;
d9b3567c 395 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
396
397 break;
871aee7a 398 }
399
400 case kBeamTest:
401 {
402 //From beam test, Alexei's results, for different ZS thresholds
403 // th=30 MeV; th = 45 MeV; th = 75 MeV
96957075 404 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
871aee7a 405 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
406 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
96957075 407 //Rescale the param[0] with 1.03
871aee7a 408 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
409
410 break;
411 }
4b58ac4f 412 case kBeamTestCorrected:
413 {
414 //From beam test, corrected for material between beam and EMCAL
415 //fNonLinearityParams[0] = 0.983
416 //fNonLinearityParams[1] = 9.83529e-01;
417 //fNonLinearityParams[2] = -1.84235e+02;
418 //fNonLinearityParams[3] = -2.05019e+00;
419 //fNonLinearityParams[4] = -5.89423e+00;
420 energy *= fNonLinearityParams[0]*( ( fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]) ) +
421 ( (fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())) *
422 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3])) ) );
423
424 break;
425 }
d9b3567c 426
427 case kNoCorrection:
428 AliDebug(2,"No correction on the energy\n");
429 break;
430
431 }
432
433 return energy;
434
435}
436
437//__________________________________________________
094786cc 438Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
439{
440 //Calculate shower depth for a given cluster energy and particle type
441
442 // parameters
cb231979 443 Float_t x0 = 1.31;
094786cc 444 Float_t ecr = 8;
445 Float_t depth = 0;
446
447 switch ( iParticle )
448 {
449 case kPhoton:
fd6df01c 450 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 451 break;
452
453 case kElectron:
fd6df01c 454 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 455 break;
456
457 case kHadron:
458 // hadron
459 // boxes anc. here
460 if(gGeoManager){
461 gGeoManager->cd("ALIC_1/XEN1_1");
462 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
463 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
fd6df01c 464 if(geoSM){
465 TGeoVolume *geoSMVol = geoSM->GetVolume();
466 TGeoShape *geoSMShape = geoSMVol->GetShape();
467 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
468 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
469 else AliFatal("Null GEANT box");
470 }else AliFatal("NULL GEANT node matrix");
094786cc 471 }
472 else{//electron
fd6df01c 473 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 474 }
475
476 break;
477
478 default://photon
fd6df01c 479 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 480 }
481
482 return depth;
483
484}
485
486//__________________________________________________
cb231979 487void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
488 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
d9b3567c 489{
490 //For a given CaloCluster gets the absId of the cell
491 //with maximum energy deposit.
492
493 Double_t eMax = -1.;
494 Double_t eCell = -1.;
094786cc 495 Float_t fraction = 1.;
496 Float_t recalFactor = 1.;
d9b3567c 497 Int_t cellAbsId = -1 ;
094786cc 498
d9b3567c 499 Int_t iTower = -1;
500 Int_t iIphi = -1;
501 Int_t iIeta = -1;
cb231979 502 Int_t iSupMod0= -1;
83bfd77a 503 //printf("---Max?\n");
d9b3567c 504 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 505 cellAbsId = clu->GetCellAbsId(iDig);
506 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 507 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 508 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 509 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
510 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
511 if(iDig==0) iSupMod0=iSupMod;
512 else if(iSupMod0!=iSupMod) {
513 shared = kTRUE;
514 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
515 }
094786cc 516 if(IsRecalibrationOn()) {
094786cc 517 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
518 }
519 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 520 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
094786cc 521 if(eCell > eMax) {
d9b3567c 522 eMax = eCell;
523 absId = cellAbsId;
524 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
525 }
526 }// cell loop
527
528 //Get from the absid the supermodule, tower and eta/phi numbers
529 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
530 //Gives SuperModule and Tower numbers
531 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 532 iIphi, iIeta,iphi,ieta);
533 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
534 //printf("Max end---\n");
d9b3567c 535
536}
537
094786cc 538//________________________________________________________________
539void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
540 //Init EMCAL recalibration factors
541 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
542 //In order to avoid rewriting the same histograms
543 Bool_t oldStatus = TH1::AddDirectoryStatus();
544 TH1::AddDirectory(kFALSE);
545
cb231979 546 fEMCALRecalibrationFactors = new TObjArray(10);
094786cc 547 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
548 //Init the histograms with 1
549 for (Int_t sm = 0; sm < 12; sm++) {
550 for (Int_t i = 0; i < 48; i++) {
551 for (Int_t j = 0; j < 24; j++) {
552 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
553 }
554 }
555 }
556 fEMCALRecalibrationFactors->SetOwner(kTRUE);
557 fEMCALRecalibrationFactors->Compress();
558
559 //In order to avoid rewriting the same histograms
560 TH1::AddDirectory(oldStatus);
561}
562
563
fd6df01c 564//________________________________________________________________
565void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
566 //Init EMCAL bad channels map
567 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
568 //In order to avoid rewriting the same histograms
569 Bool_t oldStatus = TH1::AddDirectoryStatus();
570 TH1::AddDirectory(kFALSE);
571
cb231979 572 fEMCALBadChannelMap = new TObjArray(10);
fd6df01c 573 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
6fe0e6d0 574 for (int i = 0; i < 10; i++) {
fd6df01c 575 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
576 }
577
578 //delete hTemp;
579
580 fEMCALBadChannelMap->SetOwner(kTRUE);
581 fEMCALBadChannelMap->Compress();
582
583 //In order to avoid rewriting the same histograms
584 TH1::AddDirectory(oldStatus);
585}
586
094786cc 587//________________________________________________________________
588void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
589 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
590
591 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
592 UShort_t * index = cluster->GetCellsAbsId() ;
593 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
594 Int_t ncells = cluster->GetNCells();
595
596 //Initialize some used variables
597 Float_t energy = 0;
598 Int_t absId = -1;
599 Int_t icol = -1, irow = -1, imod=1;
600 Float_t factor = 1, frac = 0;
601
602 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
603 for(Int_t icell = 0; icell < ncells; icell++){
604 absId = index[icell];
605 frac = fraction[icell];
606 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
607 Int_t iTower = -1, iIphi = -1, iIeta = -1;
608 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
609 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
610 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
611 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
612 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
613 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
614
615 energy += cells->GetCellAmplitude(absId)*factor*frac;
616 }
617
618
619 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
620
621 cluster->SetE(energy);
622
623}
624
625
d9b3567c 626//__________________________________________________
094786cc 627void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 628{
629 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
630
094786cc 631 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
632 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 633 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 634
635}
636
637//__________________________________________________
638void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
639{
640 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
641 // The algorithm is a copy of what is done in AliEMCALRecPoint
642
643 Double_t eCell = 0.;
644 Float_t fraction = 1.;
645 Float_t recalFactor = 1.;
646
647 Int_t absId = -1;
648 Int_t iTower = -1, iIphi = -1, iIeta = -1;
649 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
650 Float_t weight = 0., totalWeight=0.;
651 Float_t newPos[3] = {0,0,0};
652 Double_t pLocal[3], pGlobal[3];
cb231979 653 Bool_t shared = kFALSE;
654
094786cc 655 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
cb231979 656 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 657 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
658
83bfd77a 659 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
660
094786cc 661 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
662 absId = clu->GetCellAbsId(iDig);
663 fraction = clu->GetCellAmplitudeFraction(iDig);
664 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
665 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
666 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
667
668 if(IsRecalibrationOn()) {
669 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
670 }
671 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
672
673 weight = GetCellWeight(eCell,clEnergy);
83bfd77a 674 //printf("cell energy %f, weight %f\n",eCell,weight);
094786cc 675 totalWeight += weight;
676 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 677 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 678 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 679 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
680
094786cc 681 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
682
683 }// cell loop
684
685 if(totalWeight>0){
686 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
687 }
688
094786cc 689 //Float_t pos[]={0,0,0};
690 //clu->GetPosition(pos);
691 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 692 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 693
694 if(iSupModMax > 1) {//sector 1
695 newPos[0] +=fMisalTransShift[3];//-=3.093;
696 newPos[1] +=fMisalTransShift[4];//+=6.82;
697 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 698 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
699
094786cc 700 }
701 else {//sector 0
702 newPos[0] +=fMisalTransShift[0];//+=1.134;
703 newPos[1] +=fMisalTransShift[1];//+=8.2;
704 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 705 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
706
094786cc 707 }
83bfd77a 708 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
709
094786cc 710 clu->SetPosition(newPos);
711
094786cc 712}
713
714//__________________________________________________
715void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
716{
717 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
718 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
719
720 Double_t eCell = 1.;
721 Float_t fraction = 1.;
722 Float_t recalFactor = 1.;
723
724 Int_t absId = -1;
d9b3567c 725 Int_t iTower = -1;
094786cc 726 Int_t iIphi = -1, iIeta = -1;
727 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 728 Int_t iphi = -1, ieta =-1;
cb231979 729 Bool_t shared = kFALSE;
730
d9b3567c 731 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
cb231979 732 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 733 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
734
d9b3567c 735 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 736 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
737 Int_t startingSM = -1;
d9b3567c 738
739 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 740 absId = clu->GetCellAbsId(iDig);
741 fraction = clu->GetCellAmplitudeFraction(iDig);
742 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
743 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
d9b3567c 744 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 745
d9b3567c 746 if (iDig==0) startingSM = iSupMod;
747 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 748
749 eCell = cells->GetCellAmplitude(absId);
d9b3567c 750
094786cc 751 if(IsRecalibrationOn()) {
752 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
753 }
754 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 755
094786cc 756 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 757 if(weight < 0) weight = 0;
758 totalWeight += weight;
759 weightedCol += ieta*weight;
760 weightedRow += iphi*weight;
761
762 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
763
094786cc 764 }// cell loop
765
d9b3567c 766 Float_t xyzNew[]={0.,0.,0.};
767 if(areInSameSM == kTRUE) {
768 //printf("In Same SM\n");
769 weightedCol = weightedCol/totalWeight;
770 weightedRow = weightedRow/totalWeight;
094786cc 771 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 772 }
773 else {
774 //printf("In Different SM\n");
094786cc 775 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 776 }
d9b3567c 777
094786cc 778 clu->SetPosition(xyzNew);
d9b3567c 779
780}
781
cb231979 782//____________________________________________________________________________
783void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
784
785 //re-evaluate distance to bad channel with updated bad map
786
78467229 787 if(!fRecalDistToBadChannels) return;
cb231979 788
789 //Get channels map of the supermodule where the cluster is.
cb231979 790 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
791 Bool_t shared = kFALSE;
792 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
793 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
794
795 Int_t dRrow, dRcol;
796 Float_t minDist = 10000.;
797 Float_t dist = 0.;
798
799 //Loop on tower status map
800 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
801 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
802 //Check if tower is bad.
803 if(hMap->GetBinContent(icol,irow)==0) continue;
804 //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 805 // iSupMod,icol, irow, icolM,irowM);
cb231979 806
807 dRrow=TMath::Abs(irowM-irow);
808 dRcol=TMath::Abs(icolM-icol);
809 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
810 if(dist < minDist){
811 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
812 minDist = dist;
813 }
814
815 }
816 }
817
818 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
819 if (shared) {
820 TH2D* hMap2 = 0;
821 Int_t iSupMod2 = -1;
822
823 //The only possible combinations are (0,1), (2,3) ... (8,9)
824 if(iSupMod%2) iSupMod2 = iSupMod-1;
825 else iSupMod2 = iSupMod+1;
826 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
827
828 //Loop on tower status map of second super module
829 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
830 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
831 //Check if tower is bad.
832 if(hMap2->GetBinContent(icol,irow)==0) continue;
833 //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",
834 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
835
836 dRrow=TMath::Abs(irow-irowM);
837
838 if(iSupMod%2) {
839 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
840 }
841 else {
842 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
843 }
844
845 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
846 if(dist < minDist) minDist = dist;
847
848 }
849 }
850
851 }// shared cluster in 2 SuperModules
78467229 852
6fe0e6d0 853 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
854 cluster->SetDistanceToBadChannel(minDist);
cb231979 855
856}
857
83bfd77a 858//____________________________________________________________________________
859void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
860
861 //re-evaluate identification parameters with bayesian
862
863 if ( cluster->GetM02() != 0)
864 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
865
866 Float_t pidlist[AliPID::kSPECIESN+1];
867 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
868
869 cluster->SetPID(pidlist);
870
871}
872
873//____________________________________________________________________________
874void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
875{
876 // Calculates new center of gravity in the local EMCAL-module coordinates
877 // and tranfers into global ALICE coordinates
878 // Calculates Dispersion and main axis
879
880 Int_t nstat = 0;
881 Float_t wtot = 0. ;
882 Double_t eCell = 0.;
883 Float_t fraction = 1.;
884 Float_t recalFactor = 1.;
885
886 Int_t iSupMod = -1;
887 Int_t iTower = -1;
888 Int_t iIphi = -1;
889 Int_t iIeta = -1;
890 Int_t iphi = -1;
891 Int_t ieta = -1;
892 Double_t etai = -1.;
893 Double_t phii = -1.;
894
895 Double_t w = 0.;
896 Double_t d = 0.;
897 Double_t dxx = 0.;
898 Double_t dzz = 0.;
899 Double_t dxz = 0.;
900 Double_t xmean = 0.;
901 Double_t zmean = 0.;
902
903 //Loop on cells
904 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
905
906 //Get from the absid the supermodule, tower and eta/phi numbers
907 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
908 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
909
910 //Get the cell energy, if recalibration is on, apply factors
911 fraction = cluster->GetCellAmplitudeFraction(iDigit);
912 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
913 if(IsRecalibrationOn()) {
914 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
915 }
916 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
917
918 if(cluster->E() > 0 && eCell > 0){
919
920 w = GetCellWeight(eCell,cluster->E());
921
922 etai=(Double_t)ieta;
923 phii=(Double_t)iphi;
924 if(w > 0.0) {
925 wtot += w ;
926 nstat++;
927 //Shower shape
928 dxx += w * etai * etai ;
929 xmean+= w * etai ;
930 dzz += w * phii * phii ;
931 zmean+= w * phii ;
932 dxz += w * etai * phii ;
933 }
934 }
935 else
936 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
937 }//cell loop
938
939 //Normalize to the weight
940 if (wtot > 0) {
941 xmean /= wtot ;
942 zmean /= wtot ;
943 }
944 else
945 AliError(Form("Wrong weight %f\n", wtot));
946
947 //Calculate dispersion
948 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
949
950 //Get from the absid the supermodule, tower and eta/phi numbers
951 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
952 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
953
954 //Get the cell energy, if recalibration is on, apply factors
955 fraction = cluster->GetCellAmplitudeFraction(iDigit);
956 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
957 if(IsRecalibrationOn()) {
958 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
959 }
960 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
961
962 if(cluster->E() > 0 && eCell > 0){
963
964 w = GetCellWeight(eCell,cluster->E());
965
966 etai=(Double_t)ieta;
967 phii=(Double_t)iphi;
968 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
969 }
970 else
971 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
972 }// cell loop
973
974 //Normalize to the weigth and set shower shape parameters
975 if (wtot > 0 && nstat > 1) {
976 d /= wtot ;
977 dxx /= wtot ;
978 dzz /= wtot ;
979 dxz /= wtot ;
980 dxx -= xmean * xmean ;
981 dzz -= zmean * zmean ;
982 dxz -= xmean * zmean ;
983 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
984 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
985 }
986 else{
987 d=0. ;
988 cluster->SetM20(0.) ;
989 cluster->SetM02(0.) ;
990 }
991
992 if (d>=0)
993 cluster->SetDispersion(TMath::Sqrt(d)) ;
994 else
995 cluster->SetDispersion(0) ;
83bfd77a 996}
997
b540d03f 998//____________________________________________________________________________
bb6f5f0b 999void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr, TString dataType)
bd8c7aef 1000{
bb6f5f0b 1001 //Use dataType to indicate the input event is AOD or ESD
bd8c7aef 1002 //This function should be called before the cluster loop
1003 //Before call this function, please recalculate the cluster positions
1004 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1005 //Store matched cluster indexes and residuals
bb6f5f0b 1006
b540d03f 1007 fMatchedTrackIndex ->Reset();
bd8c7aef 1008 fMatchedClusterIndex->Reset();
b540d03f 1009 fResidualZ ->Reset();
1010 fResidualR ->Reset();
bd8c7aef 1011
67241b7e 1012 fMatchedTrackIndex ->Set(500);
1013 fMatchedClusterIndex->Set(500);
1014 fResidualZ ->Set(500);
1015 fResidualR ->Set(500);
bd8c7aef 1016
1017 Int_t matched=0;
bb6f5f0b 1018 Double_t cv[21];
1019 for (Int_t i=0; i<21;i++) cv[i]=0;
bd8c7aef 1020 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1021 {
bb6f5f0b 1022 AliExternalTrackParam *trackParam=0;
1023
1024 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1025 if(dataType.Contains("ESD"))
1026 {
1027 AliESDtrack *esdTrack = ((AliESDEvent*)event)->GetTrack(itr);
1028 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1029 const AliESDfriendTrack* friendTrack = esdTrack->GetFriendTrack();
1030 if(friendTrack && friendTrack->GetTPCOut())
1031 {
1032 //Use TPC Out as starting point if it is available
1033 trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1034 }
1035 else
1036 {
1037 //Otherwise use TPC inner
1038 trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
1039 }
1040 }
1041
1042 //If the input event is AOD, the starting point for extrapolation is at vertex
1043 //AOD tracks are selected according to its bit.
1044 else if(dataType.Contains("AOD"))
1045 {
1046 AliAODTrack *aodTrack = ((AliAODEvent*)event)->GetTrack(itr);
1047 if(!aodTrack) continue;
1048 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1049 Double_t pos[3],mom[3];
1050 aodTrack->GetXYZ(pos);
1051 aodTrack->GetPxPyPz(mom);
1052 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()));
1053 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1054 }
bd8c7aef 1055
bb6f5f0b 1056 //Return if the input data is not "AOD" or "ESD"
1057 else
1058 {
1059 printf("Wrong input data type %s! Should be \"AOD\" or \"ESD\"\n",dataType.Data());
1060 return;
1061 }
1062
1063 if(!trackParam) continue;
1064
1065 Float_t dRMax = fCutR, dZMax=fCutZ;
bd8c7aef 1066 Int_t index = -1;
b540d03f 1067 if(!clusterArr){// get clusters from event
1068 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1069 {
bb6f5f0b 1070 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
b540d03f 1071 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
bb6f5f0b 1072 if(!cluster->IsEMCAL()) continue;
1073 Float_t tmpR=-1, tmpZ=-1;
1074 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
b540d03f 1075 if(tmpR<dRMax)
1076 {
1077 dRMax=tmpR;
1078 dZMax=tmpZ;
1079 index=icl;
1080 }
bb6f5f0b 1081 delete trkPamTmp;
b540d03f 1082 }//cluster loop
1083 } else { // external cluster array, not from ESD event
1084 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1085 {
bb6f5f0b 1086 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
b540d03f 1087 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1088 if(!cluster->IsEMCAL()) continue;
bb6f5f0b 1089 Float_t tmpR=-1, tmpZ=-1;
1090 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
b540d03f 1091 if(tmpR<dRMax)
1092 {
1093 dRMax=tmpR;
1094 dZMax=tmpZ;
1095 index=icl;
1096 }
bb6f5f0b 1097 delete trkPamTmp;
b540d03f 1098 }//cluster loop
1099 }// external list of clusters
1100
bd8c7aef 1101 if(index>-1)
1102 {
b540d03f 1103 fMatchedTrackIndex ->AddAt(itr,matched);
bd8c7aef 1104 fMatchedClusterIndex->AddAt(index,matched);
b540d03f 1105 fResidualZ ->AddAt(dZMax,matched);
1106 fResidualR ->AddAt(dRMax,matched);
bd8c7aef 1107 matched++;
1108 }
bb6f5f0b 1109 delete trackParam;
bd8c7aef 1110 }//track loop
b540d03f 1111
1112 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1113
1114 fMatchedTrackIndex ->Set(matched);
bd8c7aef 1115 fMatchedClusterIndex->Set(matched);
b540d03f 1116 fResidualZ ->Set(matched);
1117 fResidualR ->Set(matched);
bd8c7aef 1118}
1119
b540d03f 1120//________________________________________________________________________________
bb6f5f0b 1121Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event)
1122{
1123 //
1124 // This function returns the index of matched cluster to input track
1125 // Cut on match is dR<10cm by default. Returns -1 if no match is found
1126
1127
1128 Float_t dRMax = fCutR;
1129 Int_t index = -1;
1130
1131 AliExternalTrackParam *trackParam=0;
1132 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
1133 if(friendTrack && friendTrack->GetTPCOut())
1134 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1135 else
1136 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1137
1138 if(!trackParam) return index;
1139 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1140 {
1141 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1142 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1143 if(!cluster->IsEMCAL()) continue;
1144 Float_t tmpR=-1, tmpZ=-1;
1145 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1146 if(tmpR>-1 && tmpR<dRMax)
1147 {
1148 dRMax=tmpR;
1149 index=icl;
1150 }
1151 delete trkPamTmp;
1152 }//cluster loop
1153 return index;
1154}
1155
1156//________________________________________________________________________________
1157Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpR, Float_t &tmpZ)
1158{
1159 //
1160 //Return the residual by extrapolating a track to a cluster
1161 //
1162 if(!trkParam || !cluster) return kFALSE;
1163 Float_t clsPos[3];
1164 Double_t trkPos[3];
1165 cluster->GetPosition(clsPos); //Has been recalculated
1166 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1167 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1168 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1169 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1170 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE;
1171 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1172 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1173 tmpZ = clsPos[2]-trkPos[2];
1174 return kTRUE;
1175}
1176
1177//________________________________________________________________________________
1178void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dR, Float_t &dZ)
bd8c7aef 1179{
bb6f5f0b 1180 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1181 //Get the residuals dR and dZ for this cluster to the closest track
1182 //Works with ESDs and AODs
bd8c7aef 1183
bb6f5f0b 1184 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 1185 {
1186 AliDebug(2,"No matched tracks found!\n");
1187 dR=999.;
1188 dZ=999.;
1189 return;
1190 }
bb6f5f0b 1191 dR = fResidualR->At(FindMatchedPosForCluster(clsIndex));
1192 dZ = fResidualZ->At(FindMatchedPosForCluster(clsIndex));
1193}
1194//________________________________________________________________________________
1195void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dR, Float_t &dZ)
1196{
1197 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1198 //Get the residuals dR and dZ for this track to the closest cluster
1199 //Works with ESDs and AODs
1200
1201 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1202 {
1203 AliDebug(2,"No matched cluster found!\n");
1204 dR=999.;
1205 dZ=999.;
1206 return;
1207 }
1208 dR = fResidualR->At(FindMatchedPosForTrack(trkIndex));
1209 dZ = fResidualZ->At(FindMatchedPosForTrack(trkIndex));
1210}
1211
1212//__________________________________________________________
1213Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1214{
1215 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1216 //Get the index of matched track to this cluster
1217 //Works with ESDs and AODs
1218
1219 if(IsClusterMatched(clsIndex))
1220 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1221 else
1222 return -1;
bd8c7aef 1223}
1224
b540d03f 1225//__________________________________________________________
bb6f5f0b 1226Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 1227{
bb6f5f0b 1228 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1229 //Get the index of matched cluster to this track
1230 //Works with ESDs and AODs
b540d03f 1231
bb6f5f0b 1232 if(IsTrackMatched(trkIndex))
1233 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 1234 else
1235 return -1;
1236}
1237
bb6f5f0b 1238//__________________________________________________
1239Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1240{
1241 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1242 //Returns if the cluster has a match
1243 if(FindMatchedPosForCluster(clsIndex) < 999)
1244 return kTRUE;
1245 else
1246 return kFALSE;
1247}
b540d03f 1248
bd8c7aef 1249//__________________________________________________
bb6f5f0b 1250Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
bd8c7aef 1251{
bb6f5f0b 1252 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1253 //Returns if the track has a match
1254 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 1255 return kTRUE;
bd8c7aef 1256 else
1257 return kFALSE;
1258}
bb6f5f0b 1259
b540d03f 1260//__________________________________________________________
bb6f5f0b 1261UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 1262{
bb6f5f0b 1263 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 1264 //Returns the position of the match in the fMatchedClusterIndex array
1265 Float_t tmpR = fCutR;
81efb149 1266 UInt_t pos = 999;
b540d03f 1267
bd8c7aef 1268 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1269 {
bb6f5f0b 1270 if(fMatchedClusterIndex->At(i)==clsIndex && fResidualR->At(i)<tmpR)
1271 {
1272 pos=i;
1273 tmpR=fResidualR->At(i);
1274 AliDebug(3,Form("Matched cluster index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1275 }
1276 }
1277 return pos;
1278}
1279
1280//__________________________________________________________
1281UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1282{
1283 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1284 //Returns the position of the match in the fMatchedTrackIndex array
1285 Float_t tmpR = fCutR;
1286 UInt_t pos = 999;
1287
1288 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1289 {
1290 if(fMatchedTrackIndex->At(i)==trkIndex && fResidualR->At(i)<tmpR)
bd8c7aef 1291 {
b540d03f 1292 pos=i;
1293 tmpR=fResidualR->At(i);
bb6f5f0b 1294 AliDebug(3,Form("Matched track index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
bd8c7aef 1295 }
b540d03f 1296 }
bd8c7aef 1297 return pos;
1298}
1299
b540d03f 1300//__________________________________________________________
bd8c7aef 1301Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1302{
1303 // Given a esd track, return whether the track survive all the cuts
1304
1305 // The different quality parameter are first
1306 // retrieved from the track. then it is found out what cuts the
1307 // track did not survive and finally the cuts are imposed.
1308
1309 UInt_t status = esdTrack->GetStatus();
1310
1311 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1312 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1313
1314 Float_t chi2PerClusterITS = -1;
1315 Float_t chi2PerClusterTPC = -1;
1316 if (nClustersITS!=0)
1317 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1318 if (nClustersTPC!=0)
1319 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1320
1321
1322 //DCA cuts
1323 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1324 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1325 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1326
1327
bd8c7aef 1328 Float_t b[2];
1329 Float_t bCov[3];
1330 esdTrack->GetImpactParameters(b,bCov);
1331 if (bCov[0]<=0 || bCov[2]<=0) {
1332 AliDebug(1, "Estimated b resolution lower or equal zero!");
1333 bCov[0]=0; bCov[2]=0;
1334 }
1335
1336 Float_t dcaToVertexXY = b[0];
1337 Float_t dcaToVertexZ = b[1];
1338 Float_t dcaToVertex = -1;
1339
1340 if (fCutDCAToVertex2D)
1341 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1342 else
1343 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1344
1345 // cut the track?
1346
1347 Bool_t cuts[kNCuts];
1348 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1349
1350 // track quality cuts
1351 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1352 cuts[0]=kTRUE;
1353 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1354 cuts[1]=kTRUE;
1355 if (nClustersTPC<fCutMinNClusterTPC)
1356 cuts[2]=kTRUE;
1357 if (nClustersITS<fCutMinNClusterITS)
1358 cuts[3]=kTRUE;
1359 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1360 cuts[4]=kTRUE;
1361 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1362 cuts[5]=kTRUE;
1363 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1364 cuts[6]=kTRUE;
1365 if (fCutDCAToVertex2D && dcaToVertex > 1)
1366 cuts[7] = kTRUE;
1367 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1368 cuts[8] = kTRUE;
1369 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1370 cuts[9] = kTRUE;
1371
82d09e74 1372 //Require at least one SPD point + anything else in ITS
1373 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1374 cuts[10] = kTRUE;
1375
bd8c7aef 1376 Bool_t cut=kFALSE;
1377 for (Int_t i=0; i<kNCuts; i++)
1378 if (cuts[i]) {cut = kTRUE;}
1379
1380 // cut the track
1381 if (cut)
1382 return kFALSE;
1383 else
1384 return kTRUE;
1385}
1386//__________________________________________________
1387void AliEMCALRecoUtils::InitTrackCuts()
1388{
1389 //Intilize the track cut criteria
82d09e74 1390 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1391 //Also you can customize the cuts using the setters
82d09e74 1392
1393 //TPC
1394 SetMinNClustersTPC(70);
bd8c7aef 1395 SetMaxChi2PerClusterTPC(4);
1396 SetAcceptKinkDaughters(kFALSE);
82d09e74 1397 SetRequireTPCRefit(kTRUE);
1398
1399 //ITS
1400 SetRequireITSRefit(kTRUE);
1401 SetMaxDCAToVertexZ(2);
1402 SetDCAToVertex2D(kFALSE);
171be15c 1403 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1404 SetMinNClustersITS();
bd8c7aef 1405}
83bfd77a 1406
b540d03f 1407//___________________________________________________
d9b3567c 1408void AliEMCALRecoUtils::Print(const Option_t *) const
1409{
1410 // Print Parameters
1411
1412 printf("AliEMCALRecoUtils Settings: \n");
1413 printf("Misalignment shifts\n");
2a71e873 1414 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,
1415 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1416 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1417 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1418 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1419
1420 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1421
1422 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
bb6f5f0b 1423 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
bd8c7aef 1424
1425 printf("Track cuts: \n");
bb6f5f0b 1426 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 1427 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1428 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1429 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1430 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1431 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1432
d9b3567c 1433}
96957075 1434
b540d03f 1435//_____________________________________________________________________
96957075 1436void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1437 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1438 //Do it only once and only if it is requested
1439
1440 if(!fUseTimeCorrectionFactors) return;
1441 if(fTimeCorrectionFactorsSet) return;
1442
1443 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1444
1445 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1446 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1447
1448 SwitchOnRecalibration();
1449 for(Int_t ism = 0; ism < 4; ism++){
1450 for(Int_t icol = 0; icol < 48; icol++){
1451 for(Int_t irow = 0; irow < 24; irow++){
1452 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1453 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1454 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1455 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1456 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1457 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1458 }
1459 }
1460 }
1461 fTimeCorrectionFactorsSet = kTRUE;
1462}
1463