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