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