VertexFitRefine.cxx

Go to the documentation of this file.

/* <===<===<===<===<===<===<===<===<===~===>===>===>===>===>===>===>===>===>
 * File Name:    VertexFitRefine.cxx
 * Author:       Hao-Kai SUN
 * Created:      2021-09-08 Wed 00:15:40 CST
 * <<=====================================>>
 * Last Updated: 2023-09-29 Thu 14:35:11 CST
 *           By: Hao-Kai SUN
 *     Update #: 153
 * ============================== CODES ==============================>>> */
#include "VertexFitRefine/VertexFitRefine.h"
 
namespace VTXPDGM {
  const double& electron = 0.0005109989461;
  const double& muon     = 0.1056583745;
  const double& pion     = 0.13957039;
  const double& kaon     = 0.493677;
  const double& proton   = 0.938272081;
 
  const double empikp[5]  = { electron, muon, pion, kaon, proton };
  const double empikp2[5] = { electron * electron, muon* muon, pion* pion, kaon* kaon,
                              proton* proton };
} // namespace VTXPDGM
 
VertexFitRefine* VertexFitRefine::m_instance = nullptr;
 
// const double VertexFitRefine::beampipe  = 3.111;
const double VertexFitRefine::obeampipe = 3.370; // + 0.0130 * 1;
// 6.3000 + 0.0130*3 cm inner mdc chamber wall
// const double VertexFitRefine::innerwall  = 6.290 - 0.0130 * 1;
// const double VertexFitRefine::oinnerwall = 6.425 + 0.0130 * 1;
 
VertexFitRefine::VertexFitRefine() {}
 
VertexFitRefine* VertexFitRefine::instance() {
  if ( m_instance == nullptr ) m_instance = new VertexFitRefine();
  return m_instance;
}
 
void VertexFitRefine::init() {
  m_trkIdxOrigin.clear();
  m_tracksOrigin.clear();
  m_trkPidOrigin.clear();
  m_wtrkInfit.clear();
  m_p4Infit.clear();
  m_x3Infit.clear();
  m_helices.clear();
  vtxfit = VertexFit::instance();
  vtxfit->init();
  vtxext = VertexExtrapolate::instance();
  m_vtxsOrigin.clear();
  thePath = 0;
}
 
VertexFitRefine::~VertexFitRefine() { vtxfit->init(); }
 
void VertexFitRefine::AddTrack( const int index, const WTrackParameter wtrk ) {
  m_trkIdxOrigin.push_back( index );
  m_tracksOrigin.push_back( nullptr );
  m_trkPidOrigin.push_back( RecMdcKalTrack::null );
  int idx = m_trkIdxOrigin.size() - 1;
  if ( index != idx )
  {
    std::cerr << "TrackPool: wrong track index " << index << ", " << m_trkIdxOrigin.size()
              << std::endl;
  }
 
  m_tracksOrigin[idx]->setPidType( m_trkPidOrigin[idx] );
  m_wtrkInfit.push_back( wtrk );
  m_helices.push_back( HepVector( 5, 0 ) );
  vtxfit->AddTrack( m_trkIdxOrigin[idx], m_wtrkInfit[idx] );
}
 
void VertexFitRefine::AddTrack( const int index, RecMdcKalTrack* p,
                                const RecMdcKalTrack::PidType pid ) {
  m_trkIdxOrigin.push_back( index );
  m_tracksOrigin.push_back( p );
  m_trkPidOrigin.push_back( pid );
  int idx = m_trkIdxOrigin.size() - 1;
  if ( index != idx )
  {
    std::cerr << "TrackPool: wrong track index " << index << ", " << m_trkIdxOrigin.size()
              << std::endl;
  }
 
  m_tracksOrigin[idx]->setPidType( m_trkPidOrigin[idx] );
  m_wtrkInfit.push_back( WTrackParameter( VTXPDGM::empikp[m_trkPidOrigin[idx]],
                                          m_tracksOrigin[idx]->helix(),
                                          m_tracksOrigin[idx]->err() ) );
  vtxfit->AddTrack( m_trkIdxOrigin[idx], m_wtrkInfit[idx] );
}
 
