/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Analysis/VertexFit/VertexFit-00-02-78/src/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;
        m_factor     = 1.000;
}

bool SecondVertexFit::Fit() 
{
        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 afield = m_factor * 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) 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;
}

---