1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
18 ///////////////////////////////////////////////////////////////////////////////
20 // Class AliEMCALRecoUtils
21 // Some utilities to recalculate the cluster position or energy linearity
24 // Author: Gustavo Conesa (LPSC- Grenoble)
25 // Track matching part: Rongrong Ma (Yale)
27 ///////////////////////////////////////////////////////////////////////////////
31 // standard C++ includes
32 //#include <Riostream.h>
35 #include <TGeoManager.h>
36 #include <TGeoMatrix.h>
40 #include "AliVCluster.h"
41 #include "AliVCaloCells.h"
42 #include "AliVEvent.h"
45 #include "AliESDEvent.h"
46 #include "AliESDtrack.h"
49 #include "AliEMCALRecoUtils.h"
50 #include "AliEMCALGeometry.h"
51 #include "AliEMCALTrack.h"
52 #include "AliEMCALCalibTimeDepCorrection.h"
53 #include "AliEMCALPIDUtils.h"
55 ClassImp(AliEMCALRecoUtils)
57 //______________________________________________
58 AliEMCALRecoUtils::AliEMCALRecoUtils():
59 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
60 fPosAlgo(kUnchanged),fW0(4.),
61 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
62 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
63 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
64 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
65 fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
66 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
67 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
68 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
69 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE)
73 // Initialize all constant values which have to be used
74 // during Reco algorithm execution
77 //Misalignment matrices
78 for(Int_t i = 0; i < 15 ; i++) {
79 fMisalTransShift[i] = 0.;
80 fMisalRotShift[i] = 0.;
84 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
85 //For kPi0GammaGamma case, but default is no correction
86 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
87 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
88 fNonLinearityParams[2] = 1.046;
91 fMatchedTrackIndex = new TArrayI();
92 fMatchedClusterIndex = new TArrayI();
93 fResidualZ = new TArrayF();
94 fResidualR = new TArrayF();
98 fPIDUtils = new AliEMCALPIDUtils();
103 //______________________________________________________________________
104 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
105 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
106 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
107 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
108 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
109 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
110 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
111 fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
112 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
113 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
114 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
115 fCutR(reco.fCutR),fCutZ(reco.fCutZ),
116 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
117 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
118 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
119 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
120 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
121 fPIDUtils(reco.fPIDUtils),
122 fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
126 for(Int_t i = 0; i < 15 ; i++) {
127 fMisalRotShift[i] = reco.fMisalRotShift[i];
128 fMisalTransShift[i] = reco.fMisalTransShift[i];
130 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
135 //______________________________________________________________________
136 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
138 //Assignment operator
140 if(this == &reco)return *this;
141 ((TNamed *)this)->operator=(reco);
143 fNonLinearityFunction = reco.fNonLinearityFunction;
144 fParticleType = reco.fParticleType;
145 fPosAlgo = reco.fPosAlgo;
147 fRecalibration = reco.fRecalibration;
148 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
149 fRemoveBadChannels = reco.fRemoveBadChannels;
150 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
151 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
152 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
153 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
156 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
157 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
162 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
163 fCutMinNClusterITS = reco.fCutMinNClusterITS;
164 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
165 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
166 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
167 fCutRequireITSRefit = reco.fCutRequireITSRefit;
168 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
169 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
170 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
171 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
173 fPIDUtils = reco.fPIDUtils;
175 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
176 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
180 // assign or copy construct
182 *fResidualR = *reco.fResidualR;
184 else fResidualR = new TArrayF(*reco.fResidualR);
187 if(fResidualR)delete fResidualR;
192 // assign or copy construct
194 *fResidualZ = *reco.fResidualZ;
196 else fResidualZ = new TArrayF(*reco.fResidualZ);
199 if(fResidualZ)delete fResidualZ;
203 if(reco.fMatchedTrackIndex){
204 // assign or copy construct
205 if(fMatchedTrackIndex){
206 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
208 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
211 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
212 fMatchedTrackIndex = 0;
215 if(reco.fMatchedClusterIndex){
216 // assign or copy construct
217 if(fMatchedClusterIndex){
218 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
220 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
223 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
224 fMatchedClusterIndex = 0;
232 //__________________________________________________
233 AliEMCALRecoUtils::~AliEMCALRecoUtils()
237 if(fEMCALRecalibrationFactors) {
238 fEMCALRecalibrationFactors->Clear();
239 delete fEMCALRecalibrationFactors;
242 if(fEMCALBadChannelMap) {
243 fEMCALBadChannelMap->Clear();
244 delete fEMCALBadChannelMap;
247 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
248 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
249 if(fResidualR) {delete fResidualR; fResidualR=0;}
250 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
254 //_______________________________________________________________
255 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
257 // Given the list of AbsId of the cluster, get the maximum cell and
258 // check if there are fNCellsFromBorder from the calorimeter border
260 //If the distance to the border is 0 or negative just exit accept all clusters
261 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
263 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
264 Bool_t shared = kFALSE;
265 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
267 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
268 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
270 if(absIdMax==-1) return kFALSE;
272 //Check if the cell is close to the borders:
273 Bool_t okrow = kFALSE;
274 Bool_t okcol = kFALSE;
276 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
277 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
283 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
286 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
290 if(!fNoEMCALBorderAtEta0){
291 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
295 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
298 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
302 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
303 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
305 if (okcol && okrow) {
306 //printf("Accept\n");
310 //printf("Reject\n");
311 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
318 //_________________________________________________________________________________________________________
319 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
320 // Check that in the cluster cells, there is no bad channel of those stored
321 // in fEMCALBadChannelMap or fPHOSBadChannelMap
323 if(!fRemoveBadChannels) return kFALSE;
324 if(!fEMCALBadChannelMap) return kFALSE;
329 for(Int_t iCell = 0; iCell<nCells; iCell++){
331 //Get the column and row
332 Int_t iTower = -1, iIphi = -1, iIeta = -1;
333 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
334 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
335 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
336 if(GetEMCALChannelStatus(imod, icol, irow)){
337 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
341 }// cell cluster loop
347 //__________________________________________________
348 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
349 // Correct cluster energy from non linearity functions
350 Float_t energy = cluster->E();
352 switch (fNonLinearityFunction) {
356 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
357 //Double_t fNonLinearityParams[0] = 1.001;
358 //Double_t fNonLinearityParams[1] = -0.01264;
359 //Double_t fNonLinearityParams[2] = -0.03632;
360 //Double_t fNonLinearityParams[3] = 0.1798;
361 //Double_t fNonLinearityParams[4] = -0.522;
362 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
363 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
364 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
370 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
371 //Double_t fNonLinearityParams[0] = 1.04;
372 //Double_t fNonLinearityParams[1] = -0.1445;
373 //Double_t fNonLinearityParams[2] = 1.046;
374 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
378 case kPi0GammaConversion:
380 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
381 //fNonLinearityParams[0] = 0.139393/0.1349766;
382 //fNonLinearityParams[1] = 0.0566186;
383 //fNonLinearityParams[2] = 0.982133;
384 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
391 //From beam test, Alexei's results, for different ZS thresholds
392 // th=30 MeV; th = 45 MeV; th = 75 MeV
393 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
394 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
395 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
396 //Rescale the param[0] with 1.03
397 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
401 case kBeamTestCorrected:
403 //From beam test, corrected for material between beam and EMCAL
404 //fNonLinearityParams[0] = 0.983
405 //fNonLinearityParams[1] = 9.83529e-01;
406 //fNonLinearityParams[2] = -1.84235e+02;
407 //fNonLinearityParams[3] = -2.05019e+00;
408 //fNonLinearityParams[4] = -5.89423e+00;
409 energy *= fNonLinearityParams[0]*( ( fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]) ) +
410 ( (fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())) *
411 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3])) ) );
417 AliDebug(2,"No correction on the energy\n");
426 //__________________________________________________
427 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
429 //Calculate shower depth for a given cluster energy and particle type
439 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
443 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
450 gGeoManager->cd("ALIC_1/XEN1_1");
451 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
452 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
454 TGeoVolume *geoSMVol = geoSM->GetVolume();
455 TGeoShape *geoSMShape = geoSMVol->GetShape();
456 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
457 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
458 else AliFatal("Null GEANT box");
459 }else AliFatal("NULL GEANT node matrix");
462 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
468 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
475 //__________________________________________________
476 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
477 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
479 //For a given CaloCluster gets the absId of the cell
480 //with maximum energy deposit.
483 Double_t eCell = -1.;
484 Float_t fraction = 1.;
485 Float_t recalFactor = 1.;
486 Int_t cellAbsId = -1 ;
492 //printf("---Max?\n");
493 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
494 cellAbsId = clu->GetCellAbsId(iDig);
495 fraction = clu->GetCellAmplitudeFraction(iDig);
496 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
497 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
498 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
499 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
500 if(iDig==0) iSupMod0=iSupMod;
501 else if(iSupMod0!=iSupMod) {
503 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
505 if(IsRecalibrationOn()) {
506 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
508 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
509 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
513 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
517 //Get from the absid the supermodule, tower and eta/phi numbers
518 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
519 //Gives SuperModule and Tower numbers
520 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
521 iIphi, iIeta,iphi,ieta);
522 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
523 //printf("Max end---\n");
527 //________________________________________________________________
528 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
529 //Init EMCAL recalibration factors
530 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
531 //In order to avoid rewriting the same histograms
532 Bool_t oldStatus = TH1::AddDirectoryStatus();
533 TH1::AddDirectory(kFALSE);
535 fEMCALRecalibrationFactors = new TObjArray(10);
536 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
537 //Init the histograms with 1
538 for (Int_t sm = 0; sm < 12; sm++) {
539 for (Int_t i = 0; i < 48; i++) {
540 for (Int_t j = 0; j < 24; j++) {
541 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
545 fEMCALRecalibrationFactors->SetOwner(kTRUE);
546 fEMCALRecalibrationFactors->Compress();
548 //In order to avoid rewriting the same histograms
549 TH1::AddDirectory(oldStatus);
553 //________________________________________________________________
554 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
555 //Init EMCAL bad channels map
556 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
557 //In order to avoid rewriting the same histograms
558 Bool_t oldStatus = TH1::AddDirectoryStatus();
559 TH1::AddDirectory(kFALSE);
561 fEMCALBadChannelMap = new TObjArray(10);
562 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
563 for (int i = 0; i < 10; i++) {
564 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
569 fEMCALBadChannelMap->SetOwner(kTRUE);
570 fEMCALBadChannelMap->Compress();
572 //In order to avoid rewriting the same histograms
573 TH1::AddDirectory(oldStatus);
576 //________________________________________________________________
577 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
578 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
580 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
581 UShort_t * index = cluster->GetCellsAbsId() ;
582 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
583 Int_t ncells = cluster->GetNCells();
585 //Initialize some used variables
588 Int_t icol = -1, irow = -1, imod=1;
589 Float_t factor = 1, frac = 0;
591 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
592 for(Int_t icell = 0; icell < ncells; icell++){
593 absId = index[icell];
594 frac = fraction[icell];
595 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
596 Int_t iTower = -1, iIphi = -1, iIeta = -1;
597 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
598 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
599 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
600 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
601 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
602 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
604 energy += cells->GetCellAmplitude(absId)*factor*frac;
608 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
610 cluster->SetE(energy);
615 //__________________________________________________
616 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
618 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
620 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
621 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
622 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
626 //__________________________________________________
627 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
629 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
630 // The algorithm is a copy of what is done in AliEMCALRecPoint
633 Float_t fraction = 1.;
634 Float_t recalFactor = 1.;
637 Int_t iTower = -1, iIphi = -1, iIeta = -1;
638 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
639 Float_t weight = 0., totalWeight=0.;
640 Float_t newPos[3] = {0,0,0};
641 Double_t pLocal[3], pGlobal[3];
642 Bool_t shared = kFALSE;
644 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
645 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
646 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
648 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
650 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
651 absId = clu->GetCellAbsId(iDig);
652 fraction = clu->GetCellAmplitudeFraction(iDig);
653 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
654 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
655 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
657 if(IsRecalibrationOn()) {
658 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
660 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
662 weight = GetCellWeight(eCell,clEnergy);
663 //printf("cell energy %f, weight %f\n",eCell,weight);
664 totalWeight += weight;
665 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
666 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
667 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
668 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
670 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
675 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
678 //Float_t pos[]={0,0,0};
679 //clu->GetPosition(pos);
680 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
681 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
683 if(iSupModMax > 1) {//sector 1
684 newPos[0] +=fMisalTransShift[3];//-=3.093;
685 newPos[1] +=fMisalTransShift[4];//+=6.82;
686 newPos[2] +=fMisalTransShift[5];//+=1.635;
687 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
691 newPos[0] +=fMisalTransShift[0];//+=1.134;
692 newPos[1] +=fMisalTransShift[1];//+=8.2;
693 newPos[2] +=fMisalTransShift[2];//+=1.197;
694 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
697 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
699 clu->SetPosition(newPos);
703 //__________________________________________________
704 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
706 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
707 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
710 Float_t fraction = 1.;
711 Float_t recalFactor = 1.;
715 Int_t iIphi = -1, iIeta = -1;
716 Int_t iSupMod = -1, iSupModMax = -1;
717 Int_t iphi = -1, ieta =-1;
718 Bool_t shared = kFALSE;
720 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
721 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
722 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
724 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
725 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
726 Int_t startingSM = -1;
728 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
729 absId = clu->GetCellAbsId(iDig);
730 fraction = clu->GetCellAmplitudeFraction(iDig);
731 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
732 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
733 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
735 if (iDig==0) startingSM = iSupMod;
736 else if(iSupMod != startingSM) areInSameSM = kFALSE;
738 eCell = cells->GetCellAmplitude(absId);
740 if(IsRecalibrationOn()) {
741 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
743 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
745 weight = GetCellWeight(eCell,clEnergy);
746 if(weight < 0) weight = 0;
747 totalWeight += weight;
748 weightedCol += ieta*weight;
749 weightedRow += iphi*weight;
751 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
755 Float_t xyzNew[]={0.,0.,0.};
756 if(areInSameSM == kTRUE) {
757 //printf("In Same SM\n");
758 weightedCol = weightedCol/totalWeight;
759 weightedRow = weightedRow/totalWeight;
760 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
763 //printf("In Different SM\n");
764 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
767 clu->SetPosition(xyzNew);
771 //____________________________________________________________________________
772 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
774 //re-evaluate distance to bad channel with updated bad map
776 if(!fRecalDistToBadChannels) return;
778 //Get channels map of the supermodule where the cluster is.
779 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
780 Bool_t shared = kFALSE;
781 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
782 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
785 Float_t minDist = 10000.;
788 //Loop on tower status map
789 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
790 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
791 //Check if tower is bad.
792 if(hMap->GetBinContent(icol,irow)==0) continue;
793 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
794 // iSupMod,icol, irow, icolM,irowM);
796 dRrow=TMath::Abs(irowM-irow);
797 dRcol=TMath::Abs(icolM-icol);
798 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
800 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
807 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
812 //The only possible combinations are (0,1), (2,3) ... (8,9)
813 if(iSupMod%2) iSupMod2 = iSupMod-1;
814 else iSupMod2 = iSupMod+1;
815 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
817 //Loop on tower status map of second super module
818 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
819 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
820 //Check if tower is bad.
821 if(hMap2->GetBinContent(icol,irow)==0) continue;
822 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
823 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
825 dRrow=TMath::Abs(irow-irowM);
828 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
831 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
834 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
835 if(dist < minDist) minDist = dist;
840 }// shared cluster in 2 SuperModules
842 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
843 cluster->SetDistanceToBadChannel(minDist);
847 //____________________________________________________________________________
848 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
850 //re-evaluate identification parameters with bayesian
852 if ( cluster->GetM02() != 0)
853 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
855 Float_t pidlist[AliPID::kSPECIESN+1];
856 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
858 cluster->SetPID(pidlist);
862 //____________________________________________________________________________
863 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
865 // Calculates new center of gravity in the local EMCAL-module coordinates
866 // and tranfers into global ALICE coordinates
867 // Calculates Dispersion and main axis
872 Float_t fraction = 1.;
873 Float_t recalFactor = 1.;
893 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
895 //Get from the absid the supermodule, tower and eta/phi numbers
896 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
897 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
899 //Get the cell energy, if recalibration is on, apply factors
900 fraction = cluster->GetCellAmplitudeFraction(iDigit);
901 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
902 if(IsRecalibrationOn()) {
903 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
905 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
907 if(cluster->E() > 0 && eCell > 0){
909 w = GetCellWeight(eCell,cluster->E());
917 dxx += w * etai * etai ;
919 dzz += w * phii * phii ;
921 dxz += w * etai * phii ;
925 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
928 //Normalize to the weight
934 AliError(Form("Wrong weight %f\n", wtot));
936 //Calculate dispersion
937 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
939 //Get from the absid the supermodule, tower and eta/phi numbers
940 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
941 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
943 //Get the cell energy, if recalibration is on, apply factors
944 fraction = cluster->GetCellAmplitudeFraction(iDigit);
945 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
946 if(IsRecalibrationOn()) {
947 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
949 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
951 if(cluster->E() > 0 && eCell > 0){
953 w = GetCellWeight(eCell,cluster->E());
957 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
960 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
963 //Normalize to the weigth and set shower shape parameters
964 if (wtot > 0 && nstat > 1) {
969 dxx -= xmean * xmean ;
970 dzz -= zmean * zmean ;
971 dxz -= xmean * zmean ;
972 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
973 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
977 cluster->SetM20(0.) ;
978 cluster->SetM02(0.) ;
982 cluster->SetDispersion(TMath::Sqrt(d)) ;
984 cluster->SetDispersion(0) ;
988 //____________________________________________________________________________
989 void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr)
991 //This function should be called before the cluster loop
992 //Before call this function, please recalculate the cluster positions
993 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
994 //Store matched cluster indexes and residuals
995 //It only works with ESDs, not AODs
997 fMatchedTrackIndex ->Reset();
998 fMatchedClusterIndex->Reset();
999 fResidualZ ->Reset();
1000 fResidualR ->Reset();
1002 fMatchedTrackIndex ->Set(500);
1003 fMatchedClusterIndex->Set(500);
1004 fResidualZ ->Set(500);
1005 fResidualR ->Set(500);
1010 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1012 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
1013 if(!track || !IsAccepted(track)) continue;
1015 Float_t dRMax = fCutR, dZMax = fCutZ;
1017 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
1018 if(!clusterArr){// get clusters from event
1019 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1021 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1022 if(!cluster->IsEMCAL()) continue;
1023 cluster->GetPosition(clsPos); //Has been recalculated
1024 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1025 emctrack->GetXYZ(trkPos);
1026 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1027 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1036 } else { // external cluster array, not from ESD event
1037 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1039 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1040 if(!cluster->IsEMCAL()) continue;
1041 cluster->GetPosition(clsPos); //Has been recalculated
1042 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1043 emctrack->GetXYZ(trkPos);
1044 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1045 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1054 }// external list of clusters
1058 fMatchedTrackIndex ->AddAt(itr,matched);
1059 fMatchedClusterIndex->AddAt(index,matched);
1060 fResidualZ ->AddAt(dZMax,matched);
1061 fResidualR ->AddAt(dRMax,matched);
1067 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1069 fMatchedTrackIndex ->Set(matched);
1070 fMatchedClusterIndex->Set(matched);
1071 fResidualZ ->Set(matched);
1072 fResidualR ->Set(matched);
1074 //printf("Number of matched pairs: %d\n",matched);
1077 //________________________________________________________________________________
1078 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1080 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1081 //Get the residuals dR and dZ for this cluster
1082 //It only works with ESDs, not AODs
1084 if( FindMatchedPos(index) >= 999 )
1086 AliDebug(2,"No matched tracks found!\n");
1091 dR = fResidualR->At(FindMatchedPos(index));
1092 dZ = fResidualZ->At(FindMatchedPos(index));
1093 //printf("dR %f, dZ %f\n",dR, dZ);
1096 //__________________________________________________________
1097 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t index)
1099 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1100 //Get the index of matched track for this cluster
1101 //It only works with ESDs, not AODs
1103 if(IsMatched(index))
1104 return fMatchedTrackIndex->At(FindMatchedPos(index));
1110 //__________________________________________________
1111 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1113 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1114 //Returns if cluster has a match
1115 if(FindMatchedPos(index) < 999)
1120 //__________________________________________________________
1121 UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
1123 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1124 //Returns the position of the match in the fMatchedClusterIndex array
1125 Float_t tmpR = fCutR;
1128 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1130 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1133 tmpR=fResidualR->At(i);
1134 AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1140 //__________________________________________________________
1141 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1143 // Given a esd track, return whether the track survive all the cuts
1145 // The different quality parameter are first
1146 // retrieved from the track. then it is found out what cuts the
1147 // track did not survive and finally the cuts are imposed.
1149 UInt_t status = esdTrack->GetStatus();
1151 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1152 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1154 Float_t chi2PerClusterITS = -1;
1155 Float_t chi2PerClusterTPC = -1;
1156 if (nClustersITS!=0)
1157 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1158 if (nClustersTPC!=0)
1159 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1163 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1164 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1165 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1170 esdTrack->GetImpactParameters(b,bCov);
1171 if (bCov[0]<=0 || bCov[2]<=0) {
1172 AliDebug(1, "Estimated b resolution lower or equal zero!");
1173 bCov[0]=0; bCov[2]=0;
1176 Float_t dcaToVertexXY = b[0];
1177 Float_t dcaToVertexZ = b[1];
1178 Float_t dcaToVertex = -1;
1180 if (fCutDCAToVertex2D)
1181 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1183 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1187 Bool_t cuts[kNCuts];
1188 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1190 // track quality cuts
1191 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1193 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1195 if (nClustersTPC<fCutMinNClusterTPC)
1197 if (nClustersITS<fCutMinNClusterITS)
1199 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1201 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1203 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1205 if (fCutDCAToVertex2D && dcaToVertex > 1)
1207 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1209 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1212 //Require at least one SPD point + anything else in ITS
1213 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1217 for (Int_t i=0; i<kNCuts; i++)
1218 if (cuts[i]) {cut = kTRUE;}
1226 //__________________________________________________
1227 void AliEMCALRecoUtils::InitTrackCuts()
1229 //Intilize the track cut criteria
1230 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1231 //Also you can customize the cuts using the setters
1234 SetMinNClustersTPC(70);
1235 SetMaxChi2PerClusterTPC(4);
1236 SetAcceptKinkDaughters(kFALSE);
1237 SetRequireTPCRefit(kTRUE);
1240 SetRequireITSRefit(kTRUE);
1241 SetMaxDCAToVertexZ(2);
1242 SetDCAToVertex2D(kFALSE);
1243 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1244 SetMinNClustersITS();
1247 //___________________________________________________
1248 void AliEMCALRecoUtils::Print(const Option_t *) const
1252 printf("AliEMCALRecoUtils Settings: \n");
1253 printf("Misalignment shifts\n");
1254 for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
1255 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1256 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1257 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1258 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1260 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1262 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1264 printf("Track cuts: \n");
1265 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1266 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1267 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1268 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1269 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1273 //_____________________________________________________________________
1274 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1275 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1276 //Do it only once and only if it is requested
1278 if(!fUseTimeCorrectionFactors) return;
1279 if(fTimeCorrectionFactorsSet) return;
1281 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1283 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1284 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1286 SwitchOnRecalibration();
1287 for(Int_t ism = 0; ism < 4; ism++){
1288 for(Int_t icol = 0; icol < 48; icol++){
1289 for(Int_t irow = 0; irow < 24; irow++){
1290 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1291 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1292 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1293 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1294 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1295 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1299 fTimeCorrectionFactorsSet = kTRUE;