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