SecondVertexFit.cxx

Go to the documentation of this file.

#include "VertexFit/SecondVertexFit.h"
#include "VertexFit/BField.h"
 
//----------------------------------------------------
// In Second Vertex Fit:
//
// Fitting Parameters
// (px, py, pz, E, xd, yd, zd, xp, yp, zp)
//
// Constraints (charge = 0)
//
// xp - xd + px/m * ct = 0
// yp - yd + py/m * ct = 0
// zp - zd + pz/m * ct = 0
//
//     ( 0 0 0 0 -1 0 0 1 0 0 )
// D = ( 0 0 0 0 0 -1 0 0 1 0 )
//     ( 0 0 0 0 0 0 -1 0 0 1 )
//
//     ( px/m )
// E = ( py/m )
//     ( pz/m )
//
//     ( xp - xd + px/m * ct )
// d = ( yp - yd + py/m * ct )
//     ( zp - zd + pz/m * ct )
//
// Constraints (charge != 0)
//
// xp - xd + px/a * sin(a ctau/m) + py/a(1-cos(a ctau/m)) = 0
// yp - yd + py/a * sin(a ctau/m) - px/a(1-cos(a ctau/m)) = 0
// zp - zd + pz/m ctau = 0
//     ( 0 0 0 0 -1 0 0 1 0 0 )
// D = ( 0 0 0 0 0 -1 0 0 1 0 )
//     ( 0 0 0 0 0 0 -1 0 0 1 )
//
//     ( px/m cos(a * ct/m) + py/m * sin (a * ct/m))
// E = ( py/m cos(a * ct/m) - px/m * sin (a * ct/m))
//     ( pz/m )
//
//----------------------------------------------------
 
SecondVertexFit *SecondVertexFit::m_pointer = 0;
 
SecondVertexFit *SecondVertexFit::instance() 
{
    if (m_pointer) 
        return m_pointer;
    m_pointer = new SecondVertexFit();
    return m_pointer;
}
 
SecondVertexFit::SecondVertexFit()
{
    HepVector vx(3,0);
    m_vpar_primary.setVx(vx);
    HepSymMatrix evx(3,0);
    evx[0][0] = 0.1 * 0.1;
    evx[1][1] = 0.1 * 0.1;
    evx[2][2] = 1.5 * 1.5;
    m_vpar_primary.setEvx(evx);
}
 
SecondVertexFit::~SecondVertexFit() 
{
//  if(m_pointer) delete m_pointer;
}
 
void SecondVertexFit::init() 
{
    clearWTrackOrigin();
    clearWTrackInfit();
    clearWTrackList();
    m_vpar_secondary = VertexParameter();
    m_lxyz       = 0;
    m_lxyz_error = 0;
    m_p4par      = HepLorentzVector(0, 0, 0, 0);
    m_crxyz      = HepVector(3, 0);
    m_chisq      = 9999;
    m_wtrk       = WTrackParameter();
    m_niter      = 10;
    m_chicut     = 500;
    m_chiter     = 1.0e-2;
}
 
