00001 #include "VertexFit/FastVertexFit.h"
00002 #include "VertexFit/BField.h"
00003
00004 const double alpha = -0.00299792458;
00005 using namespace std;
00006
00007 FastVertexFit * FastVertexFit::m_pointer = 0;
00008
00009 FastVertexFit * FastVertexFit::instance() {
00010 if(m_pointer) return m_pointer;
00011 m_pointer = new FastVertexFit();
00012 return m_pointer;
00013 }
00014
00015 FastVertexFit::~FastVertexFit() {;}
00016
00017 FastVertexFit::FastVertexFit() {;}
00018
00019 void FastVertexFit::init() {
00020 m_D.clear();
00021 m_W.clear();
00022 m_DTWD.clear();
00023 m_xp.clear();
00024 m_chisq = 999;
00025 m_Vx = HepVector(3, 0);
00026 m_Evx = HepSymMatrix(3,0);
00027 }
00028
00029 void FastVertexFit::addTrack(const int number, const HepVector helix, const HepSymMatrix err) {
00030 int ifail;
00031 int ierr;
00032
00033
00034 HepSymMatrix Wc(2, 0);
00035 Wc[0][0] = err[0][0];
00036 Wc[0][1] = Wc[1][0] = err[0][3];
00037 Wc[1][1] = err[3][3];
00038 Wc = Wc.inverse(ifail);
00039
00040 HepVector p0(3, 0), v0(3, 0);
00041 double pxy = 1./fabs(helix[2]);
00042 p0[0] = 0 - pxy*sin(helix[1]);
00043 p0[1] = pxy*cos(helix[1]);
00044 p0[2] = pxy * helix[4];
00045 v0[0] = helix[0]*cos(helix[1]);
00046 v0[1] = helix[0]*sin(helix[1]);
00047 v0[2] = helix[3];
00048
00049 HepPoint3D vv0(v0[0], v0[1], v0[2]);
00050 double bField = VertexFitBField::instance()->getBFieldZRef();
00051
00052 int charge = (helix[2]>0 ? +1 :-1);
00053 double a = alpha * bField * charge;
00054 double T0 = sqrt((p0[0]+a*p0[1])*(p0[0]+a*p0[1])+(p0[1]-a*p0[0])*(p0[1]-a*p0[0]));
00055
00056 HepMatrix Dc(2, 3, 0);
00057 Dc[0][0] = (p0[1] - a*v0[0])/T0;
00058 Dc[0][1] = 0 - (p0[0] + a*v0[1])/T0;
00059 Dc[1][0] = 0 - (p0[2]/T0) * (p0[0] + a*v0[1])/T0;
00060 Dc[1][1] = 0 - (p0[2]/T0) * (p0[1] - a*v0[0])/T0;
00061 Dc[1][2] = 1;
00062
00063 m_D.push_back(Dc);
00064 m_W.push_back(Wc);
00065
00066 HepSymMatrix DTWD(3, 0);
00067 DTWD = Wc.similarity(Dc.T());
00068 HepVector qTrk(2, 0);
00069 qTrk[0] = helix[0];
00070 qTrk[1] = helix[3];
00071
00072
00073 m_DTWD.push_back(DTWD);
00074 m_xp.push_back(v0);
00075 }
00076
00077 bool FastVertexFit::Fit() {
00078 bool fitResult = false;
00079
00080 int ifail;
00081 HepSymMatrix total_DTWD(3, 0);
00082 HepVector total_xp(3, 0);
00083
00084 for(int i = 0; i < m_DTWD.size(); i++) {
00085 total_DTWD += m_DTWD[i];
00086 total_xp += m_DTWD[i]*m_xp[i];
00087 }
00088 m_Vx = total_DTWD.inverse(ifail) * total_xp;
00089 m_Evx = total_DTWD.inverse(ifail);
00090
00091 double chisq = 0;
00092 for(int i = 0; i < m_xp.size(); i++) {
00093 double chi2 = (m_DTWD[i].similarity((m_xp[i] - m_Vx).T()))[0][0];
00094 chisq += chi2;
00095 }
00096 m_chisq = chisq;
00097
00098 fitResult = true;
00099 return fitResult;
00100 }
00101
00102 HepVector FastVertexFit::Pull() const {
00103 HepVector pull(3, 0);
00104 pull[0] = m_Vx[0]/sqrt(m_Evx[0][0]);
00105 pull[1] = m_Vx[1]/sqrt(m_Evx[1][1]);
00106 pull[2] = m_Vx[2]/sqrt(m_Evx[2][2]);
00107 return pull;
00108 }
00109
00110
00111
00112