Fixing coding rule violations (Sang-Un)
[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();
fa4287a2 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
261//__________________________________________________
094786cc 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
094786cc 388//__________________________________________________
d9b3567c 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}
d9b3567c 467//__________________________________________________
7e0ecb89 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
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
656//________________________________________________________________
fd6df01c 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
679//________________________________________________________________
094786cc 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
83bfd77a 874//____________________________________________________________________________
cb231979 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
950//____________________________________________________________________________
83bfd77a 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 {
bb6f5f0b 1116 AliExternalTrackParam *trackParam=0;
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
1128 trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1129 }
1130 else
1131 {
1132 //Otherwise use TPC inner
1133 trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
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 {
1166 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1167 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1168 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1169 Float_t tmpEta=-999, tmpPhi=-999;
1170 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
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");
1196 return;
1197 }
1198 delete trkPamTmp;
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 {
1203 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1204 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1205 if(!cluster->IsEMCAL()) continue;
1206 Float_t tmpEta=-999, tmpPhi=-999;
1207 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
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");
1233 return;
1234 }
1235 delete trkPamTmp;
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 }
bb6f5f0b 1247 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 {
1279 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1280 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
b5078f5d 1281 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
fa4287a2 1282 Float_t tmpEta=-999, tmpPhi=-999;
1283 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
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 }
1311 delete trkPamTmp;
1312 }//cluster loop
bb6f5f0b 1313 return index;
1314}
1315
1316//________________________________________________________________________________
fa4287a2 1317Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
bb6f5f0b 1318{
1319 //
1320 //Return the residual by extrapolating a track to a cluster
1321 //
1322 if(!trkParam || !cluster) return kFALSE;
1323 Float_t clsPos[3];
1324 Double_t trkPos[3];
1325 cluster->GetPosition(clsPos); //Has been recalculated
1326 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1327 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1328 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1329 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1330 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE;
1331 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
fa4287a2 1332
1333 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1334 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1335
1336 Float_t clsPhi = (Float_t)clsPosVec.Phi();
1337 if(clsPhi<0) clsPhi+=2*TMath::Pi();
1338 Float_t trkPhi = (Float_t)trkPosVec.Phi();
1339 if(trkPhi<0) trkPhi+=2*TMath::Pi();
1340 tmpPhi = clsPhi-trkPhi; // track cluster matching
1341 tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
1342
bb6f5f0b 1343 return kTRUE;
1344}
1345
1346//________________________________________________________________________________
fa4287a2 1347void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
bd8c7aef 1348{
bb6f5f0b 1349 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 1350 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 1351 //Works with ESDs and AODs
bd8c7aef 1352
bb6f5f0b 1353 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 1354 {
1355 AliDebug(2,"No matched tracks found!\n");
fa4287a2 1356 dEta=999.;
1357 dPhi=999.;
bd8c7aef 1358 return;
1359 }
fa4287a2 1360 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1361 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 1362}
1363//________________________________________________________________________________
fa4287a2 1364void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1365{
1366 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 1367 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 1368 //Works with ESDs and AODs
1369
1370 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1371 {
1372 AliDebug(2,"No matched cluster found!\n");
fa4287a2 1373 dEta=999.;
1374 dPhi=999.;
bb6f5f0b 1375 return;
1376 }
fa4287a2 1377 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1378 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 1379}
1380
1381//__________________________________________________________
1382Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1383{
1384 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1385 //Get the index of matched track to this cluster
1386 //Works with ESDs and AODs
1387
1388 if(IsClusterMatched(clsIndex))
1389 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1390 else
1391 return -1;
bd8c7aef 1392}
1393
b540d03f 1394//__________________________________________________________
bb6f5f0b 1395Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 1396{
bb6f5f0b 1397 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1398 //Get the index of matched cluster to this track
1399 //Works with ESDs and AODs
b540d03f 1400
bb6f5f0b 1401 if(IsTrackMatched(trkIndex))
1402 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 1403 else
1404 return -1;
1405}
1406
bb6f5f0b 1407//__________________________________________________
1408Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1409{
1410 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1411 //Returns if the cluster has a match
1412 if(FindMatchedPosForCluster(clsIndex) < 999)
1413 return kTRUE;
1414 else
1415 return kFALSE;
1416}
b540d03f 1417
bd8c7aef 1418//__________________________________________________
bb6f5f0b 1419Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
bd8c7aef 1420{
bb6f5f0b 1421 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1422 //Returns if the track has a match
1423 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 1424 return kTRUE;
bd8c7aef 1425 else
1426 return kFALSE;
1427}
bb6f5f0b 1428
b540d03f 1429//__________________________________________________________
bb6f5f0b 1430UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 1431{
bb6f5f0b 1432 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 1433 //Returns the position of the match in the fMatchedClusterIndex array
1434 Float_t tmpR = fCutR;
81efb149 1435 UInt_t pos = 999;
b540d03f 1436
bd8c7aef 1437 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1438 {
fa4287a2 1439 if(fMatchedClusterIndex->At(i)==clsIndex)
1440 {
1441 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1442 if(r<tmpR)
1443 {
1444 pos=i;
1445 tmpR=r;
1446 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1447 }
1448 }
bb6f5f0b 1449 }
1450 return pos;
1451}
1452
1453//__________________________________________________________
1454UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1455{
1456 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1457 //Returns the position of the match in the fMatchedTrackIndex array
1458 Float_t tmpR = fCutR;
1459 UInt_t pos = 999;
1460
1461 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1462 {
fa4287a2 1463 if(fMatchedTrackIndex->At(i)==trkIndex)
1464 {
1465 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1466 if(r<tmpR)
1467 {
1468 pos=i;
1469 tmpR=r;
1470 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1471 }
1472 }
b540d03f 1473 }
bd8c7aef 1474 return pos;
1475}
1476
b540d03f 1477//__________________________________________________________
b5078f5d 1478Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1479{
1480 // check if the cluster survives some quality cut
1481 //
1482 //
1483 Bool_t isGood=kTRUE;
08ff636b 1484 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
fa4287a2 1485 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1486 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1487 if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
b5078f5d 1488
1489 return isGood;
1490}
1491
1492//__________________________________________________________
bd8c7aef 1493Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1494{
1495 // Given a esd track, return whether the track survive all the cuts
1496
1497 // The different quality parameter are first
1498 // retrieved from the track. then it is found out what cuts the
1499 // track did not survive and finally the cuts are imposed.
1500
1501 UInt_t status = esdTrack->GetStatus();
1502
1503 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1504 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1505
1506 Float_t chi2PerClusterITS = -1;
1507 Float_t chi2PerClusterTPC = -1;
1508 if (nClustersITS!=0)
1509 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1510 if (nClustersTPC!=0)
1511 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1512
1513
1514 //DCA cuts
1515 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1516 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1517 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1518
1519
bd8c7aef 1520 Float_t b[2];
1521 Float_t bCov[3];
1522 esdTrack->GetImpactParameters(b,bCov);
1523 if (bCov[0]<=0 || bCov[2]<=0) {
1524 AliDebug(1, "Estimated b resolution lower or equal zero!");
1525 bCov[0]=0; bCov[2]=0;
1526 }
1527
1528 Float_t dcaToVertexXY = b[0];
1529 Float_t dcaToVertexZ = b[1];
1530 Float_t dcaToVertex = -1;
1531
1532 if (fCutDCAToVertex2D)
1533 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1534 else
1535 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1536
1537 // cut the track?
1538
1539 Bool_t cuts[kNCuts];
1540 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1541
1542 // track quality cuts
1543 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1544 cuts[0]=kTRUE;
1545 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1546 cuts[1]=kTRUE;
1547 if (nClustersTPC<fCutMinNClusterTPC)
1548 cuts[2]=kTRUE;
1549 if (nClustersITS<fCutMinNClusterITS)
1550 cuts[3]=kTRUE;
1551 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1552 cuts[4]=kTRUE;
1553 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1554 cuts[5]=kTRUE;
1555 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1556 cuts[6]=kTRUE;
1557 if (fCutDCAToVertex2D && dcaToVertex > 1)
1558 cuts[7] = kTRUE;
1559 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1560 cuts[8] = kTRUE;
1561 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1562 cuts[9] = kTRUE;
1563
82d09e74 1564 //Require at least one SPD point + anything else in ITS
1565 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1566 cuts[10] = kTRUE;
1567
bd8c7aef 1568 Bool_t cut=kFALSE;
1569 for (Int_t i=0; i<kNCuts; i++)
1570 if (cuts[i]) {cut = kTRUE;}
1571
1572 // cut the track
1573 if (cut)
1574 return kFALSE;
1575 else
1576 return kTRUE;
1577}
1578//__________________________________________________
1579void AliEMCALRecoUtils::InitTrackCuts()
1580{
1581 //Intilize the track cut criteria
82d09e74 1582 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1583 //Also you can customize the cuts using the setters
82d09e74 1584
1585 //TPC
1586 SetMinNClustersTPC(70);
bd8c7aef 1587 SetMaxChi2PerClusterTPC(4);
1588 SetAcceptKinkDaughters(kFALSE);
82d09e74 1589 SetRequireTPCRefit(kTRUE);
1590
1591 //ITS
1592 SetRequireITSRefit(kTRUE);
1593 SetMaxDCAToVertexZ(2);
1594 SetDCAToVertex2D(kFALSE);
171be15c 1595 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1596 SetMinNClustersITS();
bd8c7aef 1597}
83bfd77a 1598
b540d03f 1599//___________________________________________________
d9b3567c 1600void AliEMCALRecoUtils::Print(const Option_t *) const
1601{
1602 // Print Parameters
1603
1604 printf("AliEMCALRecoUtils Settings: \n");
1605 printf("Misalignment shifts\n");
2a71e873 1606 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,
1607 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1608 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1609 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1610 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1611
1612 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1613
fa4287a2 1614 printf("Matching criteria: ");
1615 if(fCutEtaPhiSum)
1616 {
1617 printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1618 }
1619 else if(fCutEtaPhiSeparate)
1620 {
1621 printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1622 }
1623 else
1624 {
1625 printf("Error\n");
1626 printf("please specify your cut criteria\n");
1627 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1628 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1629 }
1630
bb6f5f0b 1631 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
bd8c7aef 1632
1633 printf("Track cuts: \n");
fa4287a2 1634 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
bb6f5f0b 1635 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 1636 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1637 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1638 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1639 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1640 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1641
d9b3567c 1642}
96957075 1643
b540d03f 1644//_____________________________________________________________________
96957075 1645void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1646 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1647 //Do it only once and only if it is requested
1648
1649 if(!fUseTimeCorrectionFactors) return;
1650 if(fTimeCorrectionFactorsSet) return;
1651
1652 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1653
1654 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1655 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1656
1657 SwitchOnRecalibration();
1658 for(Int_t ism = 0; ism < 4; ism++){
1659 for(Int_t icol = 0; icol < 48; icol++){
1660 for(Int_t irow = 0; irow < 24; irow++){
1661 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1662 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1663 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1664 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1665 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1666 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1667 }
1668 }
1669 }
1670 fTimeCorrectionFactorsSet = kTRUE;
1671}
1672