MdcTrackUtil.cxx

Go to the documentation of this file.

#include "TrackUtil/MdcTrackUtil.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include "EventModel/EventHeader.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "RawEvent/RawDataUtil.h"
// Mointe-Carlo data
#include "Identifier/MdcID.h"
#include "MdcRawEvent/MdcDigi.h"
 
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
//  backwards compatblty wll be enabled ONLY n CLHEP 1.9
typedef HepGeom::Point3D<double> HepPoint3D;
#endif
 
// MDC reconstructed data
#include "MdcRecEvent/RecMdcHit.h"
#include "MdcRecEvent/RecMdcTrack.h"
 
// MDC Geometory
#include "MdcGeomSvc/IMdcGeomSvc.h"
 
using namespace std;
using namespace Event;
 
MdcTrackUtil* MdcTrackUtil::_myself = 0;
 
MdcTrackUtil* MdcTrackUtil::instance() {
  if ( 0 == _myself ) { _myself = new MdcTrackUtil(); }
  return _myself;
}
 
MdcTrackUtil::MdcTrackUtil() {
  // Initalze magnetic flied
  m_pIMF = service( "MagneticFieldSvc" );
  if ( !m_pIMF ) { std::cout << " ERROR Unable to open Magnetic field service " << std::endl; }
  // get Bz for Check TEMP, Bz may be changed by run
  double gaussToTesla = 1000.;
  Bz                  = m_pIMF->getReferField() * gaussToTesla;
 
  // Initialize MdcGeomSvc
  m_mdcGeomSvc = service( "MdcGeomSvc" );
  if ( !m_mdcGeomSvc ) { std::cout << " FATAL Could not load MdcGeomSvc!  " << std::endl; }
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
int MdcTrackUtil::nLayerTrackPassed( const HepVector helix ) {
  double helixParam[5];
  for ( int i = 0; i < 5; ++i ) helixParam[i] = helix[i];
 
  return nLayerTrackPassed( helixParam );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
int MdcTrackUtil::nLayerTrackPassed( const double helix[5] ) {
  int nLayer = 0;
 
  for ( unsigned iLayer = 0; iLayer < 43; iLayer++ )
  {
    // flightLength is the path length of track in the x-y plane
    // guess flightLength by the radius in middle of the wire.
    double rMidLayer    = m_mdcGeomSvc->Layer( iLayer )->Radius();
    double flightLength = rMidLayer;
 
    HepPoint3D pivot( 0., 0., 0. );
    double     dz    = helix[3];
    double     c     = CLHEP::c_light * 100.; // unit from m/s to cm/s
    double     alpha = 1 / ( c * Bz );        //~333.567
    double     kappa = helix[2];
    double     rc    = ( -1. ) * alpha / kappa;
    // std::cout<<"MdcTrackUtil alpha   "<<alpha<<std::endl;
    double tanl = helix[4];
    double phi0 = helix[1];
    double phi  = flightLength / rc + phi0; // turning angle
    double z    = pivot.z() + dz - ( alpha / kappa ) * tanl * phi;
 
    double layerHalfLength = m_mdcGeomSvc->Layer( iLayer )->Length() / 2.;
 
    // std::cout<<"MdcTrackUtil length  "<<layerHalfLength<<std::endl;
 
    if ( fabs( z ) < fabs( layerHalfLength ) ) ++nLayer;
  }
 
  return nLayer;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
HepVector MdcTrackUtil::patRecPar2BesPar( const HepVector& helixPar ) {
  HepVector helix( 5, 0 );
  double    d0    = -helixPar[0]; // cm
  double    phi0  = helixPar[1] + CLHEP::halfpi;
  double    omega = Bz * helixPar[2] / -333.567;
  double    z0    = helixPar[3]; // cm
  double    tanl  = helixPar[4];
  helix[0]        = d0;
  helix[1]        = phi0;
  helix[2]        = omega;
  helix[3]        = z0;
  helix[4]        = tanl;
  // std::cout<<"helix   "<<helix<<std::endl;
  return helix;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
HepSymMatrix MdcTrackUtil::patRecErr2BesErr( const HepSymMatrix& err ) {
  // error matrix
  // std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
  // V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
  // int n = err.num_row();
  HepSymMatrix mS( err.num_row(), 0 );
  mS[0][0]         = -1.; // dr0=-d0
  mS[1][1]         = 1.;
  mS[2][2]         = Bz / -333.567; // pxy -> cpa
  mS[3][3]         = 1.;
  mS[4][4]         = 1.;
  HepSymMatrix mVy = err.similarity( mS );
  // std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
  return mVy;
}