]> 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
Driver macro for calculating Et em corrections
[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
61160f1f 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);
61160f1f 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;
61160f1f 1121
bb6f5f0b 1122 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1c7a2bf4 1123 if(esdevent)
61160f1f 1124 {
1125 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1126 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1127 if(esdTrack->Pt()<fCutMinTrackPt) continue;
1128 const AliESDfriendTrack* friendTrack = esdTrack->GetFriendTrack();
1129 if(friendTrack && friendTrack->GetTPCOut())
bb6f5f0b 1130 {
61160f1f 1131 //Use TPC Out as starting point if it is available
1132 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
bb6f5f0b 1133 }
61160f1f 1134 else
1135 {
1136 //Otherwise use TPC inner
1137 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1138 }
1139 }
bb6f5f0b 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)
61160f1f 1144 {
1145 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1146 if(!aodTrack) continue;
1147 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1148 if(aodTrack->Pt()<fCutMinTrackPt) continue;
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
61160f1f 1158 {
1159 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1160 return;
1161 }
1162
bb6f5f0b 1163 if(!trackParam) continue;
61160f1f 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++)
61160f1f 1169 {
1170 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1171 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1172 AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation
1173 Float_t tmpEta=-999, tmpPhi=-999;
1174 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
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");
1200 if(aodevent && trackParam) delete trackParam;
1201 return;
1202 }
1203 }//cluster loop
1204 }
1205 else { // external cluster array, not from ESD event
b540d03f 1206 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
61160f1f 1207 {
1208 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1209 if(!cluster){
1210 AliInfo("Cluster not found!!!");
1211 continue;
1212 }
1213 if(!cluster->IsEMCAL()) continue;
1214 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1215 Float_t tmpEta=-999, tmpPhi=-999;
1216 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1217 if(fCutEtaPhiSum)
1218 {
1219 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1220 if(tmpR<dRMax)
1221 {
1222 dRMax=tmpR;
1223 dEtaMax=tmpEta;
1224 dPhiMax=tmpPhi;
1225 index=icl;
1226 }
1227 }
1228 else if(fCutEtaPhiSeparate)
1229 {
1230 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1231 {
1232 dEtaMax = tmpEta;
1233 dPhiMax = tmpPhi;
1234 index=icl;
1235 }
1236 }
1237 else
1238 {
1239 printf("Error: please specify your cut criteria\n");
1240 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1241 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1242 if(aodevent && trackParam) delete trackParam;
1243 return;
1244 }
b540d03f 1245 }//cluster loop
1246 }// external list of clusters
61160f1f 1247
bd8c7aef 1248 if(index>-1)
1249 {
b540d03f 1250 fMatchedTrackIndex ->AddAt(itr,matched);
bd8c7aef 1251 fMatchedClusterIndex->AddAt(index,matched);
fa4287a2 1252 fResidualEta ->AddAt(dEtaMax,matched);
1253 fResidualPhi ->AddAt(dPhiMax,matched);
bd8c7aef 1254 matched++;
1255 }
456126ad 1256 if(aodevent && trackParam) delete trackParam;
bd8c7aef 1257 }//track loop
b540d03f 1258
1259 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1260
1261 fMatchedTrackIndex ->Set(matched);
bd8c7aef 1262 fMatchedClusterIndex->Set(matched);
fa4287a2 1263 fResidualPhi ->Set(matched);
1264 fResidualEta ->Set(matched);
bd8c7aef 1265}
1266
b540d03f 1267//________________________________________________________________________________
b5078f5d 1268Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
bb6f5f0b 1269{
1270 //
1271 // This function returns the index of matched cluster to input track
fa4287a2 1272 // Returns -1 if no match is found
61160f1f 1273
1274
fa4287a2 1275 Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
bb6f5f0b 1276 Int_t index = -1;
61160f1f 1277
bb6f5f0b 1278 AliExternalTrackParam *trackParam=0;
1279 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
1280 if(friendTrack && friendTrack->GetTPCOut())
1281 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1282 else
1283 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
61160f1f 1284
bb6f5f0b 1285 if(!trackParam) return index;
1286 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 1287 {
1288 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1289 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1290 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1291 Float_t tmpEta=-999, tmpPhi=-999;
1292 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1293 if(fCutEtaPhiSum)
bb6f5f0b 1294 {
61160f1f 1295 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1296 if(tmpR<dRMax)
fa4287a2 1297 {
1298 dRMax=tmpR;
1299 dEtaMax=tmpEta;
1300 dPhiMax=tmpPhi;
1301 index=icl;
1302 }
61160f1f 1303 }
1304 else if(fCutEtaPhiSeparate)
1305 {
1306 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
fa4287a2 1307 {
1308 dEtaMax = tmpEta;
1309 dPhiMax = tmpPhi;
1310 index=icl;
1311 }
61160f1f 1312 }
1313 else
1314 {
1315 printf("Error: please specify your cut criteria\n");
1316 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1317 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1318 return -1;
1319 }
1320 }//cluster loop
bb6f5f0b 1321 return index;
1322}
1323
1324//________________________________________________________________________________
fa4287a2 1325Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
bb6f5f0b 1326{
1327 //
1328 //Return the residual by extrapolating a track to a cluster
1329 //
1330 if(!trkParam || !cluster) return kFALSE;
1331 Float_t clsPos[3];
1332 Double_t trkPos[3];
1333 cluster->GetPosition(clsPos); //Has been recalculated
1334 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1335 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1336 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1337 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
bd36717e 1338 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return kFALSE;
bb6f5f0b 1339 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
fa4287a2 1340
1341 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1342 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1343
1344 Float_t clsPhi = (Float_t)clsPosVec.Phi();
1345 if(clsPhi<0) clsPhi+=2*TMath::Pi();
1346 Float_t trkPhi = (Float_t)trkPosVec.Phi();
1347 if(trkPhi<0) trkPhi+=2*TMath::Pi();
1348 tmpPhi = clsPhi-trkPhi; // track cluster matching
1349 tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
1350
bb6f5f0b 1351 return kTRUE;
1352}
1353
1354//________________________________________________________________________________
fa4287a2 1355void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
bd8c7aef 1356{
bb6f5f0b 1357 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 1358 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 1359 //Works with ESDs and AODs
bd8c7aef 1360
bb6f5f0b 1361 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 1362 {
1363 AliDebug(2,"No matched tracks found!\n");
fa4287a2 1364 dEta=999.;
1365 dPhi=999.;
bd8c7aef 1366 return;
1367 }
fa4287a2 1368 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1369 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 1370}
1371//________________________________________________________________________________
fa4287a2 1372void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1373{
1374 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 1375 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 1376 //Works with ESDs and AODs
1377
1378 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1379 {
1380 AliDebug(2,"No matched cluster found!\n");
fa4287a2 1381 dEta=999.;
1382 dPhi=999.;
bb6f5f0b 1383 return;
1384 }
fa4287a2 1385 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1386 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 1387}
1388
1389//__________________________________________________________
1390Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1391{
1392 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1393 //Get the index of matched track to this cluster
1394 //Works with ESDs and AODs
1395
1396 if(IsClusterMatched(clsIndex))
1397 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1398 else
1399 return -1;
bd8c7aef 1400}
1401
b540d03f 1402//__________________________________________________________
bb6f5f0b 1403Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 1404{
bb6f5f0b 1405 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1406 //Get the index of matched cluster to this track
1407 //Works with ESDs and AODs
b540d03f 1408
bb6f5f0b 1409 if(IsTrackMatched(trkIndex))
1410 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 1411 else
1412 return -1;
1413}
1414
bb6f5f0b 1415//__________________________________________________
7cdec71f 1416Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 1417{
1418 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1419 //Returns if the cluster has a match
1420 if(FindMatchedPosForCluster(clsIndex) < 999)
1421 return kTRUE;
1422 else
1423 return kFALSE;
1424}
b540d03f 1425
bd8c7aef 1426//__________________________________________________
7cdec71f 1427Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 1428{
bb6f5f0b 1429 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1430 //Returns if the track has a match
1431 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 1432 return kTRUE;
bd8c7aef 1433 else
1434 return kFALSE;
1435}
bb6f5f0b 1436
b540d03f 1437//__________________________________________________________
bb6f5f0b 1438UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 1439{
bb6f5f0b 1440 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 1441 //Returns the position of the match in the fMatchedClusterIndex array
1442 Float_t tmpR = fCutR;
81efb149 1443 UInt_t pos = 999;
b540d03f 1444
bd8c7aef 1445 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1446 {
fa4287a2 1447 if(fMatchedClusterIndex->At(i)==clsIndex)
1448 {
1449 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1450 if(r<tmpR)
1451 {
1452 pos=i;
1453 tmpR=r;
1454 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1455 }
1456 }
bb6f5f0b 1457 }
1458 return pos;
1459}
1460
1461//__________________________________________________________
1462UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1463{
1464 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1465 //Returns the position of the match in the fMatchedTrackIndex array
1466 Float_t tmpR = fCutR;
1467 UInt_t pos = 999;
1468
1469 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1470 {
fa4287a2 1471 if(fMatchedTrackIndex->At(i)==trkIndex)
1472 {
1473 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1474 if(r<tmpR)
1475 {
1476 pos=i;
1477 tmpR=r;
1478 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1479 }
1480 }
b540d03f 1481 }
bd8c7aef 1482 return pos;
1483}
1484
b5078f5d 1485//__________________________________________________________
1486Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1487{
1488 // check if the cluster survives some quality cut
1489 //
1490 //
1491 Bool_t isGood=kTRUE;
08ff636b 1492 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
fa4287a2 1493 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1494 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1495 if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
b5078f5d 1496
1497 return isGood;
1498}
1499
b540d03f 1500//__________________________________________________________
bd8c7aef 1501Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1502{
1503 // Given a esd track, return whether the track survive all the cuts
1504
1505 // The different quality parameter are first
1506 // retrieved from the track. then it is found out what cuts the
1507 // track did not survive and finally the cuts are imposed.
1508
1509 UInt_t status = esdTrack->GetStatus();
1510
1511 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1512 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1513
1514 Float_t chi2PerClusterITS = -1;
1515 Float_t chi2PerClusterTPC = -1;
1516 if (nClustersITS!=0)
1517 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1518 if (nClustersTPC!=0)
1519 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1520
1521
1522 //DCA cuts
7cdec71f 1523 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
82d09e74 1524 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
7cdec71f 1525 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
82d09e74 1526
1527
bd8c7aef 1528 Float_t b[2];
1529 Float_t bCov[3];
1530 esdTrack->GetImpactParameters(b,bCov);
1531 if (bCov[0]<=0 || bCov[2]<=0) {
1532 AliDebug(1, "Estimated b resolution lower or equal zero!");
1533 bCov[0]=0; bCov[2]=0;
1534 }
1535
1536 Float_t dcaToVertexXY = b[0];
1537 Float_t dcaToVertexZ = b[1];
1538 Float_t dcaToVertex = -1;
1539
1540 if (fCutDCAToVertex2D)
1541 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1542 else
1543 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1544
1545 // cut the track?
1546
1547 Bool_t cuts[kNCuts];
1548 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1549
1550 // track quality cuts
1551 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1552 cuts[0]=kTRUE;
1553 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1554 cuts[1]=kTRUE;
1555 if (nClustersTPC<fCutMinNClusterTPC)
1556 cuts[2]=kTRUE;
1557 if (nClustersITS<fCutMinNClusterITS)
1558 cuts[3]=kTRUE;
1559 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1560 cuts[4]=kTRUE;
1561 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1562 cuts[5]=kTRUE;
1563 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1564 cuts[6]=kTRUE;
1565 if (fCutDCAToVertex2D && dcaToVertex > 1)
1566 cuts[7] = kTRUE;
1567 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1568 cuts[8] = kTRUE;
1569 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1570 cuts[9] = kTRUE;
1571
82d09e74 1572 //Require at least one SPD point + anything else in ITS
1573 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1574 cuts[10] = kTRUE;
1575
bd8c7aef 1576 Bool_t cut=kFALSE;
1577 for (Int_t i=0; i<kNCuts; i++)
1578 if (cuts[i]) {cut = kTRUE;}
1579
1580 // cut the track
1581 if (cut)
1582 return kFALSE;
1583 else
1584 return kTRUE;
1585}
1586//__________________________________________________
1587void AliEMCALRecoUtils::InitTrackCuts()
1588{
1589 //Intilize the track cut criteria
82d09e74 1590 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1591 //Also you can customize the cuts using the setters
82d09e74 1592
1593 //TPC
1594 SetMinNClustersTPC(70);
bd8c7aef 1595 SetMaxChi2PerClusterTPC(4);
1596 SetAcceptKinkDaughters(kFALSE);
82d09e74 1597 SetRequireTPCRefit(kTRUE);
1598
1599 //ITS
1600 SetRequireITSRefit(kTRUE);
1601 SetMaxDCAToVertexZ(2);
1602 SetDCAToVertex2D(kFALSE);
171be15c 1603 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1604 SetMinNClustersITS();
bd8c7aef 1605}
83bfd77a 1606
b540d03f 1607//___________________________________________________
d9b3567c 1608void AliEMCALRecoUtils::Print(const Option_t *) const
1609{
1610 // Print Parameters
1611
1612 printf("AliEMCALRecoUtils Settings: \n");
1613 printf("Misalignment shifts\n");
2a71e873 1614 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,
1615 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1616 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1617 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1618 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1619
1620 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1621
fa4287a2 1622 printf("Matching criteria: ");
1623 if(fCutEtaPhiSum)
1624 {
1625 printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1626 }
1627 else if(fCutEtaPhiSeparate)
1628 {
1629 printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1630 }
1631 else
1632 {
1633 printf("Error\n");
1634 printf("please specify your cut criteria\n");
1635 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1636 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1637 }
1638
bb6f5f0b 1639 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
bd8c7aef 1640
1641 printf("Track cuts: \n");
fa4287a2 1642 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
bb6f5f0b 1643 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 1644 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1645 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1646 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1647 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1648 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1649
d9b3567c 1650}
96957075 1651
b540d03f 1652//_____________________________________________________________________
96957075 1653void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1654 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1655 //Do it only once and only if it is requested
1656
1657 if(!fUseTimeCorrectionFactors) return;
1658 if(fTimeCorrectionFactorsSet) return;
1659
1660 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1661
1662 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1663 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1664
1665 SwitchOnRecalibration();
1666 for(Int_t ism = 0; ism < 4; ism++){
1667 for(Int_t icol = 0; icol < 48; icol++){
1668 for(Int_t irow = 0; irow < 24; irow++){
1669 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1670 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1671 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1672 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1673 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1674 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1675 }
1676 }
1677 }
1678 fTimeCorrectionFactorsSet = kTRUE;
1679}
1680