bool SecondVertexFit::Fit() 
{
//  const double alpha = -0.00299792458;
//  const double field = 1.0;
//  double field_p = VertexFitBField::instance()->getBFieldZ(m_vpar_primary.vx());
//  double field_s = VertexFitBField::instance()->getBFieldZ(m_vpar_secondary.vx());
//  double field = (field_p + field_s) * 0.5 ;
//  std::cout << "field in SeconVertexFit= " << field << std::endl;
 
    bool okfit = false;
 
    HepVector aOrigin(10, 0);
    HepVector aInfit(10, 0);
    HepSymMatrix VaOrigin(10, 0);
    HepSymMatrix VaInfit(10, 0);
    aOrigin.sub(1, wTrackOrigin(0).w());
    aOrigin.sub(8, m_vpar_primary.Vx());
    VaOrigin.sub(1, wTrackOrigin(0).Ew());
    VaOrigin.sub(8, m_vpar_primary.Evx());
    HepVector ctOrigin(1, 0);
    HepVector ctInfit(1, 0);
    HepSymMatrix Vct(1, 0);
    aInfit  = aOrigin;
    ctInfit = ctOrigin;
 
    std::vector<double> chisq;
    chisq.clear();
    double chi2 = 999;
    for(int it = 0; it < m_niter; it++)
    {
        HepMatrix D(3, 10, 0);
        HepLorentzVector p4par = HepLorentzVector(aInfit[0], aInfit[1], aInfit[2], aInfit[3]);
        HepMatrix E(3,1,0);
        HepVector d(3, 0);
        if (wTrackOrigin(0).charge() == 0) 
        {
            D[0][4] = -1.0;
            D[0][7] = 1.0;
            D[1][5] = -1.0;
            D[1][8] = 1.0;
            D[2][6] = -1.0;
            D[2][9] = 1.0;
 
            E[0][0] = p4par.px()/p4par.m();
            E[1][0] = p4par.py()/p4par.m();
            E[2][0] = p4par.pz()/p4par.m();
 
            d[0] = aInfit[7]-aInfit[4]+ctInfit[0]*p4par.px()/p4par.m();
            d[1] = aInfit[8]-aInfit[5]+ctInfit[0]*p4par.py()/p4par.m();
            d[2] = aInfit[9]-aInfit[6]+ctInfit[0]*p4par.pz()/p4par.m();
        }
        else 
        {
            double afield = VertexFitBField::instance()->getCBz(m_vpar_primary.Vx(), m_vpar_secondary.Vx());
            double a = afield * wTrackOrigin(0).charge();
            D[0][4] = -1.0;
            D[0][7] = 1.0;
            D[1][5] = -1.0;
            D[1][8] = 1.0;
            D[2][6] = -1.0;
            D[2][9] = 1.0;
 
            E[0][0] = p4par.px() / p4par.m() * cos(a * ctInfit[0] / p4par.m()) + p4par.py() / p4par.m() * sin(a * ctInfit[0] / p4par.m());
            E[1][0] = p4par.py() / p4par.m() * cos(a * ctInfit[0] / p4par.m()) - p4par.px() / p4par.m() * sin(a * ctInfit[0] / p4par.m());
            E[2][0] = p4par.pz() / p4par.m();
 
            d[0] = aInfit[7] - aInfit[4]+p4par.px()/a * sin(a * ctInfit[0]/p4par.m()) + p4par.py()/a*(1-cos(a*ctInfit[0]/p4par.m()));
            d[1] = aInfit[8] - aInfit[5]+p4par.py()/a * sin(a * ctInfit[0]/p4par.m()) - p4par.px()/a*(1-cos(a*ctInfit[0]/p4par.m()));
            d[2] = aInfit[9] - aInfit[6]+ctInfit[0]*p4par.pz()/p4par.m();
        }
        
        HepSymMatrix VD(3, 0);
        HepVector dela0(10, 0);
        HepVector lambda0(3, 0);
        HepVector delct(1, 0);
        HepVector lambda(3, 0);
        int ifail;
 
        VD      = (VaOrigin.similarity(D)).inverse(ifail);
        dela0   = aOrigin - aInfit;
        lambda0 = VD*(D*dela0 + d);
        Vct     = (VD.similarity(E.T())).inverse(ifail);
        delct   = -(Vct * E.T()) * lambda0;
        ctInfit = ctInfit + delct;
        lambda  = lambda0 + (VD * E) * delct;
        aInfit  = aOrigin - (VaOrigin * D.T()) * lambda;
        chi2    = dot(lambda, D*dela0 + d);
        VaInfit = VaOrigin - (VD.similarity(D.T())).similarity(VaOrigin);
        VaInfit = VaInfit + (((Vct.similarity(E)).similarity(VD)).similarity(D.T())).similarity(VaOrigin);
 
        chisq.push_back(chi2);
 
        if(it > 0) 
        {
            double delchi = chisq[it] - chisq[it-1];
            if (fabs(delchi) < m_chiter) 
            {
//              std::cout <<"Fit(): pass all " << it <<" , " << delchi << std::endl;
                break;
            }
        }
    }
    if (chi2 < 0 || chi2 > m_chicut)
        return okfit;
 
    HepLorentzVector p4par = HepLorentzVector(aInfit[0], aInfit[1], aInfit[2], aInfit[3]);
    m_ctau       = ctInfit[0];
    m_ctau_error = sqrt(Vct[0][0]);
    m_lxyz       = ctInfit[0] * p4par.rho() / p4par.m();
    m_lxyz_error = sqrt(Vct[0][0]) * p4par.rho() / p4par.m();
    m_chisq      = chi2;
    m_p4par      = p4par;
    for(int i = 0; i < 3; i++)
        m_crxyz[i] = aInfit[4+i];
    HepVector w(7, 0);
    HepSymMatrix Ew(7, 0);
    for(int i = 0; i < 7; i++) 
    {
        w[i] = aInfit[i];
        for(int j = 0; j < 7; j++) 
        {
            Ew[i][j] = VaInfit[i][j];
        }
    }
    m_wtrk.setW(w);
    m_wtrk.setEw(Ew);
    m_wtrk.setCharge(wTrackOrigin(0).charge());
    okfit = true;
    return okfit;
}