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 ///////////////////////////////////////////////////////////////////////////////
29 // standard C++ includes
30 //#include <Riostream.h>
33 #include <TGeoManager.h>
34 #include <TGeoMatrix.h>
38 #include "AliEMCALRecoUtils.h"
39 #include "AliEMCALGeometry.h"
40 #include "AliVCluster.h"
41 #include "AliVCaloCells.h"
42 #include "AliVEvent.h"
44 #include "AliEMCALPIDUtils.h"
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliEMCALTrack.h"
50 ClassImp(AliEMCALRecoUtils)
52 //______________________________________________
53 AliEMCALRecoUtils::AliEMCALRecoUtils():
54 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
55 fPosAlgo(kUnchanged),fW0(4.),
56 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
57 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
58 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
59 fMatchedClusterIndex(0x0), fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
60 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
61 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
62 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),
67 // Initialize all constant values which have to be used
68 // during Reco algorithm execution
71 for(Int_t i = 0; i < 15 ; i++) {
72 fMisalTransShift[i] = 0.;
73 fMisalRotShift[i] = 0.;
75 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
76 //For kPi0GammaGamma case, but default is no correction
77 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
78 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
79 fNonLinearityParams[2] = 1.046;
81 fMatchedClusterIndex = new TArrayI();
82 fResidualZ = new TArrayF();
83 fResidualR = new TArrayF();
85 fPIDUtils = new AliEMCALPIDUtils();
91 //______________________________________________________________________
92 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
93 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
94 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
95 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
96 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
97 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
98 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
99 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
100 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
101 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
102 fCutR(reco.fCutR),fCutZ(reco.fCutZ),
103 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
104 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
105 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
106 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
107 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
108 fPIDUtils(reco.fPIDUtils)
113 for(Int_t i = 0; i < 15 ; i++) {
114 fMisalRotShift[i] = reco.fMisalRotShift[i];
115 fMisalTransShift[i] = reco.fMisalTransShift[i];
117 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
122 //______________________________________________________________________
123 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
125 //Assignment operator
127 if(this == &reco)return *this;
128 ((TNamed *)this)->operator=(reco);
130 fNonLinearityFunction = reco.fNonLinearityFunction;
131 fParticleType = reco.fParticleType;
132 fPosAlgo = reco.fPosAlgo;
134 fRecalibration = reco.fRecalibration;
135 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
136 fRemoveBadChannels = reco.fRemoveBadChannels;
137 fRecalDistToBadChannels= reco.fRecalDistToBadChannels;
138 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
139 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
140 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
143 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
144 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
149 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
150 fCutMinNClusterITS = reco.fCutMinNClusterITS;
151 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
152 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
153 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
154 fCutRequireITSRefit = reco.fCutRequireITSRefit;
155 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
156 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
157 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
158 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
160 fPIDUtils = reco.fPIDUtils;
164 // assign or copy construct
166 *fResidualR = *reco.fResidualR;
168 else fResidualR = new TArrayF(*reco.fResidualR);
171 if(fResidualR)delete fResidualR;
176 // assign or copy construct
178 *fResidualZ = *reco.fResidualZ;
180 else fResidualZ = new TArrayF(*reco.fResidualZ);
183 if(fResidualZ)delete fResidualZ;
188 if(reco.fMatchedClusterIndex){
189 // assign or copy construct
190 if(fMatchedClusterIndex){
191 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
193 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
196 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
197 fMatchedClusterIndex = 0;
205 //__________________________________________________
206 AliEMCALRecoUtils::~AliEMCALRecoUtils()
210 if(fEMCALRecalibrationFactors) {
211 fEMCALRecalibrationFactors->Clear();
212 delete fEMCALRecalibrationFactors;
215 if(fEMCALBadChannelMap) {
216 fEMCALBadChannelMap->Clear();
217 delete fEMCALBadChannelMap;
220 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
221 if(fResidualR) {delete fResidualR; fResidualR=0;}
222 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
226 //_______________________________________________________________
227 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
229 // Given the list of AbsId of the cluster, get the maximum cell and
230 // check if there are fNCellsFromBorder from the calorimeter border
232 //If the distance to the border is 0 or negative just exit accept all clusters
233 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
235 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
236 Bool_t shared = kFALSE;
237 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
239 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
240 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
242 if(absIdMax==-1) return kFALSE;
244 //Check if the cell is close to the borders:
245 Bool_t okrow = kFALSE;
246 Bool_t okcol = kFALSE;
248 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
249 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
255 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
258 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
262 if(!fNoEMCALBorderAtEta0){
263 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
267 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
270 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
274 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
275 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
277 if (okcol && okrow) {
278 //printf("Accept\n");
282 //printf("Reject\n");
283 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
290 //_________________________________________________________________________________________________________
291 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
292 // Check that in the cluster cells, there is no bad channel of those stored
293 // in fEMCALBadChannelMap or fPHOSBadChannelMap
295 if(!fRemoveBadChannels) return kFALSE;
296 if(!fEMCALBadChannelMap) return kFALSE;
301 for(Int_t iCell = 0; iCell<nCells; iCell++){
303 //Get the column and row
304 Int_t iTower = -1, iIphi = -1, iIeta = -1;
305 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
306 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
307 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
308 if(GetEMCALChannelStatus(imod, icol, irow)){
309 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
313 }// cell cluster loop
319 //__________________________________________________
320 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
321 // Correct cluster energy from non linearity functions
322 Float_t energy = cluster->E();
324 switch (fNonLinearityFunction) {
328 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
329 //Double_t fNonLinearityParams[0] = 1.001;
330 //Double_t fNonLinearityParams[1] = -0.01264;
331 //Double_t fNonLinearityParams[2] = -0.03632;
332 //Double_t fNonLinearityParams[3] = 0.1798;
333 //Double_t fNonLinearityParams[4] = -0.522;
334 energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
335 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
336 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
342 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
343 //Double_t fNonLinearityParams[0] = 0.1457;
344 //Double_t fNonLinearityParams[1] = -0.02024;
345 //Double_t fNonLinearityParams[2] = 1.046;
346 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
350 case kPi0GammaConversion:
352 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
353 //fNonLinearityParams[0] = 0.139393/0.1349766;
354 //fNonLinearityParams[1] = 0.0566186;
355 //fNonLinearityParams[2] = 0.982133;
356 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
363 //From beam test, Alexei's results, for different ZS thresholds
364 // th=30 MeV; th = 45 MeV; th = 75 MeV
365 //fNonLinearityParams[0] = 0.107; 1.003; 1.002
366 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
367 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
368 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
374 AliDebug(2,"No correction on the energy\n");
383 //__________________________________________________
384 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
386 //Calculate shower depth for a given cluster energy and particle type
396 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
400 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
407 gGeoManager->cd("ALIC_1/XEN1_1");
408 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
409 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
411 TGeoVolume *geoSMVol = geoSM->GetVolume();
412 TGeoShape *geoSMShape = geoSMVol->GetShape();
413 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
414 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
415 else AliFatal("Null GEANT box");
416 }else AliFatal("NULL GEANT node matrix");
419 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
425 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
432 //__________________________________________________
433 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
434 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
436 //For a given CaloCluster gets the absId of the cell
437 //with maximum energy deposit.
440 Double_t eCell = -1.;
441 Float_t fraction = 1.;
442 Float_t recalFactor = 1.;
443 Int_t cellAbsId = -1 ;
449 //printf("---Max?\n");
450 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
451 cellAbsId = clu->GetCellAbsId(iDig);
452 fraction = clu->GetCellAmplitudeFraction(iDig);
453 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
454 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
455 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
456 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
457 if(iDig==0) iSupMod0=iSupMod;
458 else if(iSupMod0!=iSupMod) {
460 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
462 if(IsRecalibrationOn()) {
463 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
465 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
466 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
470 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
474 //Get from the absid the supermodule, tower and eta/phi numbers
475 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
476 //Gives SuperModule and Tower numbers
477 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
478 iIphi, iIeta,iphi,ieta);
479 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
480 //printf("Max end---\n");
484 //________________________________________________________________
485 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
486 //Init EMCAL recalibration factors
487 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
488 //In order to avoid rewriting the same histograms
489 Bool_t oldStatus = TH1::AddDirectoryStatus();
490 TH1::AddDirectory(kFALSE);
492 fEMCALRecalibrationFactors = new TObjArray(10);
493 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
494 //Init the histograms with 1
495 for (Int_t sm = 0; sm < 12; sm++) {
496 for (Int_t i = 0; i < 48; i++) {
497 for (Int_t j = 0; j < 24; j++) {
498 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
502 fEMCALRecalibrationFactors->SetOwner(kTRUE);
503 fEMCALRecalibrationFactors->Compress();
505 //In order to avoid rewriting the same histograms
506 TH1::AddDirectory(oldStatus);
510 //________________________________________________________________
511 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
512 //Init EMCAL bad channels map
513 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
514 //In order to avoid rewriting the same histograms
515 Bool_t oldStatus = TH1::AddDirectoryStatus();
516 TH1::AddDirectory(kFALSE);
518 fEMCALBadChannelMap = new TObjArray(10);
519 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
520 for (int i = 0; i < 10; i++) {
521 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
526 fEMCALBadChannelMap->SetOwner(kTRUE);
527 fEMCALBadChannelMap->Compress();
529 //In order to avoid rewriting the same histograms
530 TH1::AddDirectory(oldStatus);
533 //________________________________________________________________
534 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
535 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
537 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
538 UShort_t * index = cluster->GetCellsAbsId() ;
539 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
540 Int_t ncells = cluster->GetNCells();
542 //Initialize some used variables
545 Int_t icol = -1, irow = -1, imod=1;
546 Float_t factor = 1, frac = 0;
548 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
549 for(Int_t icell = 0; icell < ncells; icell++){
550 absId = index[icell];
551 frac = fraction[icell];
552 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
553 Int_t iTower = -1, iIphi = -1, iIeta = -1;
554 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
555 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
556 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
557 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
558 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
559 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
561 energy += cells->GetCellAmplitude(absId)*factor*frac;
565 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
567 cluster->SetE(energy);
572 //__________________________________________________
573 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
575 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
577 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
578 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
579 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
583 //__________________________________________________
584 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
586 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
587 // The algorithm is a copy of what is done in AliEMCALRecPoint
590 Float_t fraction = 1.;
591 Float_t recalFactor = 1.;
594 Int_t iTower = -1, iIphi = -1, iIeta = -1;
595 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
596 Float_t weight = 0., totalWeight=0.;
597 Float_t newPos[3] = {0,0,0};
598 Double_t pLocal[3], pGlobal[3];
599 Bool_t shared = kFALSE;
601 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
602 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
603 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
605 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
607 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
608 absId = clu->GetCellAbsId(iDig);
609 fraction = clu->GetCellAmplitudeFraction(iDig);
610 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
611 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
612 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
614 if(IsRecalibrationOn()) {
615 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
617 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
619 weight = GetCellWeight(eCell,clEnergy);
620 //printf("cell energy %f, weight %f\n",eCell,weight);
621 totalWeight += weight;
622 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
623 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
624 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
625 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
627 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
632 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
635 //Float_t pos[]={0,0,0};
636 //clu->GetPosition(pos);
637 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
638 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
640 if(iSupModMax > 1) {//sector 1
641 newPos[0] +=fMisalTransShift[3];//-=3.093;
642 newPos[1] +=fMisalTransShift[4];//+=6.82;
643 newPos[2] +=fMisalTransShift[5];//+=1.635;
644 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
648 newPos[0] +=fMisalTransShift[0];//+=1.134;
649 newPos[1] +=fMisalTransShift[1];//+=8.2;
650 newPos[2] +=fMisalTransShift[2];//+=1.197;
651 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
654 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
656 clu->SetPosition(newPos);
660 //__________________________________________________
661 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
663 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
664 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
667 Float_t fraction = 1.;
668 Float_t recalFactor = 1.;
672 Int_t iIphi = -1, iIeta = -1;
673 Int_t iSupMod = -1, iSupModMax = -1;
674 Int_t iphi = -1, ieta =-1;
675 Bool_t shared = kFALSE;
677 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
678 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
679 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
681 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
682 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
683 Int_t startingSM = -1;
685 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
686 absId = clu->GetCellAbsId(iDig);
687 fraction = clu->GetCellAmplitudeFraction(iDig);
688 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
689 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
690 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
692 if (iDig==0) startingSM = iSupMod;
693 else if(iSupMod != startingSM) areInSameSM = kFALSE;
695 eCell = cells->GetCellAmplitude(absId);
697 if(IsRecalibrationOn()) {
698 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
700 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
702 weight = GetCellWeight(eCell,clEnergy);
703 if(weight < 0) weight = 0;
704 totalWeight += weight;
705 weightedCol += ieta*weight;
706 weightedRow += iphi*weight;
708 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
712 Float_t xyzNew[]={0.,0.,0.};
713 if(areInSameSM == kTRUE) {
714 //printf("In Same SM\n");
715 weightedCol = weightedCol/totalWeight;
716 weightedRow = weightedRow/totalWeight;
717 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
720 //printf("In Different SM\n");
721 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
724 clu->SetPosition(xyzNew);
728 //____________________________________________________________________________
729 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
731 //re-evaluate distance to bad channel with updated bad map
733 if(!fRecalDistToBadChannels) return;
735 //Get channels map of the supermodule where the cluster is.
736 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
737 Bool_t shared = kFALSE;
738 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
739 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
742 Float_t minDist = 10000.;
745 //Loop on tower status map
746 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
747 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
748 //Check if tower is bad.
749 if(hMap->GetBinContent(icol,irow)==0) continue;
750 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
751 // iSupMod,icol, irow, icolM,irowM);
753 dRrow=TMath::Abs(irowM-irow);
754 dRcol=TMath::Abs(icolM-icol);
755 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
757 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
764 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
769 //The only possible combinations are (0,1), (2,3) ... (8,9)
770 if(iSupMod%2) iSupMod2 = iSupMod-1;
771 else iSupMod2 = iSupMod+1;
772 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
774 //Loop on tower status map of second super module
775 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
776 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
777 //Check if tower is bad.
778 if(hMap2->GetBinContent(icol,irow)==0) continue;
779 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
780 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
782 dRrow=TMath::Abs(irow-irowM);
785 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
788 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
791 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
792 if(dist < minDist) minDist = dist;
797 }// shared cluster in 2 SuperModules
799 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
800 cluster->SetDistanceToBadChannel(minDist);
804 //____________________________________________________________________________
805 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
807 //re-evaluate identification parameters with bayesian
809 if ( cluster->GetM02() != 0)
810 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
812 Float_t pidlist[AliPID::kSPECIESN+1];
813 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
815 cluster->SetPID(pidlist);
819 //____________________________________________________________________________
820 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
822 // Calculates new center of gravity in the local EMCAL-module coordinates
823 // and tranfers into global ALICE coordinates
824 // Calculates Dispersion and main axis
829 Float_t fraction = 1.;
830 Float_t recalFactor = 1.;
850 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
852 //Get from the absid the supermodule, tower and eta/phi numbers
853 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
854 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
856 //Get the cell energy, if recalibration is on, apply factors
857 fraction = cluster->GetCellAmplitudeFraction(iDigit);
858 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
859 if(IsRecalibrationOn()) {
860 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
862 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
864 if(cluster->E() > 0 && eCell > 0){
866 w = GetCellWeight(eCell,cluster->E());
874 dxx += w * etai * etai ;
876 dzz += w * phii * phii ;
878 dxz += w * etai * phii ;
882 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
885 //Normalize to the weight
891 AliError(Form("Wrong weight %f\n", wtot));
893 //Calculate dispersion
894 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
896 //Get from the absid the supermodule, tower and eta/phi numbers
897 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
898 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
900 //Get the cell energy, if recalibration is on, apply factors
901 fraction = cluster->GetCellAmplitudeFraction(iDigit);
902 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
903 if(IsRecalibrationOn()) {
904 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
906 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
908 if(cluster->E() > 0 && eCell > 0){
910 w = GetCellWeight(eCell,cluster->E());
914 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
917 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
920 //Normalize to the weigth and set shower shape parameters
921 if (wtot > 0 && nstat > 1) {
926 dxx -= xmean * xmean ;
927 dzz -= zmean * zmean ;
928 dxz -= xmean * zmean ;
929 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
930 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
934 cluster->SetM20(0.) ;
935 cluster->SetM02(0.) ;
939 cluster->SetDispersion(TMath::Sqrt(d)) ;
941 cluster->SetDispersion(0) ;
945 //__________________________________________________
946 void AliEMCALRecoUtils::FindMatches(AliVEvent *event)
948 //This function should be called before the cluster loop
949 //Before call this function, please recalculate the cluster positions
950 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
951 //Store matched cluster indexes and residuals
952 //It only works with ESDs, not AODs
954 fMatchedClusterIndex->Reset();
958 fMatchedClusterIndex->Set(100);
959 fResidualZ->Set(100);
960 fResidualR->Set(100);
965 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
967 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
968 if(!track || !IsAccepted(track)) continue;
970 Float_t dRMax = fCutR, dZMax = fCutZ;
972 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
973 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
975 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
976 if(!cluster->IsEMCAL()) continue;
977 cluster->GetPosition(clsPos); //Has been recalculated
978 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
979 emctrack->GetXYZ(trkPos);
980 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
981 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
994 fMatchedClusterIndex->AddAt(index,matched);
995 fResidualZ->AddAt(dZMax,matched);
996 fResidualR->AddAt(dRMax,matched);
1001 fMatchedClusterIndex->Set(matched);
1002 fResidualZ->Set(matched);
1003 fResidualR->Set(matched);
1005 //printf("Number of matched pairs: %d\n",matched);
1008 //__________________________________________________
1009 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1011 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1012 //Get the residuals dR and dZ for this cluster
1013 //It only works with ESDs, not AODs
1015 if( FindMatchedPos(index) >= 999 )
1017 AliDebug(2,"No matched tracks found!\n");
1022 dR = fResidualR->At(FindMatchedPos(index));
1023 dZ = fResidualZ->At(FindMatchedPos(index));
1026 //__________________________________________________
1027 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1029 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1030 //Returns if cluster has a match
1031 if(FindMatchedPos(index) < 999)
1036 //__________________________________________________
1037 UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
1039 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1040 //Returns the position of the match in the fMatchedClusterIndex array
1041 Float_t tmpR = fCutR;
1044 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1046 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1049 tmpR=fResidualR->At(i);
1051 AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1056 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1058 // Given a esd track, return whether the track survive all the cuts
1060 // The different quality parameter are first
1061 // retrieved from the track. then it is found out what cuts the
1062 // track did not survive and finally the cuts are imposed.
1064 UInt_t status = esdTrack->GetStatus();
1066 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1067 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1069 Float_t chi2PerClusterITS = -1;
1070 Float_t chi2PerClusterTPC = -1;
1071 if (nClustersITS!=0)
1072 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1073 if (nClustersTPC!=0)
1074 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1078 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1079 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1080 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1085 esdTrack->GetImpactParameters(b,bCov);
1086 if (bCov[0]<=0 || bCov[2]<=0) {
1087 AliDebug(1, "Estimated b resolution lower or equal zero!");
1088 bCov[0]=0; bCov[2]=0;
1091 Float_t dcaToVertexXY = b[0];
1092 Float_t dcaToVertexZ = b[1];
1093 Float_t dcaToVertex = -1;
1095 if (fCutDCAToVertex2D)
1096 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1098 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1102 Bool_t cuts[kNCuts];
1103 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1105 // track quality cuts
1106 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1108 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1110 if (nClustersTPC<fCutMinNClusterTPC)
1112 if (nClustersITS<fCutMinNClusterITS)
1114 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1116 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1118 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1120 if (fCutDCAToVertex2D && dcaToVertex > 1)
1122 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1124 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1127 //Require at least one SPD point + anything else in ITS
1128 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1132 for (Int_t i=0; i<kNCuts; i++)
1133 if (cuts[i]) {cut = kTRUE;}
1141 //__________________________________________________
1142 void AliEMCALRecoUtils::InitTrackCuts()
1144 //Intilize the track cut criteria
1145 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1146 //Also you can customize the cuts using the setters
1149 SetMinNClustersTPC(70);
1150 SetMaxChi2PerClusterTPC(4);
1151 SetAcceptKinkDaughters(kFALSE);
1152 SetRequireTPCRefit(kTRUE);
1155 SetRequireITSRefit(kTRUE);
1156 SetMaxDCAToVertexZ(2);
1157 SetDCAToVertex2D(kFALSE);
1160 //__________________________________________________
1161 void AliEMCALRecoUtils::Print(const Option_t *) const
1165 printf("AliEMCALRecoUtils Settings: \n");
1166 printf("Misalignment shifts\n");
1167 for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
1168 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1169 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1170 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1171 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1173 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1175 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1177 printf("Track cuts: \n");
1178 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1179 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1180 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1181 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1182 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);