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