/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Analysis/VertexFit/VertexFit-00-02-78/src/VertexConstraints.cxx

Go to the documentation of this file.

#include "VertexFit/VertexConstraints.h"
#include "VertexFit/TrackPool.h"
#include "VertexFit/BField.h"

const double VertexConstraints :: Alpha = -0.00299792458;
const int VertexConstraints::FixedVertex = 1;
const int VertexConstraints::CommonVertex = 2;

VertexConstraints::VertexConstraints()
{
        m_Ec.clear();
        m_Dc.clear();
        m_dc.clear();
        m_VD.clear();
        m_EVDE.clear();
        m_lambda0.clear();
        m_lambda.clear();
        m_covax.clear();
        m_ltrk.clear();
        m_type = 0;
}

void VertexConstraints::CommonVertexConstraints(std::vector<int> tlis) 
{
        m_ltrk = tlis;
        setType(commonVertex());
}

void VertexConstraints::FixedVertexConstraints(std::vector<int> tlis) 
{
        m_ltrk = tlis;
        setType(fixedVertex());
}

void VertexConstraints::UpdateConstraints(VertexParameter vpar, std::vector<WTrackParameter> wlis) 
{
        m_Ec.clear();
        m_Dc.clear();
        m_dc.clear();
        m_VD.clear();
        m_EVDE.clear();
        m_lambda0.clear();
        m_lambda.clear();
        m_covax.clear();

        switch(m_type) 
        {
                case FixedVertex: 
                {
                        WTrackParameter wtrk;
                        for (unsigned int i = 0; i < wlis.size(); i++) 
                        {
                                wtrk = wlis[i];
                                HepLorentzVector p = wtrk.p();
                                HepPoint3D x = wtrk.x();
                                HepPoint3D delx = vpar.vx() - x;

                                double afield = VertexFitBField::instance()->getCBz(vpar.Vx(), wtrk.X());
                                double a = afield * (0.0+wtrk.charge());
                                double J = a * (delx.x()*p.px() + delx.y()*p.py())/p.perp2();

                                J = std::min(J, 1-1e-4);
                                J = std::max(J, -1+1e-4);
                                double Rx = delx.x() - 2*p.px() * (delx.x()*p.px() + delx.y()*p.py()) / p.perp2();
                                double Ry = delx.y() - 2*p.py() * (delx.x()*p.px() + delx.y()*p.py()) / p.perp2();
                                double S = 1.0/sqrt(1-J*J)/p.perp2();
             
                                // dc
                                HepVector dc(2, 0);
                                dc[0] = delx.y() * p.px() - delx.x() * p.py() - 0.5 * a * (delx.x() * delx.x() + delx.y() * delx.y());
                                dc[1] = delx.z() - p.pz()/a*asin(J);
                                m_dc.push_back(dc);
                                // Ec
                                HepMatrix Ec(2, 3, 0);
                                m_Ec.push_back(Ec);
                                // Dc
                                HepMatrix Dc(2, 6, 0);
                                Dc[0][0] = delx.y();
                                Dc[0][1] = -delx.x();
                                Dc[0][2] = 0;
                                Dc[0][3] = p.py() + a * delx.x();
                                Dc[0][4] = -p.px() + a * delx.y();
                                Dc[0][5] = 0;
                                Dc[1][0] = -p.pz() * S * Rx;
                                Dc[1][1] = -p.pz() * S * Ry;
                                Dc[1][2] = - asin(J) / a;
                                Dc[1][3] = p.px() * p.pz() * S;
                                Dc[1][4] = p.py() * p.pz() * S;
                                Dc[1][5] = -1.0;
                                m_Dc.push_back(Dc);
                                // VD
                                HepSymMatrix VD(2, 0);
                                // int ifail;
                                // VD = ((wtrk.Ew()).similarity(Dc)).inverse(ifail); // 2x2 matrix (Dc Ew Dc.T())
                                m_VD.push_back(VD);
                                // EVDE
                                HepSymMatrix EVDE(3, 0);
                                m_EVDE.push_back(EVDE);
                                // lambda0
                                HepVector lambda0(2, 0);
                                m_lambda0.push_back(lambda0);
                                // lambda
                                HepVector lambda(2, 0);
                                m_lambda.push_back(lambda);
                                // covax
                                HepMatrix covax(6, 3, 0);
                                m_covax.push_back(covax);
                                break;
                        }
                }
                case CommonVertex:
                default: 
                {
                        WTrackParameter wtrk;
                        for (unsigned int i = 0; i < wlis.size(); i++)
                        {
                                wtrk = wlis[i];
                                HepLorentzVector p = wtrk.p();
                                HepPoint3D x = wtrk.x();
                                HepPoint3D delx = vpar.vx() - x;
                                if (wtrk.charge() == 0) 
                                {
                                        //dc
                                        HepVector dc(2,0);
                                        dc[0] = p.pz()*delx.x() - p.px()*delx.z();
                                        dc[1] = p.pz()*delx.y() - p.py()*delx.z();
                                        m_dc.push_back(dc);
                                        //Ec
                                        HepMatrix Ec(2,3,0);
                                        Ec[0][0] = p.pz();
                                        Ec[0][1] = 0;
                                        Ec[0][2] = -p.px();
                                        Ec[1][0] = 0;
                                        Ec[1][1] = p.pz();
                                        Ec[1][2] = -p.py();
                                        m_Ec.push_back(Ec);
                                        //Dc
                                        HepMatrix Dc(2,6,0);
                                        Dc[0][0] = -delx.z();
                                        Dc[0][1] = 0;
                                        Dc[0][2] = delx.x();
                                        Dc[0][3] = -p.pz();
                                        Dc[0][4] = 0;
                                        Dc[0][5] = p.px();
                                        Dc[1][0] = 0;
                                        Dc[1][1] = -delx.z();
                                        Dc[1][2] = delx.y();
                                        Dc[1][3] = 0;
                                        Dc[1][4] = -p.pz();
                                        Dc[1][5] = p.py();
                                        m_Dc.push_back(Dc);
                                } 
                                else 
                                {
                                        double afield = VertexFitBField::instance()->getCBz(vpar.Vx(), wtrk.X());
                                        double a = afield * (0.0+wtrk.charge());
                                        double J = a * (delx.x()*p.px() + delx.y()*p.py())/p.perp2();
                                        J = std::min(J, 1-1e-4);
                                        J = std::max(J, -1+1e-4);
                                        double Rx = delx.x() - 2*p.px() * (delx.x()*p.px() + delx.y()*p.py()) / p.perp2();
                                        double Ry = delx.y() - 2*p.py() * (delx.x()*p.px() + delx.y()*p.py()) / p.perp2();
                                        double S = 1.0/sqrt(1-J*J)/p.perp2();
                                        // dc
                                        HepVector dc(2, 0);
                                        dc[0] = delx.y() * p.px() - delx.x() * p.py() - 0.5 * a * (delx.x() * delx.x() + delx.y() * delx.y());
                                        dc[1] = delx.z() - p.pz()/a*asin(J);
                                        m_dc.push_back(dc);
                                        // Ec
                                        HepMatrix Ec(2, 3, 0);
                                        Ec[0][0] = -p.py()- a * delx.x();
                                        Ec[0][1] = p.px() - a * delx.y();
                                        Ec[0][2] = 0;
                                        Ec[1][0] = -p.px() * p.pz() * S ;
                                        Ec[1][1] = -p.py() * p.pz() * S ;
                                        Ec[1][2] = 1.0;
                                        m_Ec.push_back(Ec);
                                        // Dc
                                        HepMatrix Dc(2, 6, 0);
                                        Dc[0][0] = delx.y();
                                        Dc[0][1] = -delx.x();
                                        Dc[0][2] = 0;
                                        Dc[0][3] = p.py() + a * delx.x();
                                        Dc[0][4] = -p.px() + a * delx.y();
                                        Dc[0][5] = 0;
                                        Dc[1][0] = -p.pz() * S * Rx;
                                        Dc[1][1] = -p.pz() * S * Ry;
                                        Dc[1][2] = - asin(J) / a;
                                        Dc[1][3] = p.px() * p.pz() * S;
                                        Dc[1][4] = p.py() * p.pz() * S;
                                        Dc[1][5] = -1.0;
                                        m_Dc.push_back(Dc);
                                }
                                // VD
                                HepSymMatrix VD(2, 0);
                                m_VD.push_back(VD);
                                // EVDE
                                HepSymMatrix EVDE(3, 0);
                                m_EVDE.push_back(EVDE);
                                // lambda0
                                HepVector lambda0(2, 0);
                                m_lambda0.push_back(lambda0);
                                // lambda
                                HepVector lambda(2, 0);
                                m_lambda.push_back(lambda);
                                // covax
                                HepMatrix covax(6, 3, 0);
                                m_covax.push_back(covax);
                        }
                        break;
                }
        }
}

---