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