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