]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALRecoUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Verbose printout commented out
[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///////////////////////////////////////////////////////////////////////////////
28
29// --- standard c ---
30
31// standard C++ includes
32//#include <Riostream.h>
33
34// ROOT includes
094786cc 35#include <TGeoManager.h>
36#include <TGeoMatrix.h>
37#include <TGeoBBox.h>
d9b3567c 38
39// STEER includes
d9b3567c 40#include "AliVCluster.h"
41#include "AliVCaloCells.h"
bd8c7aef 42#include "AliVEvent.h"
d9b3567c 43#include "AliLog.h"
83bfd77a 44#include "AliPID.h"
bd8c7aef 45#include "AliESDEvent.h"
46#include "AliESDtrack.h"
b540d03f 47
48// EMCAL includes
49#include "AliEMCALRecoUtils.h"
50#include "AliEMCALGeometry.h"
bd8c7aef 51#include "AliEMCALTrack.h"
96957075 52#include "AliEMCALCalibTimeDepCorrection.h"
b540d03f 53#include "AliEMCALPIDUtils.h"
d9b3567c 54
55ClassImp(AliEMCALRecoUtils)
56
57//______________________________________________
58AliEMCALRecoUtils::AliEMCALRecoUtils():
fd6df01c 59 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
60 fPosAlgo(kUnchanged),fW0(4.),
61 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
78467229 62 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
bd8c7aef 63 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
b540d03f 64 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
65 fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
bd8c7aef 66 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
67 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
96957075 68 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
69 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE)
d9b3567c 70{
71//
72 // Constructor.
73 // Initialize all constant values which have to be used
74 // during Reco algorithm execution
75 //
76
b540d03f 77 //Misalignment matrices
fd6df01c 78 for(Int_t i = 0; i < 15 ; i++) {
79 fMisalTransShift[i] = 0.;
b540d03f 80 fMisalRotShift[i] = 0.;
fd6df01c 81 }
b540d03f 82
83 //Non linearity
d9b3567c 84 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
fd6df01c 85 //For kPi0GammaGamma case, but default is no correction
d9b3567c 86 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
87 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
88 fNonLinearityParams[2] = 1.046;
bd8c7aef 89
b540d03f 90 //Track matching
91 fMatchedTrackIndex = new TArrayI();
bd8c7aef 92 fMatchedClusterIndex = new TArrayI();
b540d03f 93 fResidualZ = new TArrayF();
94 fResidualR = new TArrayF();
d9b3567c 95
bd8c7aef 96 InitTrackCuts();
b540d03f 97
98 fPIDUtils = new AliEMCALPIDUtils();
99
bd8c7aef 100
d9b3567c 101}
102
103//______________________________________________________________________
104AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
094786cc 105: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
fd6df01c 106 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
107 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
78467229 108 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
109 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
83bfd77a 110 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
b540d03f 111 fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
bd8c7aef 112 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
113 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
114 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
115 fCutR(reco.fCutR),fCutZ(reco.fCutZ),
116 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
117 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
118 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
119 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
120 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
96957075 121 fPIDUtils(reco.fPIDUtils),
122 fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
d9b3567c 123{
124 //Copy ctor
125
fd6df01c 126 for(Int_t i = 0; i < 15 ; i++) {
127 fMisalRotShift[i] = reco.fMisalRotShift[i];
128 fMisalTransShift[i] = reco.fMisalTransShift[i];
129 }
d9b3567c 130 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
bd8c7aef 131
d9b3567c 132}
133
134
135//______________________________________________________________________
136AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
137{
138 //Assignment operator
139
140 if(this == &reco)return *this;
141 ((TNamed *)this)->operator=(reco);
142
96957075 143 fNonLinearityFunction = reco.fNonLinearityFunction;
144 fParticleType = reco.fParticleType;
145 fPosAlgo = reco.fPosAlgo;
146 fW0 = reco.fW0;
147 fRecalibration = reco.fRecalibration;
094786cc 148 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
96957075 149 fRemoveBadChannels = reco.fRemoveBadChannels;
150 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
151 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
152 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
153 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
bd8c7aef 154
83bfd77a 155
2a71e873 156 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
d9b3567c 157 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
158
96957075 159 fCutR = reco.fCutR;
160 fCutZ = reco.fCutZ;
bd8c7aef 161
96957075 162 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
163 fCutMinNClusterITS = reco.fCutMinNClusterITS;
164 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
165 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
166 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
167 fCutRequireITSRefit = reco.fCutRequireITSRefit;
168 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
169 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
170 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
171 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
bd8c7aef 172
96957075 173 fPIDUtils = reco.fPIDUtils;
bd8c7aef 174
96957075 175 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
176 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
177
bd8c7aef 178
179 if(reco.fResidualR){
180 // assign or copy construct
181 if(fResidualR){
182 *fResidualR = *reco.fResidualR;
183 }
184 else fResidualR = new TArrayF(*reco.fResidualR);
185 }
186 else{
187 if(fResidualR)delete fResidualR;
188 fResidualR = 0;
189 }
190
191 if(reco.fResidualZ){
192 // assign or copy construct
193 if(fResidualZ){
194 *fResidualZ = *reco.fResidualZ;
195 }
196 else fResidualZ = new TArrayF(*reco.fResidualZ);
197 }
198 else{
199 if(fResidualZ)delete fResidualZ;
200 fResidualZ = 0;
201 }
202
b540d03f 203 if(reco.fMatchedTrackIndex){
204 // assign or copy construct
205 if(fMatchedTrackIndex){
206 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
207 }
208 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
209 }
210 else{
211 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
212 fMatchedTrackIndex = 0;
213 }
bd8c7aef 214
215 if(reco.fMatchedClusterIndex){
216 // assign or copy construct
217 if(fMatchedClusterIndex){
218 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
219 }
220 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
221 }
222 else{
223 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
224 fMatchedClusterIndex = 0;
225 }
226
227
d9b3567c 228 return *this;
229}
230
231
094786cc 232//__________________________________________________
233AliEMCALRecoUtils::~AliEMCALRecoUtils()
234{
235 //Destructor.
236
237 if(fEMCALRecalibrationFactors) {
238 fEMCALRecalibrationFactors->Clear();
239 delete fEMCALRecalibrationFactors;
240 }
fd6df01c 241
242 if(fEMCALBadChannelMap) {
243 fEMCALBadChannelMap->Clear();
244 delete fEMCALBadChannelMap;
245 }
bd8c7aef 246
b540d03f 247 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
bd8c7aef 248 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
249 if(fResidualR) {delete fResidualR; fResidualR=0;}
250 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
251
094786cc 252}
253
fd6df01c 254//_______________________________________________________________
255Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
256{
257 // Given the list of AbsId of the cluster, get the maximum cell and
258 // check if there are fNCellsFromBorder from the calorimeter border
259
260 //If the distance to the border is 0 or negative just exit accept all clusters
261 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
262
cb231979 263 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
264 Bool_t shared = kFALSE;
265 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
fd6df01c 266
83bfd77a 267 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
268 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
fd6df01c 269
270 if(absIdMax==-1) return kFALSE;
271
272 //Check if the cell is close to the borders:
273 Bool_t okrow = kFALSE;
274 Bool_t okcol = kFALSE;
275
276 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
277 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
278 iSM,ieta,iphi));
279 }
280
281 //Check rows/phi
282 if(iSM < 10){
283 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
284 }
285 else{
286 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
287 }
288
289 //Check columns/eta
290 if(!fNoEMCALBorderAtEta0){
291 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
292 }
293 else{
294 if(iSM%2==0){
295 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
296 }
297 else {
298 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
299 }
300 }//eta 0 not checked
301
83bfd77a 302 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
fd6df01c 303 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
304
83bfd77a 305 if (okcol && okrow) {
306 //printf("Accept\n");
307 return kTRUE;
308 }
309 else {
310 //printf("Reject\n");
311 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
312 return kFALSE;
313 }
fd6df01c 314
315}
316
317
318//_________________________________________________________________________________________________________
319Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
320 // Check that in the cluster cells, there is no bad channel of those stored
321 // in fEMCALBadChannelMap or fPHOSBadChannelMap
322
323 if(!fRemoveBadChannels) return kFALSE;
324 if(!fEMCALBadChannelMap) return kFALSE;
325
326 Int_t icol = -1;
327 Int_t irow = -1;
328 Int_t imod = -1;
329 for(Int_t iCell = 0; iCell<nCells; iCell++){
330
331 //Get the column and row
332 Int_t iTower = -1, iIphi = -1, iIeta = -1;
333 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
334 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
335 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
83bfd77a 336 if(GetEMCALChannelStatus(imod, icol, irow)){
337 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
338 return kTRUE;
339 }
fd6df01c 340
341 }// cell cluster loop
342
343 return kFALSE;
344
345}
094786cc 346
d9b3567c 347//__________________________________________________
348Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
349// Correct cluster energy from non linearity functions
350 Float_t energy = cluster->E();
351
352 switch (fNonLinearityFunction) {
353
354 case kPi0MC:
871aee7a 355 {
d9b3567c 356 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
871aee7a 357 //Double_t fNonLinearityParams[0] = 1.001;
358 //Double_t fNonLinearityParams[1] = -0.01264;
359 //Double_t fNonLinearityParams[2] = -0.03632;
360 //Double_t fNonLinearityParams[3] = 0.1798;
361 //Double_t fNonLinearityParams[4] = -0.522;
8cdd1f1f 362 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
d9b3567c 363 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
364 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
365 break;
871aee7a 366 }
d9b3567c 367
368 case kPi0GammaGamma:
871aee7a 369 {
d9b3567c 370 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
96957075 371 //Double_t fNonLinearityParams[0] = 1.04;
372 //Double_t fNonLinearityParams[1] = -0.1445;
871aee7a 373 //Double_t fNonLinearityParams[2] = 1.046;
d9b3567c 374 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
375 break;
871aee7a 376 }
d9b3567c 377
378 case kPi0GammaConversion:
871aee7a 379 {
d9b3567c 380 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
871aee7a 381 //fNonLinearityParams[0] = 0.139393/0.1349766;
382 //fNonLinearityParams[1] = 0.0566186;
383 //fNonLinearityParams[2] = 0.982133;
d9b3567c 384 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
385
386 break;
871aee7a 387 }
388
389 case kBeamTest:
390 {
391 //From beam test, Alexei's results, for different ZS thresholds
392 // th=30 MeV; th = 45 MeV; th = 75 MeV
96957075 393 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
871aee7a 394 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
395 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
96957075 396 //Rescale the param[0] with 1.03
871aee7a 397 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
398
399 break;
400 }
4b58ac4f 401 case kBeamTestCorrected:
402 {
403 //From beam test, corrected for material between beam and EMCAL
404 //fNonLinearityParams[0] = 0.983
405 //fNonLinearityParams[1] = 9.83529e-01;
406 //fNonLinearityParams[2] = -1.84235e+02;
407 //fNonLinearityParams[3] = -2.05019e+00;
408 //fNonLinearityParams[4] = -5.89423e+00;
409 energy *= fNonLinearityParams[0]*( ( fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]) ) +
410 ( (fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())) *
411 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3])) ) );
412
413 break;
414 }
d9b3567c 415
416 case kNoCorrection:
417 AliDebug(2,"No correction on the energy\n");
418 break;
419
420 }
421
422 return energy;
423
424}
425
426//__________________________________________________
094786cc 427Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
428{
429 //Calculate shower depth for a given cluster energy and particle type
430
431 // parameters
cb231979 432 Float_t x0 = 1.31;
094786cc 433 Float_t ecr = 8;
434 Float_t depth = 0;
435
436 switch ( iParticle )
437 {
438 case kPhoton:
fd6df01c 439 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 440 break;
441
442 case kElectron:
fd6df01c 443 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 444 break;
445
446 case kHadron:
447 // hadron
448 // boxes anc. here
449 if(gGeoManager){
450 gGeoManager->cd("ALIC_1/XEN1_1");
451 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
452 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
fd6df01c 453 if(geoSM){
454 TGeoVolume *geoSMVol = geoSM->GetVolume();
455 TGeoShape *geoSMShape = geoSMVol->GetShape();
456 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
457 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
458 else AliFatal("Null GEANT box");
459 }else AliFatal("NULL GEANT node matrix");
094786cc 460 }
461 else{//electron
fd6df01c 462 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 463 }
464
465 break;
466
467 default://photon
fd6df01c 468 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 469 }
470
471 return depth;
472
473}
474
475//__________________________________________________
cb231979 476void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
477 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
d9b3567c 478{
479 //For a given CaloCluster gets the absId of the cell
480 //with maximum energy deposit.
481
482 Double_t eMax = -1.;
483 Double_t eCell = -1.;
094786cc 484 Float_t fraction = 1.;
485 Float_t recalFactor = 1.;
d9b3567c 486 Int_t cellAbsId = -1 ;
094786cc 487
d9b3567c 488 Int_t iTower = -1;
489 Int_t iIphi = -1;
490 Int_t iIeta = -1;
cb231979 491 Int_t iSupMod0= -1;
83bfd77a 492 //printf("---Max?\n");
d9b3567c 493 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 494 cellAbsId = clu->GetCellAbsId(iDig);
495 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 496 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 497 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 498 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
499 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
500 if(iDig==0) iSupMod0=iSupMod;
501 else if(iSupMod0!=iSupMod) {
502 shared = kTRUE;
503 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
504 }
094786cc 505 if(IsRecalibrationOn()) {
094786cc 506 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
507 }
508 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 509 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
094786cc 510 if(eCell > eMax) {
d9b3567c 511 eMax = eCell;
512 absId = cellAbsId;
513 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
514 }
515 }// cell loop
516
517 //Get from the absid the supermodule, tower and eta/phi numbers
518 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
519 //Gives SuperModule and Tower numbers
520 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 521 iIphi, iIeta,iphi,ieta);
522 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
523 //printf("Max end---\n");
d9b3567c 524
525}
526
094786cc 527//________________________________________________________________
528void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
529 //Init EMCAL recalibration factors
530 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
531 //In order to avoid rewriting the same histograms
532 Bool_t oldStatus = TH1::AddDirectoryStatus();
533 TH1::AddDirectory(kFALSE);
534
cb231979 535 fEMCALRecalibrationFactors = new TObjArray(10);
094786cc 536 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
537 //Init the histograms with 1
538 for (Int_t sm = 0; sm < 12; sm++) {
539 for (Int_t i = 0; i < 48; i++) {
540 for (Int_t j = 0; j < 24; j++) {
541 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
542 }
543 }
544 }
545 fEMCALRecalibrationFactors->SetOwner(kTRUE);
546 fEMCALRecalibrationFactors->Compress();
547
548 //In order to avoid rewriting the same histograms
549 TH1::AddDirectory(oldStatus);
550}
551
552
fd6df01c 553//________________________________________________________________
554void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
555 //Init EMCAL bad channels map
556 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
557 //In order to avoid rewriting the same histograms
558 Bool_t oldStatus = TH1::AddDirectoryStatus();
559 TH1::AddDirectory(kFALSE);
560
cb231979 561 fEMCALBadChannelMap = new TObjArray(10);
fd6df01c 562 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
6fe0e6d0 563 for (int i = 0; i < 10; i++) {
fd6df01c 564 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
565 }
566
567 //delete hTemp;
568
569 fEMCALBadChannelMap->SetOwner(kTRUE);
570 fEMCALBadChannelMap->Compress();
571
572 //In order to avoid rewriting the same histograms
573 TH1::AddDirectory(oldStatus);
574}
575
094786cc 576//________________________________________________________________
577void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
578 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
579
580 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
581 UShort_t * index = cluster->GetCellsAbsId() ;
582 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
583 Int_t ncells = cluster->GetNCells();
584
585 //Initialize some used variables
586 Float_t energy = 0;
587 Int_t absId = -1;
588 Int_t icol = -1, irow = -1, imod=1;
589 Float_t factor = 1, frac = 0;
590
591 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
592 for(Int_t icell = 0; icell < ncells; icell++){
593 absId = index[icell];
594 frac = fraction[icell];
595 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
596 Int_t iTower = -1, iIphi = -1, iIeta = -1;
597 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
598 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
599 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
600 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
601 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
602 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
603
604 energy += cells->GetCellAmplitude(absId)*factor*frac;
605 }
606
607
608 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
609
610 cluster->SetE(energy);
611
612}
613
614
d9b3567c 615//__________________________________________________
094786cc 616void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 617{
618 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
619
094786cc 620 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
621 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 622 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 623
624}
625
626//__________________________________________________
627void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
628{
629 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
630 // The algorithm is a copy of what is done in AliEMCALRecPoint
631
632 Double_t eCell = 0.;
633 Float_t fraction = 1.;
634 Float_t recalFactor = 1.;
635
636 Int_t absId = -1;
637 Int_t iTower = -1, iIphi = -1, iIeta = -1;
638 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
639 Float_t weight = 0., totalWeight=0.;
640 Float_t newPos[3] = {0,0,0};
641 Double_t pLocal[3], pGlobal[3];
cb231979 642 Bool_t shared = kFALSE;
643
094786cc 644 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
cb231979 645 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 646 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
647
83bfd77a 648 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
649
094786cc 650 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
651 absId = clu->GetCellAbsId(iDig);
652 fraction = clu->GetCellAmplitudeFraction(iDig);
653 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
654 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
655 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
656
657 if(IsRecalibrationOn()) {
658 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
659 }
660 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
661
662 weight = GetCellWeight(eCell,clEnergy);
83bfd77a 663 //printf("cell energy %f, weight %f\n",eCell,weight);
094786cc 664 totalWeight += weight;
665 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 666 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 667 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 668 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
669
094786cc 670 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
671
672 }// cell loop
673
674 if(totalWeight>0){
675 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
676 }
677
094786cc 678 //Float_t pos[]={0,0,0};
679 //clu->GetPosition(pos);
680 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 681 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 682
683 if(iSupModMax > 1) {//sector 1
684 newPos[0] +=fMisalTransShift[3];//-=3.093;
685 newPos[1] +=fMisalTransShift[4];//+=6.82;
686 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 687 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
688
094786cc 689 }
690 else {//sector 0
691 newPos[0] +=fMisalTransShift[0];//+=1.134;
692 newPos[1] +=fMisalTransShift[1];//+=8.2;
693 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 694 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
695
094786cc 696 }
83bfd77a 697 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
698
094786cc 699 clu->SetPosition(newPos);
700
094786cc 701}
702
703//__________________________________________________
704void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
705{
706 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
707 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
708
709 Double_t eCell = 1.;
710 Float_t fraction = 1.;
711 Float_t recalFactor = 1.;
712
713 Int_t absId = -1;
d9b3567c 714 Int_t iTower = -1;
094786cc 715 Int_t iIphi = -1, iIeta = -1;
716 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 717 Int_t iphi = -1, ieta =-1;
cb231979 718 Bool_t shared = kFALSE;
719
d9b3567c 720 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
cb231979 721 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 722 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
723
d9b3567c 724 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 725 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
726 Int_t startingSM = -1;
d9b3567c 727
728 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 729 absId = clu->GetCellAbsId(iDig);
730 fraction = clu->GetCellAmplitudeFraction(iDig);
731 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
732 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
d9b3567c 733 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 734
d9b3567c 735 if (iDig==0) startingSM = iSupMod;
736 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 737
738 eCell = cells->GetCellAmplitude(absId);
d9b3567c 739
094786cc 740 if(IsRecalibrationOn()) {
741 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
742 }
743 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 744
094786cc 745 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 746 if(weight < 0) weight = 0;
747 totalWeight += weight;
748 weightedCol += ieta*weight;
749 weightedRow += iphi*weight;
750
751 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
752
094786cc 753 }// cell loop
754
d9b3567c 755 Float_t xyzNew[]={0.,0.,0.};
756 if(areInSameSM == kTRUE) {
757 //printf("In Same SM\n");
758 weightedCol = weightedCol/totalWeight;
759 weightedRow = weightedRow/totalWeight;
094786cc 760 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 761 }
762 else {
763 //printf("In Different SM\n");
094786cc 764 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 765 }
d9b3567c 766
094786cc 767 clu->SetPosition(xyzNew);
d9b3567c 768
769}
770
cb231979 771//____________________________________________________________________________
772void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
773
774 //re-evaluate distance to bad channel with updated bad map
775
78467229 776 if(!fRecalDistToBadChannels) return;
cb231979 777
778 //Get channels map of the supermodule where the cluster is.
cb231979 779 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
780 Bool_t shared = kFALSE;
781 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
782 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
783
784 Int_t dRrow, dRcol;
785 Float_t minDist = 10000.;
786 Float_t dist = 0.;
787
788 //Loop on tower status map
789 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
790 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
791 //Check if tower is bad.
792 if(hMap->GetBinContent(icol,irow)==0) continue;
793 //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 794 // iSupMod,icol, irow, icolM,irowM);
cb231979 795
796 dRrow=TMath::Abs(irowM-irow);
797 dRcol=TMath::Abs(icolM-icol);
798 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
799 if(dist < minDist){
800 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
801 minDist = dist;
802 }
803
804 }
805 }
806
807 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
808 if (shared) {
809 TH2D* hMap2 = 0;
810 Int_t iSupMod2 = -1;
811
812 //The only possible combinations are (0,1), (2,3) ... (8,9)
813 if(iSupMod%2) iSupMod2 = iSupMod-1;
814 else iSupMod2 = iSupMod+1;
815 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
816
817 //Loop on tower status map of second super module
818 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
819 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
820 //Check if tower is bad.
821 if(hMap2->GetBinContent(icol,irow)==0) continue;
822 //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",
823 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
824
825 dRrow=TMath::Abs(irow-irowM);
826
827 if(iSupMod%2) {
828 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
829 }
830 else {
831 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
832 }
833
834 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
835 if(dist < minDist) minDist = dist;
836
837 }
838 }
839
840 }// shared cluster in 2 SuperModules
78467229 841
6fe0e6d0 842 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
843 cluster->SetDistanceToBadChannel(minDist);
cb231979 844
845}
846
83bfd77a 847//____________________________________________________________________________
848void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
849
850 //re-evaluate identification parameters with bayesian
851
852 if ( cluster->GetM02() != 0)
853 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
854
855 Float_t pidlist[AliPID::kSPECIESN+1];
856 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
857
858 cluster->SetPID(pidlist);
859
860}
861
862//____________________________________________________________________________
863void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
864{
865 // Calculates new center of gravity in the local EMCAL-module coordinates
866 // and tranfers into global ALICE coordinates
867 // Calculates Dispersion and main axis
868
869 Int_t nstat = 0;
870 Float_t wtot = 0. ;
871 Double_t eCell = 0.;
872 Float_t fraction = 1.;
873 Float_t recalFactor = 1.;
874
875 Int_t iSupMod = -1;
876 Int_t iTower = -1;
877 Int_t iIphi = -1;
878 Int_t iIeta = -1;
879 Int_t iphi = -1;
880 Int_t ieta = -1;
881 Double_t etai = -1.;
882 Double_t phii = -1.;
883
884 Double_t w = 0.;
885 Double_t d = 0.;
886 Double_t dxx = 0.;
887 Double_t dzz = 0.;
888 Double_t dxz = 0.;
889 Double_t xmean = 0.;
890 Double_t zmean = 0.;
891
892 //Loop on cells
893 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
894
895 //Get from the absid the supermodule, tower and eta/phi numbers
896 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
897 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
898
899 //Get the cell energy, if recalibration is on, apply factors
900 fraction = cluster->GetCellAmplitudeFraction(iDigit);
901 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
902 if(IsRecalibrationOn()) {
903 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
904 }
905 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
906
907 if(cluster->E() > 0 && eCell > 0){
908
909 w = GetCellWeight(eCell,cluster->E());
910
911 etai=(Double_t)ieta;
912 phii=(Double_t)iphi;
913 if(w > 0.0) {
914 wtot += w ;
915 nstat++;
916 //Shower shape
917 dxx += w * etai * etai ;
918 xmean+= w * etai ;
919 dzz += w * phii * phii ;
920 zmean+= w * phii ;
921 dxz += w * etai * phii ;
922 }
923 }
924 else
925 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
926 }//cell loop
927
928 //Normalize to the weight
929 if (wtot > 0) {
930 xmean /= wtot ;
931 zmean /= wtot ;
932 }
933 else
934 AliError(Form("Wrong weight %f\n", wtot));
935
936 //Calculate dispersion
937 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
938
939 //Get from the absid the supermodule, tower and eta/phi numbers
940 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
941 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
942
943 //Get the cell energy, if recalibration is on, apply factors
944 fraction = cluster->GetCellAmplitudeFraction(iDigit);
945 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
946 if(IsRecalibrationOn()) {
947 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
948 }
949 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
950
951 if(cluster->E() > 0 && eCell > 0){
952
953 w = GetCellWeight(eCell,cluster->E());
954
955 etai=(Double_t)ieta;
956 phii=(Double_t)iphi;
957 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
958 }
959 else
960 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
961 }// cell loop
962
963 //Normalize to the weigth and set shower shape parameters
964 if (wtot > 0 && nstat > 1) {
965 d /= wtot ;
966 dxx /= wtot ;
967 dzz /= wtot ;
968 dxz /= wtot ;
969 dxx -= xmean * xmean ;
970 dzz -= zmean * zmean ;
971 dxz -= xmean * zmean ;
972 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
973 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
974 }
975 else{
976 d=0. ;
977 cluster->SetM20(0.) ;
978 cluster->SetM02(0.) ;
979 }
980
981 if (d>=0)
982 cluster->SetDispersion(TMath::Sqrt(d)) ;
983 else
984 cluster->SetDispersion(0) ;
985
986}
987
b540d03f 988//____________________________________________________________________________
989void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr)
bd8c7aef 990{
991 //This function should be called before the cluster loop
992 //Before call this function, please recalculate the cluster positions
993 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
994 //Store matched cluster indexes and residuals
995 //It only works with ESDs, not AODs
996
b540d03f 997 fMatchedTrackIndex ->Reset();
bd8c7aef 998 fMatchedClusterIndex->Reset();
b540d03f 999 fResidualZ ->Reset();
1000 fResidualR ->Reset();
bd8c7aef 1001
67241b7e 1002 fMatchedTrackIndex ->Set(500);
1003 fMatchedClusterIndex->Set(500);
1004 fResidualZ ->Set(500);
1005 fResidualR ->Set(500);
bd8c7aef 1006
1007 Int_t matched=0;
1008 Float_t clsPos[3];
1009 Double_t trkPos[3];
1010 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1011 {
1012 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
1013 if(!track || !IsAccepted(track)) continue;
1014
1015 Float_t dRMax = fCutR, dZMax = fCutZ;
1016 Int_t index = -1;
1017 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
b540d03f 1018 if(!clusterArr){// get clusters from event
1019 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1020 {
1021 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1022 if(!cluster->IsEMCAL()) continue;
1023 cluster->GetPosition(clsPos); //Has been recalculated
1024 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1025 emctrack->GetXYZ(trkPos);
1026 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1027 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
bd8c7aef 1028
b540d03f 1029 if(tmpR<dRMax)
1030 {
1031 dRMax=tmpR;
1032 dZMax=tmpZ;
1033 index=icl;
1034 }
1035 }//cluster loop
1036 } else { // external cluster array, not from ESD event
1037 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1038 {
1039 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1040 if(!cluster->IsEMCAL()) continue;
1041 cluster->GetPosition(clsPos); //Has been recalculated
1042 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1043 emctrack->GetXYZ(trkPos);
1044 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1045 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1046
1047 if(tmpR<dRMax)
1048 {
1049 dRMax=tmpR;
1050 dZMax=tmpZ;
1051 index=icl;
1052 }
1053 }//cluster loop
1054 }// external list of clusters
1055
bd8c7aef 1056 if(index>-1)
1057 {
b540d03f 1058 fMatchedTrackIndex ->AddAt(itr,matched);
bd8c7aef 1059 fMatchedClusterIndex->AddAt(index,matched);
b540d03f 1060 fResidualZ ->AddAt(dZMax,matched);
1061 fResidualR ->AddAt(dRMax,matched);
bd8c7aef 1062 matched++;
1063 }
1064 delete emctrack;
1065 }//track loop
b540d03f 1066
1067 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1068
1069 fMatchedTrackIndex ->Set(matched);
bd8c7aef 1070 fMatchedClusterIndex->Set(matched);
b540d03f 1071 fResidualZ ->Set(matched);
1072 fResidualR ->Set(matched);
bd8c7aef 1073
1074 //printf("Number of matched pairs: %d\n",matched);
1075}
1076
b540d03f 1077//________________________________________________________________________________
bd8c7aef 1078void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1079{
1080 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1081 //Get the residuals dR and dZ for this cluster
1082 //It only works with ESDs, not AODs
1083
81efb149 1084 if( FindMatchedPos(index) >= 999 )
bd8c7aef 1085 {
1086 AliDebug(2,"No matched tracks found!\n");
1087 dR=999.;
1088 dZ=999.;
1089 return;
1090 }
1091 dR = fResidualR->At(FindMatchedPos(index));
1092 dZ = fResidualZ->At(FindMatchedPos(index));
96957075 1093 //printf("dR %f, dZ %f\n",dR, dZ);
bd8c7aef 1094}
1095
b540d03f 1096//__________________________________________________________
1097Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t index)
1098{
1099 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1100 //Get the index of matched track for this cluster
1101 //It only works with ESDs, not AODs
1102
1103 if(IsMatched(index))
1104 return fMatchedTrackIndex->At(FindMatchedPos(index));
1105 else
1106 return -1;
1107}
1108
1109
bd8c7aef 1110//__________________________________________________
1111Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1112{
1113 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1114 //Returns if cluster has a match
81efb149 1115 if(FindMatchedPos(index) < 999)
82d09e74 1116 return kTRUE;
bd8c7aef 1117 else
1118 return kFALSE;
1119}
b540d03f 1120//__________________________________________________________
81efb149 1121UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
bd8c7aef 1122{
1123 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1124 //Returns the position of the match in the fMatchedClusterIndex array
1125 Float_t tmpR = fCutR;
81efb149 1126 UInt_t pos = 999;
b540d03f 1127
bd8c7aef 1128 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1129 {
1130 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
bd8c7aef 1131 {
b540d03f 1132 pos=i;
1133 tmpR=fResidualR->At(i);
82d09e74 1134 AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
bd8c7aef 1135 }
b540d03f 1136 }
bd8c7aef 1137 return pos;
1138}
1139
b540d03f 1140//__________________________________________________________
bd8c7aef 1141Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1142{
1143 // Given a esd track, return whether the track survive all the cuts
1144
1145 // The different quality parameter are first
1146 // retrieved from the track. then it is found out what cuts the
1147 // track did not survive and finally the cuts are imposed.
1148
1149 UInt_t status = esdTrack->GetStatus();
1150
1151 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1152 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1153
1154 Float_t chi2PerClusterITS = -1;
1155 Float_t chi2PerClusterTPC = -1;
1156 if (nClustersITS!=0)
1157 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1158 if (nClustersTPC!=0)
1159 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1160
1161
1162 //DCA cuts
1163 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1164 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1165 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1166
1167
bd8c7aef 1168 Float_t b[2];
1169 Float_t bCov[3];
1170 esdTrack->GetImpactParameters(b,bCov);
1171 if (bCov[0]<=0 || bCov[2]<=0) {
1172 AliDebug(1, "Estimated b resolution lower or equal zero!");
1173 bCov[0]=0; bCov[2]=0;
1174 }
1175
1176 Float_t dcaToVertexXY = b[0];
1177 Float_t dcaToVertexZ = b[1];
1178 Float_t dcaToVertex = -1;
1179
1180 if (fCutDCAToVertex2D)
1181 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1182 else
1183 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1184
1185 // cut the track?
1186
1187 Bool_t cuts[kNCuts];
1188 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1189
1190 // track quality cuts
1191 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1192 cuts[0]=kTRUE;
1193 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1194 cuts[1]=kTRUE;
1195 if (nClustersTPC<fCutMinNClusterTPC)
1196 cuts[2]=kTRUE;
1197 if (nClustersITS<fCutMinNClusterITS)
1198 cuts[3]=kTRUE;
1199 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1200 cuts[4]=kTRUE;
1201 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1202 cuts[5]=kTRUE;
1203 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1204 cuts[6]=kTRUE;
1205 if (fCutDCAToVertex2D && dcaToVertex > 1)
1206 cuts[7] = kTRUE;
1207 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1208 cuts[8] = kTRUE;
1209 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1210 cuts[9] = kTRUE;
1211
82d09e74 1212 //Require at least one SPD point + anything else in ITS
1213 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1214 cuts[10] = kTRUE;
1215
bd8c7aef 1216 Bool_t cut=kFALSE;
1217 for (Int_t i=0; i<kNCuts; i++)
1218 if (cuts[i]) {cut = kTRUE;}
1219
1220 // cut the track
1221 if (cut)
1222 return kFALSE;
1223 else
1224 return kTRUE;
1225}
1226//__________________________________________________
1227void AliEMCALRecoUtils::InitTrackCuts()
1228{
1229 //Intilize the track cut criteria
82d09e74 1230 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1231 //Also you can customize the cuts using the setters
82d09e74 1232
1233 //TPC
1234 SetMinNClustersTPC(70);
bd8c7aef 1235 SetMaxChi2PerClusterTPC(4);
1236 SetAcceptKinkDaughters(kFALSE);
82d09e74 1237 SetRequireTPCRefit(kTRUE);
1238
1239 //ITS
1240 SetRequireITSRefit(kTRUE);
1241 SetMaxDCAToVertexZ(2);
1242 SetDCAToVertex2D(kFALSE);
171be15c 1243 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1244 SetMinNClustersITS();
bd8c7aef 1245}
83bfd77a 1246
b540d03f 1247//___________________________________________________
d9b3567c 1248void AliEMCALRecoUtils::Print(const Option_t *) const
1249{
1250 // Print Parameters
1251
1252 printf("AliEMCALRecoUtils Settings: \n");
1253 printf("Misalignment shifts\n");
2a71e873 1254 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,
1255 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1256 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1257 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1258 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1259
1260 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1261
1262 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1263
1264 printf("Track cuts: \n");
1265 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1266 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1267 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1268 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1269 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1270
d9b3567c 1271}
96957075 1272
b540d03f 1273//_____________________________________________________________________
96957075 1274void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1275 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1276 //Do it only once and only if it is requested
1277
1278 if(!fUseTimeCorrectionFactors) return;
1279 if(fTimeCorrectionFactorsSet) return;
1280
1281 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1282
1283 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1284 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1285
1286 SwitchOnRecalibration();
1287 for(Int_t ism = 0; ism < 4; ism++){
1288 for(Int_t icol = 0; icol < 48; icol++){
1289 for(Int_t irow = 0; irow < 24; irow++){
1290 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1291 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1292 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1293 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1294 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1295 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1296 }
1297 }
1298 }
1299 fTimeCorrectionFactorsSet = kTRUE;
1300}
1301