bool VertexFitRefine::Fit() {
  HepPoint3D vx( 0.0, 0.0, 0.0 );
 
  for ( uint i = 0; i != m_trkIdxOrigin.size(); ++i )
  {
    m_p4Infit.push_back( m_wtrkInfit[i].p() );
    m_x3Infit.push_back( m_wtrkInfit[i].x() );
  }
 
  //// setup the initial vertex
  HepPoint3D vWideVertex( 0., 0., 0. );
 
  HepSymMatrix evWideVertex( 3, 0 );
  evWideVertex[0][0] = 1.0E12;
  evWideVertex[1][1] = 1.0E12;
  evWideVertex[2][2] = 1.0E12;
 
  VertexParameter wideVertex;
  wideVertex.setVx( vWideVertex );
  wideVertex.setEvx( evWideVertex );
 
  if ( vtxfit->m_vpar_infit.size() == 0 )
  {
    std::cerr << "Not set Vertex?" << std::endl;
    return false;
  }
 
  HepPoint3D   ZDP  = vtxfit->vx( 0 );
  HepPoint3D   ZDPE = HepPoint3D( 9999.9, 9999.9, 9999.9 );
  bool         ZFit = false;
  HepVector    ZVx  = vtxfit->Vx( 0 );
  HepSymMatrix ZEVx = evWideVertex;
 
  //// do the fit
  if ( vtxfit->Fit( 0 ) )
  {
    vtxfit->Swim( 0 );
    vtxfit->BuildVirtualParticle( 0 );
 
    for ( uint i = 0; i != m_trkIdxOrigin.size(); ++i )
    {
      m_p4Infit[i] = vtxfit->pfit( i );
      m_x3Infit[i] = vtxfit->xfit( i );
    }
 
    ZDP  = vtxfit->vx( 0 );
    ZVx  = vtxfit->Vx( 0 );
    ZEVx = vtxfit->Evx( 0 );
    ZDPE.set( vtxfit->errorVx( 0, 0 ), vtxfit->errorVx( 0, 1 ), vtxfit->errorVx( 0, 2 ) );
    ZFit = true;
  }
 
  if ( ZFit )
  {
    if ( ZDP.perp() > obeampipe )
    {
      //// initialize the fitter
      vtxfit->init();
 
      for ( uint i = 0; i != m_trkIdxOrigin.size(); ++i )
      {
        m_tracksOrigin[i]->setPidType( m_trkPidOrigin[i] );
        if ( m_trkPidOrigin[i] != RecMdcKalTrack::null )
        {
          m_wtrkInfit[i] =
              WTrackParameter( VTXPDGM::empikp[m_trkPidOrigin[i]], m_tracksOrigin[i]->fhelix(),
                               m_tracksOrigin[i]->ferr() );
        }
        vtxfit->AddTrack( i, m_wtrkInfit[i] );
      }
 
      wideVertex.setVx( ZVx );
      // wideVertex.setEvx(ZEVx * 25.0); // 5 * sigma square, not good
      wideVertex.setEvx( evWideVertex );
 
      //// add the daughters
      vtxfit->AddVertex( 0, wideVertex, m_trkIdxOrigin );
      //// do the fit
      if ( vtxfit->Fit( 0 ) )
      {
        vtxfit->Swim( 0 );
        vtxfit->BuildVirtualParticle( 0 );
        vx      = vtxfit->vx( 0 );
        ZVx     = vtxfit->Vx( 0 );
        ZEVx    = vtxfit->Evx( 0 );
        thePath = 3;
      }
      else
      {
        vx      = ZDP;
        thePath = 2;
      }
    }
    else
    {
      vx      = ZDP;
      thePath = 1;
    }
  }
  else
  { // initial ZFit failed.
    //// initialize the fitter
    vtxfit->init();
 
    for ( uint i = 0; i != m_trkIdxOrigin.size(); ++i )
    {
      m_tracksOrigin[i]->setPidType( m_trkPidOrigin[i] );
      if ( m_trkPidOrigin[i] != RecMdcKalTrack::null )
      {
        m_wtrkInfit[i] =
            WTrackParameter( VTXPDGM::empikp[m_trkPidOrigin[i]], m_tracksOrigin[i]->fhelix(),
                             m_tracksOrigin[i]->ferr() );
      }
      vtxfit->AddTrack( i, m_wtrkInfit[i] );
    }
 
    if ( vtxfit->Fit( 0 ) )
    {
      vtxfit->Swim( 0 );
      vtxfit->BuildVirtualParticle( 0 );
      vx      = vtxfit->vx( 0 );
      ZVx     = vtxfit->Vx( 0 );
      ZEVx    = vtxfit->Evx( 0 );
      thePath = 4;
    }
    else
    {
      thePath = 5;
      return false;
    }
  }
  //// initialize the fitter
  vtxfit->init();
 
  for ( uint i = 0; i != m_trkIdxOrigin.size(); ++i )
  {
    m_tracksOrigin[i]->setPidType( m_trkPidOrigin[i] );
    if ( m_trkPidOrigin[i] != RecMdcKalTrack::null )
    {
      vtxext->KalFitExt( vx, m_tracksOrigin[i], m_trkPidOrigin[i] );
      m_wtrkInfit[i] = WTrackParameter( VTXPDGM::empikp[m_trkPidOrigin[i]],
                                        vtxext->getHelixVector(), vtxext->getErrorMatrix() );
    }
    m_p4Infit[i] = m_wtrkInfit[i].p();
    m_x3Infit[i] = m_wtrkInfit[i].x();
 
    vtxfit->AddTrack( i, m_wtrkInfit[i] );
  }
 
  wideVertex.setVx( vx );
  // wideVertex.setEvx(ZEVx * 25.0); // 5 * sigma square, not good, why?
  wideVertex.setEvx( evWideVertex );
  //// add the daughters
  vtxfit->AddVertex( 0, wideVertex, m_trkIdxOrigin );
 
  //// do the fit
  if ( vtxfit->Fit( 0 ) )
  {
    vtxfit->Swim( 0 );
    vtxfit->BuildVirtualParticle( 0 );
  }
  else { thePath = 6; }
 
  return true;
}
 
/* ===================================================================<<< */
/* =================== VertexFitRefine.cxx ends here ==================== */