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]));
403 AliDebug(2,"No correction on the energy\n");
412 //__________________________________________________
413 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
415 //Calculate shower depth for a given cluster energy and particle type
425 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
429 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
436 gGeoManager->cd("ALIC_1/XEN1_1");
437 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
438 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
440 TGeoVolume *geoSMVol = geoSM->GetVolume();
441 TGeoShape *geoSMShape = geoSMVol->GetShape();
442 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
443 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
444 else AliFatal("Null GEANT box");
445 }else AliFatal("NULL GEANT node matrix");
448 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
454 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
461 //__________________________________________________
462 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
463 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
465 //For a given CaloCluster gets the absId of the cell
466 //with maximum energy deposit.
469 Double_t eCell = -1.;
470 Float_t fraction = 1.;
471 Float_t recalFactor = 1.;
472 Int_t cellAbsId = -1 ;
478 //printf("---Max?\n");
479 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
480 cellAbsId = clu->GetCellAbsId(iDig);
481 fraction = clu->GetCellAmplitudeFraction(iDig);
482 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
483 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
484 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
485 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
486 if(iDig==0) iSupMod0=iSupMod;
487 else if(iSupMod0!=iSupMod) {
489 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
491 if(IsRecalibrationOn()) {
492 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
494 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
495 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
499 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
503 //Get from the absid the supermodule, tower and eta/phi numbers
504 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
505 //Gives SuperModule and Tower numbers
506 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
507 iIphi, iIeta,iphi,ieta);
508 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
509 //printf("Max end---\n");
513 //________________________________________________________________
514 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
515 //Init EMCAL recalibration factors
516 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
517 //In order to avoid rewriting the same histograms
518 Bool_t oldStatus = TH1::AddDirectoryStatus();
519 TH1::AddDirectory(kFALSE);
521 fEMCALRecalibrationFactors = new TObjArray(10);
522 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
523 //Init the histograms with 1
524 for (Int_t sm = 0; sm < 12; sm++) {
525 for (Int_t i = 0; i < 48; i++) {
526 for (Int_t j = 0; j < 24; j++) {
527 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
531 fEMCALRecalibrationFactors->SetOwner(kTRUE);
532 fEMCALRecalibrationFactors->Compress();
534 //In order to avoid rewriting the same histograms
535 TH1::AddDirectory(oldStatus);
539 //________________________________________________________________
540 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
541 //Init EMCAL bad channels map
542 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
543 //In order to avoid rewriting the same histograms
544 Bool_t oldStatus = TH1::AddDirectoryStatus();
545 TH1::AddDirectory(kFALSE);
547 fEMCALBadChannelMap = new TObjArray(10);
548 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
549 for (int i = 0; i < 10; i++) {
550 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
555 fEMCALBadChannelMap->SetOwner(kTRUE);
556 fEMCALBadChannelMap->Compress();
558 //In order to avoid rewriting the same histograms
559 TH1::AddDirectory(oldStatus);
562 //________________________________________________________________
563 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
564 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
566 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
567 UShort_t * index = cluster->GetCellsAbsId() ;
568 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
569 Int_t ncells = cluster->GetNCells();
571 //Initialize some used variables
574 Int_t icol = -1, irow = -1, imod=1;
575 Float_t factor = 1, frac = 0;
577 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
578 for(Int_t icell = 0; icell < ncells; icell++){
579 absId = index[icell];
580 frac = fraction[icell];
581 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
582 Int_t iTower = -1, iIphi = -1, iIeta = -1;
583 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
584 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
585 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
586 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
587 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
588 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
590 energy += cells->GetCellAmplitude(absId)*factor*frac;
594 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
596 cluster->SetE(energy);
601 //__________________________________________________
602 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
604 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
606 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
607 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
608 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
612 //__________________________________________________
613 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
615 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
616 // The algorithm is a copy of what is done in AliEMCALRecPoint
619 Float_t fraction = 1.;
620 Float_t recalFactor = 1.;
623 Int_t iTower = -1, iIphi = -1, iIeta = -1;
624 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
625 Float_t weight = 0., totalWeight=0.;
626 Float_t newPos[3] = {0,0,0};
627 Double_t pLocal[3], pGlobal[3];
628 Bool_t shared = kFALSE;
630 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
631 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
632 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
634 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
636 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
637 absId = clu->GetCellAbsId(iDig);
638 fraction = clu->GetCellAmplitudeFraction(iDig);
639 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
640 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
641 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
643 if(IsRecalibrationOn()) {
644 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
646 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
648 weight = GetCellWeight(eCell,clEnergy);
649 //printf("cell energy %f, weight %f\n",eCell,weight);
650 totalWeight += weight;
651 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
652 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
653 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
654 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
656 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
661 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
664 //Float_t pos[]={0,0,0};
665 //clu->GetPosition(pos);
666 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
667 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
669 if(iSupModMax > 1) {//sector 1
670 newPos[0] +=fMisalTransShift[3];//-=3.093;
671 newPos[1] +=fMisalTransShift[4];//+=6.82;
672 newPos[2] +=fMisalTransShift[5];//+=1.635;
673 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
677 newPos[0] +=fMisalTransShift[0];//+=1.134;
678 newPos[1] +=fMisalTransShift[1];//+=8.2;
679 newPos[2] +=fMisalTransShift[2];//+=1.197;
680 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
683 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
685 clu->SetPosition(newPos);
689 //__________________________________________________
690 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
692 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
693 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
696 Float_t fraction = 1.;
697 Float_t recalFactor = 1.;
701 Int_t iIphi = -1, iIeta = -1;
702 Int_t iSupMod = -1, iSupModMax = -1;
703 Int_t iphi = -1, ieta =-1;
704 Bool_t shared = kFALSE;
706 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
707 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
708 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
710 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
711 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
712 Int_t startingSM = -1;
714 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
715 absId = clu->GetCellAbsId(iDig);
716 fraction = clu->GetCellAmplitudeFraction(iDig);
717 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
718 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
719 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
721 if (iDig==0) startingSM = iSupMod;
722 else if(iSupMod != startingSM) areInSameSM = kFALSE;
724 eCell = cells->GetCellAmplitude(absId);
726 if(IsRecalibrationOn()) {
727 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
729 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
731 weight = GetCellWeight(eCell,clEnergy);
732 if(weight < 0) weight = 0;
733 totalWeight += weight;
734 weightedCol += ieta*weight;
735 weightedRow += iphi*weight;
737 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
741 Float_t xyzNew[]={0.,0.,0.};
742 if(areInSameSM == kTRUE) {
743 //printf("In Same SM\n");
744 weightedCol = weightedCol/totalWeight;
745 weightedRow = weightedRow/totalWeight;
746 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
749 //printf("In Different SM\n");
750 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
753 clu->SetPosition(xyzNew);
757 //____________________________________________________________________________
758 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
760 //re-evaluate distance to bad channel with updated bad map
762 if(!fRecalDistToBadChannels) return;
764 //Get channels map of the supermodule where the cluster is.
765 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
766 Bool_t shared = kFALSE;
767 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
768 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
771 Float_t minDist = 10000.;
774 //Loop on tower status map
775 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
776 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
777 //Check if tower is bad.
778 if(hMap->GetBinContent(icol,irow)==0) continue;
779 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
780 // iSupMod,icol, irow, icolM,irowM);
782 dRrow=TMath::Abs(irowM-irow);
783 dRcol=TMath::Abs(icolM-icol);
784 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
786 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
793 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
798 //The only possible combinations are (0,1), (2,3) ... (8,9)
799 if(iSupMod%2) iSupMod2 = iSupMod-1;
800 else iSupMod2 = iSupMod+1;
801 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
803 //Loop on tower status map of second super module
804 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
805 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
806 //Check if tower is bad.
807 if(hMap2->GetBinContent(icol,irow)==0) continue;
808 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
809 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
811 dRrow=TMath::Abs(irow-irowM);
814 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
817 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
820 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
821 if(dist < minDist) minDist = dist;
826 }// shared cluster in 2 SuperModules
828 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
829 cluster->SetDistanceToBadChannel(minDist);
833 //____________________________________________________________________________
834 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
836 //re-evaluate identification parameters with bayesian
838 if ( cluster->GetM02() != 0)
839 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
841 Float_t pidlist[AliPID::kSPECIESN+1];
842 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
844 cluster->SetPID(pidlist);
848 //____________________________________________________________________________
849 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
851 // Calculates new center of gravity in the local EMCAL-module coordinates
852 // and tranfers into global ALICE coordinates
853 // Calculates Dispersion and main axis
858 Float_t fraction = 1.;
859 Float_t recalFactor = 1.;
879 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
881 //Get from the absid the supermodule, tower and eta/phi numbers
882 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
883 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
885 //Get the cell energy, if recalibration is on, apply factors
886 fraction = cluster->GetCellAmplitudeFraction(iDigit);
887 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
888 if(IsRecalibrationOn()) {
889 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
891 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
893 if(cluster->E() > 0 && eCell > 0){
895 w = GetCellWeight(eCell,cluster->E());
903 dxx += w * etai * etai ;
905 dzz += w * phii * phii ;
907 dxz += w * etai * phii ;
911 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
914 //Normalize to the weight
920 AliError(Form("Wrong weight %f\n", wtot));
922 //Calculate dispersion
923 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
925 //Get from the absid the supermodule, tower and eta/phi numbers
926 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
927 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
929 //Get the cell energy, if recalibration is on, apply factors
930 fraction = cluster->GetCellAmplitudeFraction(iDigit);
931 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
932 if(IsRecalibrationOn()) {
933 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
935 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
937 if(cluster->E() > 0 && eCell > 0){
939 w = GetCellWeight(eCell,cluster->E());
943 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
946 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
949 //Normalize to the weigth and set shower shape parameters
950 if (wtot > 0 && nstat > 1) {
955 dxx -= xmean * xmean ;
956 dzz -= zmean * zmean ;
957 dxz -= xmean * zmean ;
958 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
959 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
963 cluster->SetM20(0.) ;
964 cluster->SetM02(0.) ;
968 cluster->SetDispersion(TMath::Sqrt(d)) ;
970 cluster->SetDispersion(0) ;
974 //____________________________________________________________________________
975 void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr)
977 //This function should be called before the cluster loop
978 //Before call this function, please recalculate the cluster positions
979 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
980 //Store matched cluster indexes and residuals
981 //It only works with ESDs, not AODs
983 fMatchedTrackIndex ->Reset();
984 fMatchedClusterIndex->Reset();
985 fResidualZ ->Reset();
986 fResidualR ->Reset();
988 fMatchedTrackIndex ->Set(100);
989 fMatchedClusterIndex->Set(100);
990 fResidualZ ->Set(100);
991 fResidualR ->Set(100);
996 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
998 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
999 if(!track || !IsAccepted(track)) continue;
1001 Float_t dRMax = fCutR, dZMax = fCutZ;
1003 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
1004 if(!clusterArr){// get clusters from event
1005 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1007 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1008 if(!cluster->IsEMCAL()) continue;
1009 cluster->GetPosition(clsPos); //Has been recalculated
1010 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1011 emctrack->GetXYZ(trkPos);
1012 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1013 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1022 } else { // external cluster array, not from ESD event
1023 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1025 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1026 if(!cluster->IsEMCAL()) continue;
1027 cluster->GetPosition(clsPos); //Has been recalculated
1028 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
1029 emctrack->GetXYZ(trkPos);
1030 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1031 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1040 }// external list of clusters
1044 fMatchedTrackIndex ->AddAt(itr,matched);
1045 fMatchedClusterIndex->AddAt(index,matched);
1046 fResidualZ ->AddAt(dZMax,matched);
1047 fResidualR ->AddAt(dRMax,matched);
1053 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1055 fMatchedTrackIndex ->Set(matched);
1056 fMatchedClusterIndex->Set(matched);
1057 fResidualZ ->Set(matched);
1058 fResidualR ->Set(matched);
1060 //printf("Number of matched pairs: %d\n",matched);
1063 //________________________________________________________________________________
1064 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1066 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1067 //Get the residuals dR and dZ for this cluster
1068 //It only works with ESDs, not AODs
1070 if( FindMatchedPos(index) >= 999 )
1072 AliDebug(2,"No matched tracks found!\n");
1077 dR = fResidualR->At(FindMatchedPos(index));
1078 dZ = fResidualZ->At(FindMatchedPos(index));
1079 //printf("dR %f, dZ %f\n",dR, dZ);
1082 //__________________________________________________________
1083 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t index)
1085 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1086 //Get the index of matched track for this cluster
1087 //It only works with ESDs, not AODs
1089 if(IsMatched(index))
1090 return fMatchedTrackIndex->At(FindMatchedPos(index));
1096 //__________________________________________________
1097 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1099 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1100 //Returns if cluster has a match
1101 if(FindMatchedPos(index) < 999)
1106 //__________________________________________________________
1107 UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
1109 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1110 //Returns the position of the match in the fMatchedClusterIndex array
1111 Float_t tmpR = fCutR;
1114 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1116 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1119 tmpR=fResidualR->At(i);
1120 AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1126 //__________________________________________________________
1127 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1129 // Given a esd track, return whether the track survive all the cuts
1131 // The different quality parameter are first
1132 // retrieved from the track. then it is found out what cuts the
1133 // track did not survive and finally the cuts are imposed.
1135 UInt_t status = esdTrack->GetStatus();
1137 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1138 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1140 Float_t chi2PerClusterITS = -1;
1141 Float_t chi2PerClusterTPC = -1;
1142 if (nClustersITS!=0)
1143 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1144 if (nClustersTPC!=0)
1145 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1149 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1150 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1151 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1156 esdTrack->GetImpactParameters(b,bCov);
1157 if (bCov[0]<=0 || bCov[2]<=0) {
1158 AliDebug(1, "Estimated b resolution lower or equal zero!");
1159 bCov[0]=0; bCov[2]=0;
1162 Float_t dcaToVertexXY = b[0];
1163 Float_t dcaToVertexZ = b[1];
1164 Float_t dcaToVertex = -1;
1166 if (fCutDCAToVertex2D)
1167 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1169 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1173 Bool_t cuts[kNCuts];
1174 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1176 // track quality cuts
1177 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1179 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1181 if (nClustersTPC<fCutMinNClusterTPC)
1183 if (nClustersITS<fCutMinNClusterITS)
1185 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1187 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1189 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1191 if (fCutDCAToVertex2D && dcaToVertex > 1)
1193 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1195 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1198 //Require at least one SPD point + anything else in ITS
1199 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1203 for (Int_t i=0; i<kNCuts; i++)
1204 if (cuts[i]) {cut = kTRUE;}
1212 //__________________________________________________
1213 void AliEMCALRecoUtils::InitTrackCuts()
1215 //Intilize the track cut criteria
1216 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1217 //Also you can customize the cuts using the setters
1220 SetMinNClustersTPC(70);
1221 SetMaxChi2PerClusterTPC(4);
1222 SetAcceptKinkDaughters(kFALSE);
1223 SetRequireTPCRefit(kTRUE);
1226 SetRequireITSRefit(kTRUE);
1227 SetMaxDCAToVertexZ(2);
1228 SetDCAToVertex2D(kFALSE);
1229 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1230 SetMinNClustersITS();
1233 //___________________________________________________
1234 void AliEMCALRecoUtils::Print(const Option_t *) const
1238 printf("AliEMCALRecoUtils Settings: \n");
1239 printf("Misalignment shifts\n");
1240 for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
1241 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1242 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1243 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1244 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1246 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1248 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1250 printf("Track cuts: \n");
1251 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1252 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1253 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1254 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1255 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1259 //_____________________________________________________________________
1260 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1261 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1262 //Do it only once and only if it is requested
1264 if(!fUseTimeCorrectionFactors) return;
1265 if(fTimeCorrectionFactorsSet) return;
1267 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1269 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1270 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1272 SwitchOnRecalibration();
1273 for(Int_t ism = 0; ism < 4; ism++){
1274 for(Int_t icol = 0; icol < 48; icol++){
1275 for(Int_t irow = 0; irow < 24; irow++){
1276 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1277 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1278 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1279 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1280 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1281 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1285 fTimeCorrectionFactorsSet = kTRUE;