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