VertexFit.h

Go to the documentation of this file.

#ifndef VertexFit_VertexFit_H
#define VertexFit_VertexFit_H
 
#include <vector>
#include "VertexFit/WTrackParameter.h"
#include "VertexFit/VertexParameter.h"
#include "VertexFit/VertexConstraints.h"
#include "VertexFit/TrackPool.h"
 
// NOTE: if you want to update the parameters of daughter tracks after vertex fitting,
// you should add the following code after invoking the interface "Fit()".
// VertexFit->Swim(n); 
// here, n is vertex number.
 
////////////////////////////////////////////////////////////////////////////////////////
//
// Attention: In this new Version, in order to find out why chisquare can be negative,
// so when successfully returned from Fit(int n), you will find chisq may be negative.
// You can Add this code:
//
// if (chisq < 0) return StatusCode::SUCCESS
//
// to avoid this situation  
//
////////////////////////////////////////////////////////////////////////////////////////
class VertexFit : public TrackPool 
{
public:
    // constructor & deconstructor
    static VertexFit *instance();
    ~VertexFit();
 
    // initialization, must be called before VertexFit each time
    void init();
 
    // add methods
    void AddVertex(int number, VertexParameter vpar, std::vector<int> lis);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6, int n7);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11);
    void AddVertex(int number, VertexParameter vpar, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12);
    void AddBeamFit(int number, VertexParameter vpar, int n);
    
    // build virtual particle, parameter number must be match with vertex number
    void BuildVirtualParticle(int number);
 
    // set iteration number and chisq cut
    void setIterNumber(const int niter = 10) { m_niter = niter; }
    void setChisqCut(const double chicut = 1000, const double chiter = 1.0e-3) { m_chicut = chicut;m_chiter=chiter; }
 
    // fit function
    bool Fit();
    bool Fit(int n);
    bool BeamFit(int n);
    
    // renew track parameters and convariance
    void Swim(int n){ vertexCovMatrix(n); swimVertex(n);}
    
    //
    // Fit Results
    //
 
    // chisq of fit
    double chisq() const {return m_chi;}
    double chisq(int n) const {return m_chisq[n];}
    
    // Pull distribution for 5-helix parameters.
    // Return true on success, the resulting pull will be stored in p;
    // Return false on failure.
    bool pull(int n, int itk, HepVector& p);
    
    // updated WTrack parameter in vertex fit
    HepLorentzVector pfit(int n) const {return wTrackInfit(n).p();}
    HepPoint3D xfit(int n) const {return wTrackInfit(n).x();}
    HepVector w(int n) const {return wTrackInfit(n).w();}
    HepSymMatrix Ew(int n) const {return wTrackInfit(n).Ew();}
    WTrackParameter wtrk(int n) const {return wTrackInfit(n);}
    
    // time counter
    HepVector cpu() const {return m_cpu;}
  
    // updated Vertex Parameter in vertex fit
    HepPoint3D vx(int n) const { return m_vpar_infit[n].vx();}
    HepVector Vx(int n) const { return m_vpar_infit[n].Vx();}
    HepSymMatrix Evx(int n) const { return m_vpar_infit[n].Evx();}
    double errorVx(int n, int i) const { return sqrt((m_vpar_infit[n].Evx())[i][i]);}
    VertexParameter vpar(int n) const {return m_vpar_infit[n];}
  
    // virtual particle from Vertex Fit
    WTrackParameter wVirtualTrack(int n) const {return m_virtual_wtrk[n];}
 
    // Users need not care below
