00001 #include "TrackUtil/MdcTrackUtil.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "GaudiKernel/Bootstrap.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "EventModel/EventHeader.h"
00007 #include "RawEvent/RawDataUtil.h"
00008 #include "CLHEP/Units/PhysicalConstants.h"
00009
00010 #include "Identifier/MdcID.h"
00011 #include "MdcRawEvent/MdcDigi.h"
00012
00013 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00014
00015 typedef HepGeom::Point3D<double> HepPoint3D;
00016 #endif
00017
00018
00019 #include "MdcRecEvent/RecMdcTrack.h"
00020 #include "MdcRecEvent/RecMdcHit.h"
00021
00022
00023 #include "MdcGeomSvc/MdcGeomSvc.h"
00024
00025 using namespace std;
00026 using namespace Event;
00027
00028 MdcTrackUtil* MdcTrackUtil::_myself = 0;
00029
00030
00031 MdcTrackUtil* MdcTrackUtil::instance() {
00032 if( 0 == _myself ) {
00033 _myself = new MdcTrackUtil();
00034 }
00035 return _myself;
00036 }
00037
00038
00039 MdcTrackUtil::MdcTrackUtil(){
00040
00041 IService* svc;
00042 Gaudi::svcLocator()->getService("MagneticFieldSvc",svc);
00043 m_pIMF = dynamic_cast<IMagneticFieldSvc*> (svc);
00044 if(! m_pIMF){
00045 std::cout<<" ERROR Unable to open Magnetic field service "<<std::endl;
00046 }
00047
00048 double gaussToTesla = 1000.;
00049 Bz = m_pIMF->getReferField()*gaussToTesla;
00050
00051
00052 Gaudi::svcLocator()->getService("MdcGeomSvc",svc);
00053 m_mdcGeomSvc= dynamic_cast<MdcGeomSvc*> (svc);
00054 if(! m_mdcGeomSvc){
00055 std::cout<<" FATAL Could not load MdcGeomSvc! "<<std::endl;
00056 }
00057 }
00058
00059
00060 int MdcTrackUtil::nLayerTrackPassed(const HepVector helix){
00061 double helixParam[5];
00062 for(int i=0; i<5; ++i) helixParam[i] = helix[i];
00063
00064 return nLayerTrackPassed(helixParam);
00065 }
00066
00067
00068 int MdcTrackUtil::nLayerTrackPassed(const double helix[5]){
00069 int nLayer = 0;
00070
00071 for(unsigned iLayer=0; iLayer<43; iLayer++){
00072
00073
00074 double rMidLayer = m_mdcGeomSvc->Layer(iLayer)->Radius();
00075 double flightLength = rMidLayer;
00076
00077 HepPoint3D pivot(0.,0.,0.);
00078 double dz = helix[3];
00079 double c = CLHEP::c_light * 100.;
00080 double alpha = 1/(c * Bz);
00081 double kappa = helix[2];
00082 double rc = (-1.)*alpha/kappa;
00083
00084 double tanl = helix[4];
00085 double phi0 = helix[1];
00086 double phi = flightLength/rc + phi0;
00087 double z = pivot.z() + dz - (alpha/kappa) * tanl * phi;
00088
00089 double layerHalfLength = m_mdcGeomSvc->Layer(iLayer)->Length()/2.;
00090
00091
00092
00093 if (fabs(z) < fabs(layerHalfLength)) ++nLayer;
00094 }
00095
00096 return nLayer;
00097 }
00098
00099
00100 HepVector MdcTrackUtil::patRecPar2BesPar(const HepVector& helixPar){
00101 HepVector helix(5,0);
00102 double d0 = -helixPar[0];
00103 double phi0 = helixPar[1]+ CLHEP::halfpi;
00104 double omega = Bz*helixPar[2]/-333.567;
00105 double z0 = helixPar[3];
00106 double tanl = helixPar[4];
00107 helix[0] = d0;
00108 helix[1] = phi0;
00109 helix[2] = omega;
00110 helix[3] = z0;
00111 helix[4] = tanl;
00112
00113 return helix;
00114 }
00115
00116
00117 HepSymMatrix MdcTrackUtil::patRecErr2BesErr(const HepSymMatrix& err){
00118
00119
00120
00121
00122 HepSymMatrix mS(err.num_row(),0);
00123 mS[0][0]=-1.;
00124 mS[1][1]=1.;
00125 mS[2][2]=Bz/-333.567;
00126 mS[3][3]=1.;
00127 mS[4][4]=1.;
00128 HepSymMatrix mVy= err.similarity(mS);
00129
00130 return mVy;
00131 }
00132