FastVertexFit.cxx

Go to the documentation of this file.

#include "VertexFit/FastVertexFit.h"
#include "VertexFit/BField.h" //FIXME
 
const double alpha = -0.00299792458;
using namespace std;
 
FastVertexFit * FastVertexFit::m_pointer = 0;
 
FastVertexFit * FastVertexFit::instance() {
  if(m_pointer) return m_pointer;
  m_pointer = new FastVertexFit();
  return m_pointer;
}
 
FastVertexFit::~FastVertexFit() {;}
 
FastVertexFit::FastVertexFit() {;}
 
void FastVertexFit::init() {
  m_D.clear();
  m_W.clear();
  m_DTWD.clear();
  m_xp.clear();
  m_chisq = 999;
  m_Vx = HepVector(3, 0);
  m_Evx = HepSymMatrix(3,0);
}
 
void FastVertexFit::addTrack(const int number, const HepVector helix, const HepSymMatrix err) {
  int ifail;
  int ierr;
  //HepSymMatrix err0(5, 0); //inverse of err;
  //err0 = err.inverse(ifail);
  HepSymMatrix Wc(2, 0);
  Wc[0][0] = err[0][0];
  Wc[0][1] = Wc[1][0] = err[0][3];
  Wc[1][1] = err[3][3];
  Wc = Wc.inverse(ifail);
 
  HepVector p0(3, 0), v0(3, 0);
  double pxy = 1./fabs(helix[2]);
  p0[0] = 0 - pxy*sin(helix[1]); //px
  p0[1] = pxy*cos(helix[1]);     //py
  p0[2] = pxy * helix[4];        //pz
  v0[0] = helix[0]*cos(helix[1]);//x
  v0[1] = helix[0]*sin(helix[1]);//y
  v0[2] = helix[3];              //z  
 
  HepPoint3D vv0(v0[0], v0[1], v0[2]);
  double bField = VertexFitBField::instance()->getBFieldZRef();
 
  int charge = (helix[2]>0 ? +1 :-1);
  double a = alpha * bField * charge;
  double T0 = sqrt((p0[0]+a*p0[1])*(p0[0]+a*p0[1])+(p0[1]-a*p0[0])*(p0[1]-a*p0[0]));
 
  HepMatrix Dc(2, 3, 0);
  Dc[0][0] = (p0[1] - a*v0[0])/T0;
  Dc[0][1] = 0 - (p0[0] + a*v0[1])/T0;
  Dc[1][0] = 0 - (p0[2]/T0) * (p0[0] + a*v0[1])/T0;
  Dc[1][1] = 0 - (p0[2]/T0) * (p0[1] - a*v0[0])/T0;
  Dc[1][2] = 1;
 
  m_D.push_back(Dc);
  m_W.push_back(Wc);
 
  HepSymMatrix DTWD(3, 0);
  DTWD = Wc.similarity(Dc.T());
  HepVector qTrk(2, 0);
  qTrk[0] = helix[0];
  qTrk[1] = helix[3];
  //HepVector xp(3, 0);
  //xp = Dc.inverse(ierr) * qTrk;
  m_DTWD.push_back(DTWD);
  m_xp.push_back(v0);
}
 
bool FastVertexFit::Fit() {
  bool fitResult = false;
 
  int ifail;
  HepSymMatrix total_DTWD(3, 0);
  HepVector total_xp(3, 0);
 
  for(int i = 0; i < m_DTWD.size(); i++) {
    total_DTWD += m_DTWD[i];  
    total_xp += m_DTWD[i]*m_xp[i]; 
  }
  m_Vx = total_DTWD.inverse(ifail) * total_xp;
  m_Evx = total_DTWD.inverse(ifail);
 
  double chisq = 0;
  for(int i = 0; i < m_xp.size(); i++) {
    double chi2 = (m_DTWD[i].similarity((m_xp[i] - m_Vx).T()))[0][0];  
    chisq += chi2;
  }
  m_chisq = chisq;
 
  fitResult = true;
  return fitResult;
}
 
HepVector FastVertexFit::Pull() const {
  HepVector pull(3, 0);
  pull[0] = m_Vx[0]/sqrt(m_Evx[0][0]);
  pull[1] = m_Vx[1]/sqrt(m_Evx[1][1]);
  pull[2] = m_Vx[2]/sqrt(m_Evx[2][2]);
  return pull;
}
 
//FastVertexFit::updateMatrices(const HepVector helix, const HepSymMatrix err) {
  
//}