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