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