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