private:
    // renew vertex constraints
    void UpdateConstraints(const VertexConstraints &vc);
    // vertex fit
    void fitVertex(int n);
    // renew covariance of input track parameters
    void vertexCovMatrix(int n);
    // extrapolate input track parameters to the vertex
    void swimVertex(int n);
    
    void fitBeam(int n);
    void swimBeam(int n);
 
    // vertex parameters before fit
    std::vector<VertexParameter> m_vpar_origin;
    // vertex parameters in fit
    std::vector<VertexParameter> m_vpar_infit;
    // vertex constraints
    std::vector<VertexConstraints> m_vc;
    // chisquare counter
    std::vector<double> m_chisq;
    double m_chi;   // total chisquare
    int m_nvtrk;    // number of WTracks in VertexFit
    // virtual particle in WTrackParameter format
    std::vector<WTrackParameter> m_virtual_wtrk;
  
    // track parameters storage and access
    // origin vertex and its covariance matrix
    HepVector m_xOrigin;
    HepSymMatrix m_xcovOrigin;
    HepSymMatrix m_xcovOriginInversed;
    inline HepVector xOrigin() const {return m_xOrigin;}
    inline void setXOrigin(const HepVector &x) { m_xOrigin = x;}
    inline HepSymMatrix xcovOrigin() const {return m_xcovOrigin;}
    inline void setXCovOrigin(const HepSymMatrix &v) {m_xcovOrigin = v;}
    inline HepSymMatrix xcovOriginInversed() const {return m_xcovOriginInversed;}
    inline void setXCovOriginInversed(const HepSymMatrix &v){m_xcovOriginInversed = v;}
 
    // vertex and covariance matrix in fit
    HepVector m_xInfit;
    HepSymMatrix m_xcovInfit;
    HepSymMatrix m_xcovInfitInversed;
    inline HepVector xInfit() const {return m_xInfit;}
    inline void setXInfit(const HepVector &x) {m_xInfit = x;}
    inline HepSymMatrix xcovInfit() const {return m_xcovInfit;}
    void setXCovInfit(const HepSymMatrix &v) {m_xcovInfit = v;}
    inline HepSymMatrix xcovInfitInversed() const {return m_xcovInfitInversed;}
    void setXCovInfitInversed(const HepSymMatrix &v) {m_xcovInfitInversed = v;}
  
    // track parameters and covariance matrice at initial
    HepVector m_pOrigin;
    HepSymMatrix m_pcovOrigin;
    inline HepVector pOrigin(int i) const { return m_pOrigin.sub(i*NTRKPAR+1, (i+1)*NTRKPAR); }
    inline void setPOrigin(int i, const HepVector &p) { m_pOrigin.sub(i*NTRKPAR+1, p); }
    inline HepSymMatrix pcovOrigin(int i) const { return m_pcovOrigin.sub(i*NTRKPAR+1, (i+1)*NTRKPAR); }
    inline void setPCovOrigin(int i, const HepSymMatrix &v) { m_pcovOrigin.sub(i*NTRKPAR+1,v); }
 
    // track parameters and covariance matrice in fit
    HepVector m_pInfit;
    HepSymMatrix m_pcovInfit;
    inline HepVector pInfit(int i) const { return m_pInfit.sub(i*NTRKPAR+1, (i+1)*NTRKPAR); }
    inline void setPInfit(int i, const HepVector &p) { m_pInfit.sub(i*NTRKPAR+1, p); }
    inline HepSymMatrix pcovInfit(int i) const { return m_pcovInfit.sub(i*NTRKPAR+1, (i+1)*NTRKPAR); }
    inline void setPCovInfit(int i, const HepSymMatrix &v) { m_pcovInfit.sub(i*NTRKPAR+1, v); }
 
    // some matrices convenient for calculation
    HepMatrix m_B;  // NCONSTR x NVTXPAR -- E
    inline HepMatrix vfB(int i) const {return m_B.sub(i*NCONSTR+1, (i+1)*NCONSTR, 1, NVTXPAR);}
    inline void setB(int i, const HepMatrix &e) {m_B.sub(i*NCONSTR+1, 1, e);}
    
    HepMatrix m_BT; // NVTXPAR x NCONSTR -- E.T()
    inline HepMatrix vfBT(int i) const {return m_BT.sub(1, NVTXPAR, i*NCONSTR+1, (i+1)*NCONSTR);}
    inline void setBT(int i, const HepMatrix &e) {m_BT.sub(1, i*NCONSTR+1, e);}
    
    HepMatrix m_A;  // NCONSTR x NTRKPAR -- D
    inline HepMatrix vfA(int i) const {return m_A.sub(i*NCONSTR+1, (i+1)*NCONSTR, i*NTRKPAR+1, (i+1)*NTRKPAR);}
    inline void setA(int i, const HepMatrix &d) {m_A.sub(i*NCONSTR+1, i*NTRKPAR+1, d);}
    
    HepMatrix m_AT; // NTRKPAR x NCONSTR -- D.T()
    inline HepMatrix vfAT(int i) const {return m_AT.sub(i*NTRKPAR+1, (i+1)*NTRKPAR, i*NCONSTR+1, (i+1)*NCONSTR);}
    inline void setAT(int i, const HepMatrix &d) {m_AT.sub(i*NTRKPAR+1, i*NCONSTR+1, d);}
    
    HepMatrix m_KQ; // NVTXPAR x NCONSTR -- Vx E.T() VD 
    inline HepMatrix vfKQ(int i) const {return m_KQ.sub(1, NVTXPAR, i*NCONSTR+1, (i+1)*NCONSTR);}
    inline void setKQ(int i, const HepMatrix &d) {m_KQ.sub(1, i*NCONSTR+1, d);}
    
    HepVector m_G;  // NCONSTR x 1  -- Dda + Edx + d
    inline HepVector vfG(int i) const {return m_G.sub(i*NCONSTR+1, (i+1)*NCONSTR);}
    inline void setG(int i, const HepVector &p) {m_G.sub(i*NCONSTR+1, p);}
    
    HepSymMatrix m_W;   // NCONSTR x NCONSTR -- VD
    inline HepSymMatrix vfW(int i) const {return m_W.sub(i*NCONSTR+1, (i+1)*NCONSTR);}
    inline void setW(int i, HepSymMatrix &m) {m_W.sub(i*NCONSTR+1, m);}
    
    HepMatrix m_E;  // NTRKPAR x NVTXPAR -- -Va0 D.T() VD E Vx
    inline HepMatrix vfE(int i) const {return m_E.sub(i*NTRKPAR+1, (i+1)*NTRKPAR, 1, NVTXPAR);}
    inline void setE(int i, const HepMatrix &p){m_E.sub(i*NTRKPAR+1, 1, p);}
 
    //NEW ADDED
    HepMatrix m_TRA;    // transform matrix from 7x1 to 6x1
    HepMatrix m_TRB;    // transform matrix from 6x1 to 7x1
    // convert HepVector 6x1 to 7x1
    // (px, py, pz, x, y, z) -> (px, py, pz, e, x, y, z)
    HepVector Convert67(const double &mass, const HepVector &p);    
    // convert HepVector 7x1 to 6x1
    // (px, py, pz, e, x, y, z) -> (px, py, pz, x, y, z)
    HepVector Convert76(const HepVector &p);
    //NEW ADDED
    
    // Singleton Design
    VertexFit();
    static VertexFit *m_pointer;
 
    int m_niter;    //number of iteration for vertex fitting
    double m_chicut;    //chisquare upper limit
    double m_chiter;
    HepVector m_cpu;
  
    static const int NTRKPAR;   //number of track parameters
    static const int NVTXPAR;   //number of vertex parameters
    static const int NCONSTR;   //number of vertex constraints
};
#endif  //VertexFit_VertexFit_H