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