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