00001 #include "GaudiKernel/Bootstrap.h"
00002 #include "GaudiKernel/IJobOptionsSvc.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/PropertyMgr.h"
00005
00006 #include "EvtRecEvent/EvtRecVeeVertex.h"
00007 #include "DTagAlg/LocalKsSelector.h"
00008 #include "VertexFit/VertexFit.h"
00009 #include "VertexFit/SecondVertexFit.h"
00010 #include "VertexFit/IVertexDbSvc.h"
00011
00012 LocalKsSelector::LocalKsSelector()
00013 {
00014 IJobOptionsSvc* jobSvc;
00015 Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
00016
00017 PropertyMgr m_propMgr;
00018
00019 m_propMgr.declareProperty("KsMinMassCut", m_minMass = 0.470 );
00020 m_propMgr.declareProperty("KsMaxMassCut", m_maxMass = 0.528 );
00021 m_propMgr.declareProperty("KsMaxChisq", m_maxChisq = 100 );
00022
00023 m_propMgr.declareProperty("DoSecondaryVFit", m_doSecondaryVFit = false );
00024 m_propMgr.declareProperty("KsMaxVFitChisq", m_maxVFitChisq = 100 );
00025 m_propMgr.declareProperty("KsMinFlightSig", m_minFlightSig = 2.0 );
00026
00027 jobSvc->setMyProperties("LocalKsSelector", &m_propMgr);
00028 }
00029
00030 bool LocalKsSelector::operator() (CDKs& aKs) {
00031
00032 aKs.setUserTag(1);
00033 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( aKs.navKshort() );
00034
00035 if ( ks->vertexId() != 310 ) return false;
00036
00037 double mass = ks->mass();
00038 if ((mass <= m_minMass)||(mass >= m_maxMass)) return false;
00039 if ( ks->chi2() >= m_maxChisq ) return false;
00040
00041 if(mass < 0.487 || mass > 0.511)
00042 aKs.setUserTag(2);
00043
00044 if( !m_doSecondaryVFit ) return true;
00045
00046
00047
00048
00049
00050 EvtRecTrack* veeTrack1 = ks->pairDaughters().first;
00051 RecMdcKalTrack* veeKalTrack1 = veeTrack1->mdcKalTrack();
00052 WTrackParameter veeInitialWTrack1 = WTrackParameter(0.13957018, veeKalTrack1->getZHelix(), veeKalTrack1->getZError());
00053
00054 EvtRecTrack* veeTrack2 = ks->pairDaughters().second;
00055 RecMdcKalTrack* veeKalTrack2 = veeTrack2->mdcKalTrack();
00056 WTrackParameter veeInitialWTrack2 = WTrackParameter(0.13957018, veeKalTrack2->getZHelix(), veeKalTrack2->getZError());
00057
00058 VertexParameter wideVertex;
00059 HepPoint3D vWideVertex(0., 0., 0.);
00060 HepSymMatrix evWideVertex(3, 0);
00061
00062 evWideVertex[0][0] = 1.0e12;
00063 evWideVertex[1][1] = 1.0e12;
00064 evWideVertex[2][2] = 1.0e12;
00065
00066 wideVertex.setVx(vWideVertex);
00067 wideVertex.setEvx(evWideVertex);
00068
00069
00070 VertexFit* vtxfit = VertexFit::instance();
00071 vtxfit->init();
00072
00073
00074 vtxfit->AddTrack(0,veeInitialWTrack1);
00075 vtxfit->AddTrack(1,veeInitialWTrack2);
00076 vtxfit->AddVertex(0,wideVertex,0,1);
00077
00078
00079 vtxfit->Fit(0);
00080 vtxfit->Swim(0);
00081 vtxfit->BuildVirtualParticle(0);
00082
00083
00084 SecondVertexFit* svtxfit = SecondVertexFit::instance();
00085 svtxfit->init();
00086
00087
00088 VertexParameter beamSpot;
00089 HepPoint3D vBeamSpot(0., 0., 0.);
00090 HepSymMatrix evBeamSpot(3, 0);
00091
00092 IVertexDbSvc* vtxsvc;
00093 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00094 if(vtxsvc->isVertexValid()){
00095 double* dbv = vtxsvc->PrimaryVertex();
00096 double* vv = vtxsvc->SigmaPrimaryVertex();
00097 for (unsigned int ivx = 0; ivx < 3; ivx++){
00098 vBeamSpot[ivx] = dbv[ivx];
00099 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
00100 }
00101 }
00102 else{
00103 cout << "KSSELECTOR ERROR: Could not find a vertex from VertexDbSvc" << endl;
00104 return false;
00105 }
00106
00107 beamSpot.setVx(vBeamSpot);
00108 beamSpot.setEvx(evBeamSpot);
00109
00110 VertexParameter primaryVertex(beamSpot);
00111 svtxfit->setPrimaryVertex(primaryVertex);
00112
00113
00114 svtxfit->setVpar(vtxfit->vpar(0));
00115
00116
00117 svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
00118
00119
00120
00121 if( !svtxfit->Fit() ) return false;
00122
00123
00124 double vfitlength = svtxfit->decayLength();
00125 double vfiterror = svtxfit->decayLengthError();
00126 double vfitchi2 = svtxfit->chisq();
00127 double flightsig = 0;
00128 if( vfiterror != 0 )
00129 flightsig = vfitlength/vfiterror;
00130
00131
00132
00133 if( vfitchi2 > m_maxVFitChisq ) return false;
00134 if( flightsig < m_minFlightSig ) return false;
00135
00136 return true;
00137 }
00138
00139 LocalKsSelector ksSelector;