cleaning up HAVE_NOT_ALTRORAWSTREAMV3 preprocessor define
[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
232//__________________________________________________
094786cc 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
347//__________________________________________________
d9b3567c 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 }
d9b3567c 401
402 case kNoCorrection:
403 AliDebug(2,"No correction on the energy\n");
404 break;
405
406 }
407
408 return energy;
409
410}
411
412//__________________________________________________
094786cc 413Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
414{
415 //Calculate shower depth for a given cluster energy and particle type
416
417 // parameters
cb231979 418 Float_t x0 = 1.31;
094786cc 419 Float_t ecr = 8;
420 Float_t depth = 0;
421
422 switch ( iParticle )
423 {
424 case kPhoton:
fd6df01c 425 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 426 break;
427
428 case kElectron:
fd6df01c 429 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 430 break;
431
432 case kHadron:
433 // hadron
434 // boxes anc. here
435 if(gGeoManager){
436 gGeoManager->cd("ALIC_1/XEN1_1");
437 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
438 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
fd6df01c 439 if(geoSM){
440 TGeoVolume *geoSMVol = geoSM->GetVolume();
441 TGeoShape *geoSMShape = geoSMVol->GetShape();
442 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
443 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
444 else AliFatal("Null GEANT box");
445 }else AliFatal("NULL GEANT node matrix");
094786cc 446 }
447 else{//electron
fd6df01c 448 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 449 }
450
451 break;
452
453 default://photon
fd6df01c 454 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 455 }
456
457 return depth;
458
459}
460
461//__________________________________________________
cb231979 462void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
463 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
d9b3567c 464{
465 //For a given CaloCluster gets the absId of the cell
466 //with maximum energy deposit.
467
468 Double_t eMax = -1.;
469 Double_t eCell = -1.;
094786cc 470 Float_t fraction = 1.;
471 Float_t recalFactor = 1.;
d9b3567c 472 Int_t cellAbsId = -1 ;
094786cc 473
d9b3567c 474 Int_t iTower = -1;
475 Int_t iIphi = -1;
476 Int_t iIeta = -1;
cb231979 477 Int_t iSupMod0= -1;
83bfd77a 478 //printf("---Max?\n");
d9b3567c 479 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 480 cellAbsId = clu->GetCellAbsId(iDig);
481 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 482 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 483 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 484 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
485 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
486 if(iDig==0) iSupMod0=iSupMod;
487 else if(iSupMod0!=iSupMod) {
488 shared = kTRUE;
489 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
490 }
094786cc 491 if(IsRecalibrationOn()) {
094786cc 492 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
493 }
494 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 495 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
094786cc 496 if(eCell > eMax) {
d9b3567c 497 eMax = eCell;
498 absId = cellAbsId;
499 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
500 }
501 }// cell loop
502
503 //Get from the absid the supermodule, tower and eta/phi numbers
504 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
505 //Gives SuperModule and Tower numbers
506 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 507 iIphi, iIeta,iphi,ieta);
508 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
509 //printf("Max end---\n");
d9b3567c 510
511}
512
094786cc 513//________________________________________________________________
514void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
515 //Init EMCAL recalibration factors
516 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
517 //In order to avoid rewriting the same histograms
518 Bool_t oldStatus = TH1::AddDirectoryStatus();
519 TH1::AddDirectory(kFALSE);
520
cb231979 521 fEMCALRecalibrationFactors = new TObjArray(10);
094786cc 522 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));
523 //Init the histograms with 1
524 for (Int_t sm = 0; sm < 12; sm++) {
525 for (Int_t i = 0; i < 48; i++) {
526 for (Int_t j = 0; j < 24; j++) {
527 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
528 }
529 }
530 }
531 fEMCALRecalibrationFactors->SetOwner(kTRUE);
532 fEMCALRecalibrationFactors->Compress();
533
534 //In order to avoid rewriting the same histograms
535 TH1::AddDirectory(oldStatus);
536}
537
538
539//________________________________________________________________
fd6df01c 540void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
541 //Init EMCAL bad channels map
542 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
543 //In order to avoid rewriting the same histograms
544 Bool_t oldStatus = TH1::AddDirectoryStatus();
545 TH1::AddDirectory(kFALSE);
546
cb231979 547 fEMCALBadChannelMap = new TObjArray(10);
fd6df01c 548 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
6fe0e6d0 549 for (int i = 0; i < 10; i++) {
fd6df01c 550 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
551 }
552
553 //delete hTemp;
554
555 fEMCALBadChannelMap->SetOwner(kTRUE);
556 fEMCALBadChannelMap->Compress();
557
558 //In order to avoid rewriting the same histograms
559 TH1::AddDirectory(oldStatus);
560}
561
562//________________________________________________________________
094786cc 563void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
564 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
565
566 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
567 UShort_t * index = cluster->GetCellsAbsId() ;
568 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
569 Int_t ncells = cluster->GetNCells();
570
571 //Initialize some used variables
572 Float_t energy = 0;
573 Int_t absId = -1;
574 Int_t icol = -1, irow = -1, imod=1;
575 Float_t factor = 1, frac = 0;
576
577 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
578 for(Int_t icell = 0; icell < ncells; icell++){
579 absId = index[icell];
580 frac = fraction[icell];
581 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
582 Int_t iTower = -1, iIphi = -1, iIeta = -1;
583 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
584 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
585 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
586 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
587 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
588 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
589
590 energy += cells->GetCellAmplitude(absId)*factor*frac;
591 }
592
593
594 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
595
596 cluster->SetE(energy);
597
598}
599
600
d9b3567c 601//__________________________________________________
094786cc 602void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 603{
604 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
605
094786cc 606 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
607 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 608 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 609
610}
611
612//__________________________________________________
613void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
614{
615 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
616 // The algorithm is a copy of what is done in AliEMCALRecPoint
617
618 Double_t eCell = 0.;
619 Float_t fraction = 1.;
620 Float_t recalFactor = 1.;
621
622 Int_t absId = -1;
623 Int_t iTower = -1, iIphi = -1, iIeta = -1;
624 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
625 Float_t weight = 0., totalWeight=0.;
626 Float_t newPos[3] = {0,0,0};
627 Double_t pLocal[3], pGlobal[3];
cb231979 628 Bool_t shared = kFALSE;
629
094786cc 630 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
cb231979 631 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 632 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
633
83bfd77a 634 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
635
094786cc 636 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
637 absId = clu->GetCellAbsId(iDig);
638 fraction = clu->GetCellAmplitudeFraction(iDig);
639 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
640 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
641 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
642
643 if(IsRecalibrationOn()) {
644 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
645 }
646 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
647
648 weight = GetCellWeight(eCell,clEnergy);
83bfd77a 649 //printf("cell energy %f, weight %f\n",eCell,weight);
094786cc 650 totalWeight += weight;
651 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 652 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 653 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 654 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
655
094786cc 656 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
657
658 }// cell loop
659
660 if(totalWeight>0){
661 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
662 }
663
094786cc 664 //Float_t pos[]={0,0,0};
665 //clu->GetPosition(pos);
666 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 667 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 668
669 if(iSupModMax > 1) {//sector 1
670 newPos[0] +=fMisalTransShift[3];//-=3.093;
671 newPos[1] +=fMisalTransShift[4];//+=6.82;
672 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 673 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
674
094786cc 675 }
676 else {//sector 0
677 newPos[0] +=fMisalTransShift[0];//+=1.134;
678 newPos[1] +=fMisalTransShift[1];//+=8.2;
679 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 680 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
681
094786cc 682 }
83bfd77a 683 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
684
094786cc 685 clu->SetPosition(newPos);
686
094786cc 687}
688
689//__________________________________________________
690void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
691{
692 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
693 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
694
695 Double_t eCell = 1.;
696 Float_t fraction = 1.;
697 Float_t recalFactor = 1.;
698
699 Int_t absId = -1;
d9b3567c 700 Int_t iTower = -1;
094786cc 701 Int_t iIphi = -1, iIeta = -1;
702 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 703 Int_t iphi = -1, ieta =-1;
cb231979 704 Bool_t shared = kFALSE;
705
d9b3567c 706 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
cb231979 707 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 708 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
709
d9b3567c 710 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 711 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
712 Int_t startingSM = -1;
d9b3567c 713
714 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 715 absId = clu->GetCellAbsId(iDig);
716 fraction = clu->GetCellAmplitudeFraction(iDig);
717 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
718 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
d9b3567c 719 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 720
d9b3567c 721 if (iDig==0) startingSM = iSupMod;
722 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 723
724 eCell = cells->GetCellAmplitude(absId);
d9b3567c 725
094786cc 726 if(IsRecalibrationOn()) {
727 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
728 }
729 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 730
094786cc 731 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 732 if(weight < 0) weight = 0;
733 totalWeight += weight;
734 weightedCol += ieta*weight;
735 weightedRow += iphi*weight;
736
737 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
738
094786cc 739 }// cell loop
740
d9b3567c 741 Float_t xyzNew[]={0.,0.,0.};
742 if(areInSameSM == kTRUE) {
743 //printf("In Same SM\n");
744 weightedCol = weightedCol/totalWeight;
745 weightedRow = weightedRow/totalWeight;
094786cc 746 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 747 }
748 else {
749 //printf("In Different SM\n");
094786cc 750 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 751 }
d9b3567c 752
094786cc 753 clu->SetPosition(xyzNew);
d9b3567c 754
755}
756
83bfd77a 757//____________________________________________________________________________
cb231979 758void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
759
760 //re-evaluate distance to bad channel with updated bad map
761
78467229 762 if(!fRecalDistToBadChannels) return;
cb231979 763
764 //Get channels map of the supermodule where the cluster is.
cb231979 765 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
766 Bool_t shared = kFALSE;
767 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
768 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
769
770 Int_t dRrow, dRcol;
771 Float_t minDist = 10000.;
772 Float_t dist = 0.;
773
774 //Loop on tower status map
775 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
776 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
777 //Check if tower is bad.
778 if(hMap->GetBinContent(icol,irow)==0) continue;
779 //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 780 // iSupMod,icol, irow, icolM,irowM);
cb231979 781
782 dRrow=TMath::Abs(irowM-irow);
783 dRcol=TMath::Abs(icolM-icol);
784 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
785 if(dist < minDist){
786 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
787 minDist = dist;
788 }
789
790 }
791 }
792
793 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
794 if (shared) {
795 TH2D* hMap2 = 0;
796 Int_t iSupMod2 = -1;
797
798 //The only possible combinations are (0,1), (2,3) ... (8,9)
799 if(iSupMod%2) iSupMod2 = iSupMod-1;
800 else iSupMod2 = iSupMod+1;
801 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
802
803 //Loop on tower status map of second super module
804 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
805 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
806 //Check if tower is bad.
807 if(hMap2->GetBinContent(icol,irow)==0) continue;
808 //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",
809 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
810
811 dRrow=TMath::Abs(irow-irowM);
812
813 if(iSupMod%2) {
814 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
815 }
816 else {
817 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
818 }
819
820 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
821 if(dist < minDist) minDist = dist;
822
823 }
824 }
825
826 }// shared cluster in 2 SuperModules
78467229 827
6fe0e6d0 828 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
829 cluster->SetDistanceToBadChannel(minDist);
cb231979 830
831}
832
833//____________________________________________________________________________
83bfd77a 834void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
835
836 //re-evaluate identification parameters with bayesian
837
838 if ( cluster->GetM02() != 0)
839 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
840
841 Float_t pidlist[AliPID::kSPECIESN+1];
842 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
843
844 cluster->SetPID(pidlist);
845
846}
847
848//____________________________________________________________________________
849void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
850{
851 // Calculates new center of gravity in the local EMCAL-module coordinates
852 // and tranfers into global ALICE coordinates
853 // Calculates Dispersion and main axis
854
855 Int_t nstat = 0;
856 Float_t wtot = 0. ;
857 Double_t eCell = 0.;
858 Float_t fraction = 1.;
859 Float_t recalFactor = 1.;
860
861 Int_t iSupMod = -1;
862 Int_t iTower = -1;
863 Int_t iIphi = -1;
864 Int_t iIeta = -1;
865 Int_t iphi = -1;
866 Int_t ieta = -1;
867 Double_t etai = -1.;
868 Double_t phii = -1.;
869
870 Double_t w = 0.;
871 Double_t d = 0.;
872 Double_t dxx = 0.;
873 Double_t dzz = 0.;
874 Double_t dxz = 0.;
875 Double_t xmean = 0.;
876 Double_t zmean = 0.;
877
878 //Loop on cells
879 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
880
881 //Get from the absid the supermodule, tower and eta/phi numbers
882 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
883 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
884
885 //Get the cell energy, if recalibration is on, apply factors
886 fraction = cluster->GetCellAmplitudeFraction(iDigit);
887 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
888 if(IsRecalibrationOn()) {
889 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
890 }
891 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
892
893 if(cluster->E() > 0 && eCell > 0){
894
895 w = GetCellWeight(eCell,cluster->E());
896
897 etai=(Double_t)ieta;
898 phii=(Double_t)iphi;
899 if(w > 0.0) {
900 wtot += w ;
901 nstat++;
902 //Shower shape
903 dxx += w * etai * etai ;
904 xmean+= w * etai ;
905 dzz += w * phii * phii ;
906 zmean+= w * phii ;
907 dxz += w * etai * phii ;
908 }
909 }
910 else
911 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
912 }//cell loop
913
914 //Normalize to the weight
915 if (wtot > 0) {
916 xmean /= wtot ;
917 zmean /= wtot ;
918 }
919 else
920 AliError(Form("Wrong weight %f\n", wtot));
921
922 //Calculate dispersion
923 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
924
925 //Get from the absid the supermodule, tower and eta/phi numbers
926 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
927 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
928
929 //Get the cell energy, if recalibration is on, apply factors
930 fraction = cluster->GetCellAmplitudeFraction(iDigit);
931 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
932 if(IsRecalibrationOn()) {
933 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
934 }
935 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
936
937 if(cluster->E() > 0 && eCell > 0){
938
939 w = GetCellWeight(eCell,cluster->E());
940
941 etai=(Double_t)ieta;
942 phii=(Double_t)iphi;
943 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
944 }
945 else
946 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
947 }// cell loop
948
949 //Normalize to the weigth and set shower shape parameters
950 if (wtot > 0 && nstat > 1) {
951 d /= wtot ;
952 dxx /= wtot ;
953 dzz /= wtot ;
954 dxz /= wtot ;
955 dxx -= xmean * xmean ;
956 dzz -= zmean * zmean ;
957 dxz -= xmean * zmean ;
958 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
959 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
960 }
961 else{
962 d=0. ;
963 cluster->SetM20(0.) ;
964 cluster->SetM02(0.) ;
965 }
966
967 if (d>=0)
968 cluster->SetDispersion(TMath::Sqrt(d)) ;
969 else
970 cluster->SetDispersion(0) ;
971
972}
973
b540d03f 974//____________________________________________________________________________
975void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr)
bd8c7aef 976{
977 //This function should be called before the cluster loop
978 //Before call this function, please recalculate the cluster positions
979 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
980 //Store matched cluster indexes and residuals
981 //It only works with ESDs, not AODs
982
b540d03f 983 fMatchedTrackIndex ->Reset();
bd8c7aef 984 fMatchedClusterIndex->Reset();
b540d03f 985 fResidualZ ->Reset();
986 fResidualR ->Reset();
bd8c7aef 987
b540d03f 988 fMatchedTrackIndex ->Set(100);
bd8c7aef 989 fMatchedClusterIndex->Set(100);
b540d03f 990 fResidualZ ->Set(100);
991 fResidualR ->Set(100);
bd8c7aef 992
993 Int_t matched=0;
994 Float_t clsPos[3];
995 Double_t trkPos[3];
996 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
997 {
998 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
999 if(!track || !IsAccepted(track)) continue;
1000
1001 Float_t dRMax = fCutR, dZMax = fCutZ;
1002 Int_t index = -1;
1003 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
b540d03f 1004 if(!clusterArr){// get clusters from event
1005 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1006 {
1007 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1008 if(!cluster->IsEMCAL()) continue;
1009 cluster->GetPosition(clsPos); //Has been recalculated
1010 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1011 emctrack->GetXYZ(trkPos);
1012 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) );
1013 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
bd8c7aef 1014
b540d03f 1015 if(tmpR<dRMax)
1016 {
1017 dRMax=tmpR;
1018 dZMax=tmpZ;
1019 index=icl;
1020 }
1021 }//cluster loop
1022 } else { // external cluster array, not from ESD event
1023 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1024 {
1025 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1026 if(!cluster->IsEMCAL()) continue;
1027 cluster->GetPosition(clsPos); //Has been recalculated
1028 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1029 emctrack->GetXYZ(trkPos);
1030 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) );
1031 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1032
1033 if(tmpR<dRMax)
1034 {
1035 dRMax=tmpR;
1036 dZMax=tmpZ;
1037 index=icl;
1038 }
1039 }//cluster loop
1040 }// external list of clusters
1041
bd8c7aef 1042 if(index>-1)
1043 {
b540d03f 1044 fMatchedTrackIndex ->AddAt(itr,matched);
bd8c7aef 1045 fMatchedClusterIndex->AddAt(index,matched);
b540d03f 1046 fResidualZ ->AddAt(dZMax,matched);
1047 fResidualR ->AddAt(dRMax,matched);
bd8c7aef 1048 matched++;
1049 }
1050 delete emctrack;
1051 }//track loop
b540d03f 1052
1053 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1054
1055 fMatchedTrackIndex ->Set(matched);
bd8c7aef 1056 fMatchedClusterIndex->Set(matched);
b540d03f 1057 fResidualZ ->Set(matched);
1058 fResidualR ->Set(matched);
bd8c7aef 1059
1060 //printf("Number of matched pairs: %d\n",matched);
1061}
1062
b540d03f 1063//________________________________________________________________________________
bd8c7aef 1064void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1065{
1066 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1067 //Get the residuals dR and dZ for this cluster
1068 //It only works with ESDs, not AODs
1069
81efb149 1070 if( FindMatchedPos(index) >= 999 )
bd8c7aef 1071 {
1072 AliDebug(2,"No matched tracks found!\n");
1073 dR=999.;
1074 dZ=999.;
1075 return;
1076 }
1077 dR = fResidualR->At(FindMatchedPos(index));
1078 dZ = fResidualZ->At(FindMatchedPos(index));
96957075 1079 //printf("dR %f, dZ %f\n",dR, dZ);
bd8c7aef 1080}
1081
b540d03f 1082//__________________________________________________________
1083Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t index)
1084{
1085 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1086 //Get the index of matched track for this cluster
1087 //It only works with ESDs, not AODs
1088
1089 if(IsMatched(index))
1090 return fMatchedTrackIndex->At(FindMatchedPos(index));
1091 else
1092 return -1;
1093}
1094
1095
bd8c7aef 1096//__________________________________________________
1097Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1098{
1099 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1100 //Returns if cluster has a match
81efb149 1101 if(FindMatchedPos(index) < 999)
82d09e74 1102 return kTRUE;
bd8c7aef 1103 else
1104 return kFALSE;
1105}
b540d03f 1106//__________________________________________________________
81efb149 1107UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
bd8c7aef 1108{
1109 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1110 //Returns the position of the match in the fMatchedClusterIndex array
1111 Float_t tmpR = fCutR;
81efb149 1112 UInt_t pos = 999;
b540d03f 1113
bd8c7aef 1114 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1115 {
1116 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
bd8c7aef 1117 {
b540d03f 1118 pos=i;
1119 tmpR=fResidualR->At(i);
82d09e74 1120 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 1121 }
b540d03f 1122 }
bd8c7aef 1123 return pos;
1124}
1125
b540d03f 1126//__________________________________________________________
bd8c7aef 1127Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1128{
1129 // Given a esd track, return whether the track survive all the cuts
1130
1131 // The different quality parameter are first
1132 // retrieved from the track. then it is found out what cuts the
1133 // track did not survive and finally the cuts are imposed.
1134
1135 UInt_t status = esdTrack->GetStatus();
1136
1137 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1138 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1139
1140 Float_t chi2PerClusterITS = -1;
1141 Float_t chi2PerClusterTPC = -1;
1142 if (nClustersITS!=0)
1143 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1144 if (nClustersTPC!=0)
1145 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1146
1147
1148 //DCA cuts
1149 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1150 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1151 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1152
1153
bd8c7aef 1154 Float_t b[2];
1155 Float_t bCov[3];
1156 esdTrack->GetImpactParameters(b,bCov);
1157 if (bCov[0]<=0 || bCov[2]<=0) {
1158 AliDebug(1, "Estimated b resolution lower or equal zero!");
1159 bCov[0]=0; bCov[2]=0;
1160 }
1161
1162 Float_t dcaToVertexXY = b[0];
1163 Float_t dcaToVertexZ = b[1];
1164 Float_t dcaToVertex = -1;
1165
1166 if (fCutDCAToVertex2D)
1167 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1168 else
1169 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1170
1171 // cut the track?
1172
1173 Bool_t cuts[kNCuts];
1174 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1175
1176 // track quality cuts
1177 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1178 cuts[0]=kTRUE;
1179 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1180 cuts[1]=kTRUE;
1181 if (nClustersTPC<fCutMinNClusterTPC)
1182 cuts[2]=kTRUE;
1183 if (nClustersITS<fCutMinNClusterITS)
1184 cuts[3]=kTRUE;
1185 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1186 cuts[4]=kTRUE;
1187 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1188 cuts[5]=kTRUE;
1189 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1190 cuts[6]=kTRUE;
1191 if (fCutDCAToVertex2D && dcaToVertex > 1)
1192 cuts[7] = kTRUE;
1193 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1194 cuts[8] = kTRUE;
1195 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1196 cuts[9] = kTRUE;
1197
82d09e74 1198 //Require at least one SPD point + anything else in ITS
1199 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1200 cuts[10] = kTRUE;
1201
bd8c7aef 1202 Bool_t cut=kFALSE;
1203 for (Int_t i=0; i<kNCuts; i++)
1204 if (cuts[i]) {cut = kTRUE;}
1205
1206 // cut the track
1207 if (cut)
1208 return kFALSE;
1209 else
1210 return kTRUE;
1211}
1212//__________________________________________________
1213void AliEMCALRecoUtils::InitTrackCuts()
1214{
1215 //Intilize the track cut criteria
82d09e74 1216 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1217 //Also you can customize the cuts using the setters
82d09e74 1218
1219 //TPC
1220 SetMinNClustersTPC(70);
bd8c7aef 1221 SetMaxChi2PerClusterTPC(4);
1222 SetAcceptKinkDaughters(kFALSE);
82d09e74 1223 SetRequireTPCRefit(kTRUE);
1224
1225 //ITS
1226 SetRequireITSRefit(kTRUE);
1227 SetMaxDCAToVertexZ(2);
1228 SetDCAToVertex2D(kFALSE);
171be15c 1229 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1230 SetMinNClustersITS();
bd8c7aef 1231}
83bfd77a 1232
b540d03f 1233//___________________________________________________
d9b3567c 1234void AliEMCALRecoUtils::Print(const Option_t *) const
1235{
1236 // Print Parameters
1237
1238 printf("AliEMCALRecoUtils Settings: \n");
1239 printf("Misalignment shifts\n");
2a71e873 1240 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,
1241 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1242 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1243 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1244 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1245
1246 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1247
1248 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1249
1250 printf("Track cuts: \n");
1251 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1252 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1253 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1254 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1255 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1256
d9b3567c 1257}
96957075 1258
b540d03f 1259//_____________________________________________________________________
96957075 1260void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1261 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1262 //Do it only once and only if it is requested
1263
1264 if(!fUseTimeCorrectionFactors) return;
1265 if(fTimeCorrectionFactorsSet) return;
1266
1267 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1268
1269 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1270 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1271
1272 SwitchOnRecalibration();
1273 for(Int_t ism = 0; ism < 4; ism++){
1274 for(Int_t icol = 0; icol < 48; icol++){
1275 for(Int_t irow = 0; irow < 24; irow++){
1276 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1277 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1278 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1279 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1280 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1281 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1282 }
1283 }
1284 }
1285 fTimeCorrectionFactorsSet = kTRUE;
1286}
1287