BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
KalFitTrack Class Reference

Description of a track class (<- Helix.cc) More...

#include <KalFitTrack.h>

+ Inheritance diagram for KalFitTrack:

Public Member Functions

 KalFitTrack (const HepPoint3D &pivot, const CLHEP::HepVector &a, const CLHEP::HepSymMatrix &Ea, unsigned int m, double chiSq, unsigned int nhits)
 constructor
 
 ~KalFitTrack (void)
 destructor
 
double intersect_cylinder (double r) const
 Intersection with different geometry.
 
double intersect_yz_plane (const HepTransform3D &plane, double x) const
 
double intersect_zx_plane (const HepTransform3D &plane, double y) const
 
double intersect_xy_plane (double z) const
 
void msgasmdc (double path, int index)
 Calculate multiple scattering angle.
 
void ms (double path, const KalFitMaterial &m, int index)
 
void eloss (double path, const KalFitMaterial &m, int index)
 Calculate total energy lost in material.
 
double smoother_Mdc (KalFitHelixSeg &seg, CLHEP::Hep3Vector &meas, int &flg, int csmflag)
 Kalman smoother for Mdc.
 
double smoother_Mdc_csmalign (KalFitHelixSeg &seg, CLHEP::Hep3Vector &meas, int &flg, int csmflag)
 
double smoother_Mdc (KalFitHitMdc &HitMdc, CLHEP::Hep3Vector &meas, KalFitHelixSeg &seg, double &dchi2, int csmflag)
 
double update_hits (KalFitHitMdc &HitMdc, int inext, CLHEP::Hep3Vector &meas, int way, double &dchi2, double &dtrack, double &dtracknew, double &dtdc, int csmflag)
 Include the Mdc wire hits.
 
double update_hits (KalFitHelixSeg &HelixSeg, int inext, CLHEP::Hep3Vector &meas, int way, double &dchi2, int csmflag)
 
double update_hits_csmalign (KalFitHelixSeg &HelixSeg, int inext, CLHEP::Hep3Vector &meas, int way, double &dchi2, int csmflag)
 
double chi2_next (Helix &H, KalFitHitMdc &HitMdc, int csmflag)
 
double chi2_next (Helix &H, KalFitHitMdc &HitMdc)
 
void update_last (void)
 Record the current parameters as ..._last information :
 
void update_forMdc (void)
 
void point_last (const HepPoint3D &point)
 set and give out the last point of the track
 
const HepPoint3D & point_last (void)
 
const HepPoint3D & pivot_last (void) const
 returns helix parameters
 
const CLHEP::HepVector & a_last (void) const
 
const CLHEP::HepSymMatrix & Ea_last (void) const
 
const HepPoint3D & pivot_forMdc (void) const
 
const CLHEP::HepVector & a_forMdc (void) const
 
const CLHEP::HepSymMatrix & Ea_forMdc (void) const
 
void path_add (double path)
 Update the path length estimation.
 
void addPathSM (double path)
 
double getPathSM (void)
 
void addTofSM (double time)
 
double getTofSM (void)
 
void fiTerm (double fi)
 
double getFiTerm (void)
 
void tof (double path)
 Update the tof estimation.
 
double filter (double v_m, const CLHEP::HepVector &m_H, double v_d, double m_V)
 
double cor_tanldep (double *p, double er)
 Correct the error according the current tanl value :
 
void update_bit (int i)
 
int insist (void) const
 Extractor :
 
int type (void) const
 
int trasan_id (void) const
 
double r0 (void) const
 
double mass (void) const
 
double chiSq (void) const
 
double chiSq_back (void) const
 
int ndf_back (void) const
 
double pathip (void) const
 
double path_rd (void) const
 
double path_ab (void) const
 
double * pathl (void)
 
CLHEP::Hep3Vector * mom (void)
 
double tof (void) const
 
double tof_kaon (void) const
 
double tof_proton (void) const
 
double p_kaon (void) const
 
double p_proton (void) const
 
double dchi2_max (void) const
 
double r_max (void) const
 
unsigned int nchits (void) const
 
unsigned int nster (void) const
 
unsigned int ncath (void) const
 
int pat1 (void) const
 
int pat2 (void) const
 
int nhit_r (void) const
 
int nhit_z (void) const
 
void type (int t)
 Reinitialize (modificator)
 
void trasan_id (int t)
 
void insist (int t)
 
void pathip (double pl)
 
void p_kaon (double pl)
 
void p_proton (double pl)
 
void chiSq (double c)
 
void chiSq_back (double c)
 
void ndf_back (int n)
 
void nchits (int n)
 
void nster (int n)
 
void add_nhit_r (void)
 
void add_nhit_z (void)
 
double PathL (int layer)
 Function to calculate the path length in the layer.
 
void appendHitsMdc (KalFitHitMdc h)
 Functions for Mdc hits list.
 
void HitsMdc (vector< KalFitHitMdc > &lh)
 
vector< KalFitHitMdc > & HitsMdc (void)
 
KalFitHitMdc & HitMdc (int i)
 
void appendHelixSegs (KalFitHelixSeg s)
 
void HelixSegs (vector< KalFitHelixSeg > &vs)
 
vector< KalFitHelixSeg > & HelixSegs (void)
 
KalFitHelixSeg & HelixSeg (int i)
 
void order_wirhit (int index)
 
void order_hits (void)
 
void number_wirhit (void)
 
const HepPoint3D & pivot_numf (const HepPoint3D &newPivot)
 Sets pivot position in a given mag field.
 
const HepPoint3D & pivot_numf (const HepPoint3D &newPivot, double &pathl)
 
double radius_numf (void) const
 Estimation of the radius in a given mag field.
 
double getSigma (int layerId, double driftDist) const
 
double getSigma (KalFitHitMdc &hitmdc, double tanlam, int lr, double dist) const
 
double getDriftDist (KalFitHitMdc &hitmdc, double drifttime, int lr) const
 
double getDriftTime (KalFitHitMdc &hitmdc, double toftime) const
 
double getT0 (void) const
 
HepSymMatrix getInitMatrix (void) const
 
double getDigi () const
 
void chgmass (int i)
 
int nLayerUsed ()
 
void resetLayerUsed ()
 
void useLayer (int iLay)
 
- Public Member Functions inherited from KalmanFit::Helix
 Helix (const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
 Constructor with pivot, helix parameter a, and its error matrix.
 
 Helix (const HepPoint3D &pivot, const HepVector &a)
 Constructor without error matrix.
 
 Helix (const HepPoint3D &position, const Hep3Vector &momentum, double charge)
 Constructor with position, momentum, and charge.
 
virtual ~Helix ()
 Destructor.
 
const HepPoint3D & center (void) const
 returns position of helix center(z = 0.);
 
const HepPoint3D & pivot (void) const
 returns pivot position.
 
double radius (void) const
 returns radious of helix.
 
HepPoint3D x (double dPhi=0.) const
 returns position after rotating angle dPhi in phi direction.
 
double * x (double dPhi, double p[3]) const
 
HepPoint3D x (double dPhi, HepSymMatrix &Ex) const
 returns position and convariance matrix(Ex) after rotation.
 
Hep3Vector direction (double dPhi=0.) const
 returns direction vector after rotating angle dPhi in phi direction.
 
Hep3Vector momentum (double dPhi=0.) const
 returns momentum vector after rotating angle dPhi in phi direction.
 
Hep3Vector momentum (double dPhi, HepSymMatrix &Em) const
 returns momentum vector after rotating angle dPhi in phi direction.
 
HepLorentzVector momentum (double dPhi, double mass) const
 returns 4momentum vector after rotating angle dPhi in phi direction.
 
HepLorentzVector momentum (double dPhi, double mass, HepSymMatrix &Em) const
 returns 4momentum vector after rotating angle dPhi in phi direction.
 
HepLorentzVector momentum (double dPhi, double mass, HepPoint3D &x, HepSymMatrix &Emx) const
 returns 4momentum vector after rotating angle dPhi in phi direction.
 
double dr (void) const
 returns an element of parameters.
 
double phi0 (void) const
 
double kappa (void) const
 
double dz (void) const
 
double tanl (void) const
 
double curv (void) const
 
double sinPhi0 (void) const
 
double cosPhi0 (void) const
 
const HepVector & a (void) const
 returns helix parameters.
 
const HepSymMatrix & Ea (void) const
 returns error matrix.
 
double approach (KalFitHitMdc &hit, bool doSagCorrection) const
 
double approach (HepPoint3D pfwd, HepPoint3D pbwd, bool doSagCorrection) const
 
const HepVector & a (const HepVector &newA)
 sets helix parameters.
 
const HepSymMatrix & Ea (const HepSymMatrix &newdA)
 sets helix paramters and error matrix.
 
const HepPoint3D & pivot (const HepPoint3D &newPivot)
 sets pivot position.
 
void set (const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
 sets helix pivot position, parameters, and error matrix.
 
void ignoreErrorMatrix (void)
 unsets error matrix. Error calculations will be ignored after this function call until an error matrix be set again. 0 matrix will be return as a return value for error matrix when you call functions which returns an error matrix.
 
double bFieldZ (double)
 sets/returns z componet of the magnetic field.
 
double bFieldZ (void) const
 
double alpha (void) const
 
Helix & operator= (const Helix &)
 Copy operator.
 
HepMatrix delApDelA (const HepVector &ap) const
 
HepMatrix delXDelA (double phi) const
 
HepMatrix delMDelA (double phi) const
 
HepMatrix del4MDelA (double phi, double mass) const
 
HepMatrix del4MXDelA (double phi, double mass) const
 

Static Public Member Functions

static void setT0 (double t0)
 
static void setInitMatrix (HepSymMatrix m)
 
static void setMdcCalibFunSvc (const MdcCalibFunSvc *calibsvc)
 
static void setMagneticFieldSvc (IMagneticFieldSvc *)
 
static void setIMdcGeomSvc (IMdcGeomSvc *igeomsvc)
 
static void setMdcDigiCol (MdcDigiCol *digicol)
 
static int nmass (void)
 
static double mass (int i)
 
static void LR (int x)
 
static int lead (void)
 Magnetic field map.
 
static void lead (int i)
 
static int back (void)
 
static void back (int i)
 
static int resol (void)
 
static void resol (int i)
 
static int numf (void)
 
static void numf (int i)
 

Static Public Attributes

static double mdcGasRadlen_ = 0.
 
static int tprop_ = 1
 for signal propagation correction
 
static int debug_ = 0
 for debug

 
static double chi2_hitf_ = 1000
 Cut chi2 for each hit.
 
static double chi2_hits_ = 1000
 
static int numf_ = 0
 Flag for treatment of non-uniform mag field.
 
static int inner_steps_ = 3
 
static int outer_steps_ = 3
 
static double dchi2cutf_anal [43] = {0.}
 
static double dchi2cuts_anal [43] = {0.}
 
static double dchi2cutf_calib [43] = {0.}
 
static double dchi2cuts_calib [43] = {0.}
 
static int numfcor_ = 1
 NUMF treatment improved.
 
static double Bznom_ = 10
 
static int steplev_ = 0
 
static int Tof_correc_ = 1
 Flag for TOF correction.
 
static int strag_ = 1
 Flag to take account of energy loss straggling :
 
static double factor_strag_ = 0.4
 factor of energy loss straggling for electron
 
static int nmdc_hit2_ = 500
 Cut chi2 for each hit.
 
static double chi2mdc_hit2_
 
static int tofall_ = 1
 
static int resolflag_ = 0
 wire resoltion flag
 
static int LR_ = 1
 Use L/R decision from MdcRecHit information :
 
static int drifttime_choice_ = 0
 the drifttime choice
 
- Static Public Attributes inherited from KalmanFit::Helix
static const double ConstantAlpha = 333.564095
 Constant alpha for uniform field.
 

Detailed Description

Description of a track class (<- Helix.cc)

Definition at line 34 of file KalFitTrack.h.

Constructor & Destructor Documentation

◆ KalFitTrack()

KalFitTrack::KalFitTrack ( const HepPoint3D &  pivot,
const CLHEP::HepVector &  a,
const CLHEP::HepSymMatrix &  Ea,
unsigned int  m,
double  chiSq,
unsigned int  nhits 
)

constructor

Definition at line 81 of file KalFitTrack.cxx.

85: Helix(pivot, a, Ea), type_(0), insist_(0), chiSq_(0),
86 nchits_(0), nster_(0), ncath_(0),
87 ndf_back_(0), chiSq_back_(0), pathip_(0),
88 path_rd_(0), path_ab_(0), tof_(0), dchi2_max_(0), r_max_(0),
89 tof_kaon_(0), tof_proton_(0), pat1_(0), pat2_(0), layer_prec_(0),
90 trasan_id_(0), nhit_r_(0), nhit_z_(0),pathSM_(0),tof2_(0)
91{
92 memset(PathL_, 0, sizeof(PathL_));
93 l_mass_ = m;
94 if (m < NMASS) mass_ = MASS[m];
95 else mass_ = MASS[2];
96 r0_ = fabs(center().perp() - fabs(radius()));
97 //bFieldZ(Bznom_);
98 Bznom_=bFieldZ(); // 2012-09-13 wangll
99 update_last();
100}
KalFitTrack::Bznom_
static double Bznom_
Definition: KalFitTrack.h:299
KalFitTrack::update_last
void update_last(void)
Record the current parameters as ..._last information :
Definition: KalFitTrack.cxx:109
KalmanFit::Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/helix/Helix.h:194
KalmanFit::Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/helix/Helix.h:260
KalmanFit::Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/helix/Helix.h:254
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/helix/Helix.h:200

◆ ~KalFitTrack()

KalFitTrack::~KalFitTrack ( void  )

destructor

Definition at line 103 of file KalFitTrack.cxx.

104{
105 // delete all objects
106
107}

Member Function Documentation

◆ a_forMdc()

const CLHEP::HepVector & KalFitTrack::a_forMdc ( void  ) const
inline

Definition at line 151 of file KalFitTrack.h.

151{ return a_forMdc_; }

◆ a_last()

const CLHEP::HepVector & KalFitTrack::a_last ( void  ) const
inline

Definition at line 147 of file KalFitTrack.h.

147{ return a_last_; }

◆ add_nhit_r()

void KalFitTrack::add_nhit_r ( void  )
inline

Definition at line 220 of file KalFitTrack.h.

220{ nhit_r_++;}

◆ add_nhit_z()

void KalFitTrack::add_nhit_z ( void  )
inline

Definition at line 221 of file KalFitTrack.h.

221{ nhit_z_++;}

◆ addPathSM()

void KalFitTrack::addPathSM ( double  path)

Definition at line 1296 of file KalFitTrack.cxx.

1296 {
1297 pathSM_ += path;
1298}

Referenced by KalFitAlg::smoother_anal().

◆ addTofSM()

void KalFitTrack::addTofSM ( double  time)

Definition at line 1301 of file KalFitTrack.cxx.

1301 {
1302 tof2_ += time;
1303}

Referenced by KalFitAlg::smoother_anal().

◆ appendHelixSegs()

void KalFitTrack::appendHelixSegs ( KalFitHelixSeg  s)

◆ appendHitsMdc()

void KalFitTrack::appendHitsMdc ( KalFitHitMdc  h)

Functions for Mdc hits list.

Definition at line 1837 of file KalFitTrack.cxx.

1837{ HitsMdc_.push_back(h);}

Referenced by KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), and KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew().

◆ back() [1/2]

void KalFitTrack::back ( int  i)
static

Definition at line 1829 of file KalFitTrack.cxx.

1829{ back_ = i;}

◆ back() [2/2]

int KalFitTrack::back ( void  )
static

Definition at line 1830 of file KalFitTrack.cxx.

1830{ return back_;}

Referenced by KalFitAlg::initialize().

◆ chgmass()

void KalFitTrack::chgmass ( int  i)

Definition at line 1822 of file KalFitTrack.cxx.

1822 {
1823 mass_=MASS[i];
1824 l_mass_=i;
1825}

Referenced by KalFitAlg::complete_track().

◆ chi2_next() [1/2]

double KalFitTrack::chi2_next ( Helix &  H,
KalFitHitMdc &  HitMdc 
)

Definition at line 2805 of file KalFitTrack.cxx.

2805 {
2806
2807 double lr = HitMdc.LR();
2808 const KalFitWire& Wire = HitMdc.wire();
2809 int wire_ID = Wire.geoID();
2810 int layerid = HitMdc.wire().layer().layerId();
2811 double entrangle = HitMdc.rechitptr()->getEntra();
2812
2813 HepPoint3D fwd(Wire.fwd());
2814 HepPoint3D bck(Wire.bck());
2815 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
2816 Helix work = H;
2817 work.ignoreErrorMatrix();
2818 work.pivot((fwd + bck) * .5);
2819 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
2820 H.pivot(x0kal);
2821
2822 Hep3Vector meas = H.momentum(0).cross(wire).unit();
2823
2824 if (wire_ID<0 || wire_ID>6796){ //bes
2825 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
2826 << std::endl;
2827 return DBL_MAX;
2828 }
2829
2830 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
2831 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
2832 double tofest(0);
2833 double phi = fmod(phi0() + M_PI4, M_PI2);
2834 double csf0 = cos(phi);
2835 double snf0 = (1. - csf0) * (1. + csf0);
2836 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
2837 if(phi > M_PI) snf0 = - snf0;
2838
2839 if (Tof_correc_) {
2840 Hep3Vector ip(0, 0, 0);
2841 Helix work = *(Helix*)this;
2842 work.ignoreErrorMatrix();
2843 work.pivot(ip);
2844 double phi_ip = work.phi0();
2845 if (fabs(phi - phi_ip) > M_PI) {
2846 if (phi > phi_ip) phi -= 2 * M_PI;
2847 else phi_ip -= 2 * M_PI;
2848 }
2849 double t = tanl();
2850 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
2851 double pmag( sqrt( 1.0 + t*t ) / kappa());
2852 double mass_over_p( mass_ / pmag );
2853 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2854 tofest = l / ( 29.9792458 * beta );
2855 // if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
2856 }
2857
2858 const HepSymMatrix& ea = H.Ea();
2859 const HepVector& v_a = H.a();
2860 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
2861
2862 HepVector v_H(5, 0);
2863 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
2864 v_H[3] = -meas.z();
2865 HepMatrix v_HT = v_H.T();
2866
2867 double estim = (v_HT * v_a)[0];
2868 HepVector ea_v_H = ea * v_H;
2869 HepMatrix ea_v_HT = (ea_v_H).T();
2870 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
2871
2872 HepSymMatrix eaNewL(5), eaNewR(5);
2873 HepVector aNewL(5), aNewR(5);
2874
2875 //double time = HitMdc.tdc();
2876 //if (Tof_correc_)
2877 // time = time - tofest;
2878 double drifttime = getDriftTime(HitMdc , tofest);
2879 double ddl = getDriftDist(HitMdc, drifttime , 0 );
2880 double ddr = getDriftDist(HitMdc, drifttime , 1 );
2881 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
2882 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
2883
2884 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
2885 double er_dmeas2[2] = {0. , 0.};
2886 if(resolflag_ == 1) {
2887 er_dmeas2[0] = 0.1*erddl;
2888 er_dmeas2[1] = 0.1*erddr;
2889 }
2890 else if(resolflag_ == 0) {
2891 // int layid = HitMdc.wire().layer().layerId();
2892 // double sigma = getSigma(layid, dd);
2893 // er_dmeas2[0] = er_dmeas2[1] = sigma;
2894 }
2895
2896 if ((LR_==0 && lr != 1.0) ||
2897 (LR_==1 && lr == -1.0)) {
2898
2899 double er_dmeasL, dmeasL;
2900 if(Tof_correc_) {
2901 dmeasL = (-1.0)*fabs(dmeas2[0]);
2902 er_dmeasL = er_dmeas2[0];
2903 } else {
2904 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
2905 er_dmeasL = HitMdc.erdist()[0];
2906 }
2907
2908 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
2909 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
2910 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
2911 if(0. == RkL) RkL = 1.e-4;
2912
2913 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
2914 aNewL = v_a + diffL;
2915 double sigL = dmeasL -(v_HT * aNewL)[0];
2916 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
2917 }
2918
2919 if ((LR_==0 && lr != -1.0) ||
2920 (LR_==1 && lr == 1.0)) {
2921
2922 double er_dmeasR, dmeasR;
2923 if(Tof_correc_) {
2924 dmeasR = dmeas2[1];
2925 er_dmeasR = er_dmeas2[1];
2926 } else {
2927 dmeasR = fabs(HitMdc.dist()[1]);
2928 er_dmeasR = HitMdc.erdist()[1];
2929 }
2930
2931 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
2932 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
2933 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
2934 if(0. == RkR) RkR = 1.e-4;
2935
2936 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
2937 aNewR = v_a + diffR;
2938 double sigR = dmeasR -(v_HT * aNewR)[0];
2939 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
2940 }
2941
2942 if (dchi2R < dchi2L){
2943 H.a(aNewR); H.Ea(eaNewR);
2944 } else {
2945 H.a(aNewL); H.Ea(eaNewL);
2946 }
2947 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
2948}
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
momentum
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
Definition: BesBdkRc/BesBdkRc-00-00-02/src/fortran/HepEvt.h:23
DBL_MAX
#define DBL_MAX
Definition: KalFitTrack.h:4
M_PI
#define M_PI
Definition: TConstant.h:4
t
TTree * t
Definition: binning.cxx:23
KalFitHitMdc::LR
int LR(void) const
Definition: KalFitHitMdc.h:29
KalFitHitMdc::rechitptr
RecMdcHit * rechitptr(void)
Definition: KalFitHitMdc.h:36
KalFitHitMdc::erdist
const double * erdist(void) const
Definition: KalFitHitMdc.h:34
KalFitHitMdc::dist
const double * dist(void) const
Definition: KalFitHitMdc.h:33
KalFitHitMdc::wire
const KalFitWire & wire(void) const
Definition: KalFitHitMdc.h:35
KalFitLayer_Mdc::layerId
const int layerId(void) const
returns layer ID
Definition: KalFitLayer_Mdc.h:95
KalFitTrack::getDriftTime
double getDriftTime(KalFitHitMdc &hitmdc, double toftime) const
Definition: KalFitTrack2.cxx:74
KalFitTrack::getSigma
double getSigma(int layerId, double driftDist) const
Definition: KalFitTrack.cxx:3096
KalFitTrack::resolflag_
static int resolflag_
wire resoltion flag
Definition: KalFitTrack.h:316
KalFitTrack::Tof_correc_
static int Tof_correc_
Flag for TOF correction.
Definition: KalFitTrack.h:305
KalFitTrack::getDriftDist
double getDriftDist(KalFitHitMdc &hitmdc, double drifttime, int lr) const
Definition: KalFitTrack2.cxx:194
KalFitTrack::LR_
static int LR_
Use L/R decision from MdcRecHit information :
Definition: KalFitTrack.h:319
KalFitTrack::HitMdc
KalFitHitMdc & HitMdc(int i)
Definition: KalFitTrack.h:230
KalFitWire
Description of a Wire class.
Definition: KalFitWire.h:46
KalFitWire::geoID
unsigned int geoID(void) const
Definition: KalFitWire.h:74
KalFitWire::bck
HepPoint3D bck(void) const
Definition: KalFitWire.h:93
KalFitWire::fwd
HepPoint3D fwd(void) const
Geometry :
Definition: KalFitWire.h:92
KalFitWire::layer
const KalFitLayer_Mdc & layer(void) const
Definition: KalFitWire.h:70
KalmanFit::Helix::ignoreErrorMatrix
void ignoreErrorMatrix(void)
unsets error matrix. Error calculations will be ignored after this function call until an error matri...
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/src/helix/Helix.cxx:722
KalmanFit::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/src/helix/Helix.cxx:131
RecMdcHit::getEntra
const double getEntra(void) const
Definition: RecMdcHit.h:54
H
IMPLICIT REAL *A H
Definition: myXsection.h:1

◆ chi2_next() [2/2]

double KalFitTrack::chi2_next ( Helix &  H,
KalFitHitMdc &  HitMdc,
int  csmflag 
)

Definition at line 2950 of file KalFitTrack.cxx.

2950 {
2951
2952 double lr = HitMdc.LR();
2953 const KalFitWire& Wire = HitMdc.wire();
2954 int wire_ID = Wire.geoID();
2955 int layerid = HitMdc.wire().layer().layerId();
2956 double entrangle = HitMdc.rechitptr()->getEntra();
2957
2958 HepPoint3D fwd(Wire.fwd());
2959 HepPoint3D bck(Wire.bck());
2960 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
2961 Helix work = H;
2962 work.ignoreErrorMatrix();
2963 work.pivot((fwd + bck) * .5);
2964 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
2965 H.pivot(x0kal);
2966
2967 Hep3Vector meas = H.momentum(0).cross(wire).unit();
2968
2969 if (wire_ID<0 || wire_ID>6796){ //bes
2970 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
2971 << std::endl;
2972 return DBL_MAX;
2973 }
2974
2975 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
2976 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
2977 double tofest(0);
2978 double phi = fmod(phi0() + M_PI4, M_PI2);
2979 double csf0 = cos(phi);
2980 double snf0 = (1. - csf0) * (1. + csf0);
2981 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
2982 if(phi > M_PI) snf0 = - snf0;
2983
2984 if (Tof_correc_) {
2985 Hep3Vector ip(0, 0, 0);
2986 Helix work = *(Helix*)this;
2987 work.ignoreErrorMatrix();
2988 work.pivot(ip);
2989 double phi_ip = work.phi0();
2990 if (fabs(phi - phi_ip) > M_PI) {
2991 if (phi > phi_ip) phi -= 2 * M_PI;
2992 else phi_ip -= 2 * M_PI;
2993 }
2994 double t = tanl();
2995 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
2996 double pmag( sqrt( 1.0 + t*t ) / kappa());
2997 double mass_over_p( mass_ / pmag );
2998 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2999 tofest = l / ( 29.9792458 * beta );
3000 if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
3001 }
3002
3003 const HepSymMatrix& ea = H.Ea();
3004 const HepVector& v_a = H.a();
3005 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
3006
3007 HepVector v_H(5, 0);
3008 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
3009 v_H[3] = -meas.z();
3010 HepMatrix v_HT = v_H.T();
3011
3012 double estim = (v_HT * v_a)[0];
3013 HepVector ea_v_H = ea * v_H;
3014 HepMatrix ea_v_HT = (ea_v_H).T();
3015 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
3016
3017 HepSymMatrix eaNewL(5), eaNewR(5);
3018 HepVector aNewL(5), aNewR(5);
3019
3020 //double time = HitMdc.tdc();
3021 //if (Tof_correc_)
3022 // time = time - tofest;
3023 double drifttime = getDriftTime(HitMdc , tofest);
3024 double ddl = getDriftDist(HitMdc, drifttime , 0 );
3025 double ddr = getDriftDist(HitMdc, drifttime , 1 );
3026 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
3027 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
3028
3029 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
3030 double er_dmeas2[2] = {0. , 0.};
3031 if(resolflag_ == 1) {
3032 er_dmeas2[0] = 0.1*erddl;
3033 er_dmeas2[1] = 0.1*erddr;
3034 }
3035 else if(resolflag_ == 0) {
3036 // int layid = HitMdc.wire().layer().layerId();
3037 // double sigma = getSigma(layid, dd);
3038 // er_dmeas2[0] = er_dmeas2[1] = sigma;
3039 }
3040
3041 if ((LR_==0 && lr != 1.0) ||
3042 (LR_==1 && lr == -1.0)) {
3043
3044 double er_dmeasL, dmeasL;
3045 if(Tof_correc_) {
3046 dmeasL = (-1.0)*fabs(dmeas2[0]);
3047 er_dmeasL = er_dmeas2[0];
3048 } else {
3049 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
3050 er_dmeasL = HitMdc.erdist()[0];
3051 }
3052
3053 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
3054 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
3055 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
3056 if(0. == RkL) RkL = 1.e-4;
3057
3058 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
3059 aNewL = v_a + diffL;
3060 double sigL = dmeasL -(v_HT * aNewL)[0];
3061 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
3062 }
3063
3064 if ((LR_==0 && lr != -1.0) ||
3065 (LR_==1 && lr == 1.0)) {
3066
3067 double er_dmeasR, dmeasR;
3068 if(Tof_correc_) {
3069 dmeasR = dmeas2[1];
3070 er_dmeasR = er_dmeas2[1];
3071 } else {
3072 dmeasR = fabs(HitMdc.dist()[1]);
3073 er_dmeasR = HitMdc.erdist()[1];
3074 }
3075
3076 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
3077 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
3078 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
3079 if(0. == RkR) RkR = 1.e-4;
3080
3081 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
3082 aNewR = v_a + diffR;
3083 double sigR = dmeasR -(v_HT * aNewR)[0];
3084 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
3085 }
3086
3087 if (dchi2R < dchi2L){
3088 H.a(aNewR); H.Ea(eaNewR);
3089 } else {
3090 H.a(aNewL); H.Ea(eaNewL);
3091 }
3092 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
3093}
KalFitWire::y
double y(void) const
Definition: KalFitWire.h:101

Referenced by update_hits_csmalign().

◆ chiSq() [1/2]

void KalFitTrack::chiSq ( double  c)
inline

Definition at line 214 of file KalFitTrack.h.

214{ chiSq_ = c; }

◆ chiSq() [2/2]

double KalFitTrack::chiSq ( void  ) const
inline

Definition at line 183 of file KalFitTrack.h.

183{ return chiSq_; }

Referenced by KalFitAlg::fillTds(), KalFitAlg::fillTds_ip(), KalFitAlg::fillTds_lead(), filter(), KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), and KalFitAlg::start_seed().

◆ chiSq_back() [1/2]

void KalFitTrack::chiSq_back ( double  c)
inline

Definition at line 215 of file KalFitTrack.h.

215{ chiSq_back_ = c; }

◆ chiSq_back() [2/2]

double KalFitTrack::chiSq_back ( void  ) const
inline

Definition at line 184 of file KalFitTrack.h.

184{ return chiSq_back_; }

Referenced by KalFitAlg::fillTds_back(), and KalFitAlg::start_seed().

◆ cor_tanldep()

double KalFitTrack::cor_tanldep ( double *  p,
double  er 
)

Correct the error according the current tanl value :

◆ dchi2_max()

double KalFitTrack::dchi2_max ( void  ) const
inline

Definition at line 196 of file KalFitTrack.h.

196{ return dchi2_max_; }

◆ Ea_forMdc()

const CLHEP::HepSymMatrix & KalFitTrack::Ea_forMdc ( void  ) const
inline

Definition at line 152 of file KalFitTrack.h.

152{ return Ea_forMdc_; }

◆ Ea_last()

const CLHEP::HepSymMatrix & KalFitTrack::Ea_last ( void  ) const
inline

Definition at line 148 of file KalFitTrack.h.

148{ return Ea_last_; }

◆ eloss()

void KalFitTrack::eloss ( double  path,
const KalFitMaterial &  m,
int  index 
)

Calculate total energy lost in material.

Definition at line 276 of file KalFitTrack.cxx.

278{
279#ifdef YDEBUG
280 cout<<"eloss():ea before: "<<Ea()<<endl;
281#endif
282 HepVector v_a = a();
283 double v_a_2_2 = v_a[2] * v_a[2];
284 double v_a_4_2 = v_a[4] * v_a[4];
285 double pmag = 1 / fabs(v_a[2]) * sqrt(1 + v_a_4_2);
286 double psq = pmag * pmag;
287 double E = sqrt(mass_ * mass_ + psq);
288 double dE = material.dE(mass_, path, pmag);
289 //std::cout<<" eloss(): dE: "<<dE<<std::endl;//wangll
290
291 if (index > 0)
292 psq += dE * (dE + 2 * sqrt(mass_ * mass_ + psq));
293 else {
294 double dE_max = E - mass_;
295 if( dE<dE_max ) psq += dE * (dE - 2 * sqrt(mass_ * mass_ + psq));
296 else psq=-1.0;
297 }
298
299 if (tofall_ && index < 0){
300 // Kaon case :
301 if (p_kaon_ > 0){
302 double psq_kaon = p_kaon_ * p_kaon_;
303 double dE_kaon = material.dE(MASS[3], path, p_kaon_);
304 psq_kaon += dE_kaon * (dE_kaon -
305 2 * sqrt(MASS[3] * MASS[3] + psq_kaon));
306 if (psq_kaon < 0) psq_kaon = 0;
307 p_kaon_ = sqrt(psq_kaon);
308 }
309 // Proton case :
310 if (p_proton_ > 0){
311 double psq_proton = p_proton_ * p_proton_;
312 double dE_proton = material.dE(MASS[4], path, p_proton_);
313 psq_proton += dE_proton * (dE_proton -
314 2 * sqrt(MASS[4] * MASS[4] + psq_proton));
315 if (psq_proton < 0) psq_proton = 0;
316 p_proton_ = sqrt(psq_proton);
317 }
318 }
319
320 double dpt;
321 //cout<<"eloss(): psq = "<<psq<<endl;//wangll
322 if(psq < 0) dpt = 9999;
323 else dpt = v_a[2] * pmag / sqrt(psq);
324 //cout<<"eloss():k: "<<v_a[2]<<" k' "<<dpt<<endl;//wangll
325#ifdef YDEBUG
326 cout<<"eloss():k: "<<v_a[2]<<" k' "<<dpt<<endl;
327#endif
328 // attempt to take account of energy loss for error matrix
329
330 HepSymMatrix ea = Ea();
331 double r_E_prim = (E + dE)/E;
332
333 // 1/ Straggling in the energy loss :
334 if (strag_){
335 double del_E(0);
336 if(l_mass_==0) {
337 del_E = dE*factor_strag_;
338 } else {
339 del_E = material.del_E(mass_, path, pmag);
340 }
341
342 ea[2][2] += (v_a[2] * v_a[2]) * E * E * del_E * del_E / (psq*psq);
343 }
344
345 // Effect of the change of curvature (just variables change):
346 double dpt2 = dpt*dpt;
347 double A = dpt*dpt2/(v_a_2_2*v_a[2])*r_E_prim;
348 double B = v_a[4]/(1+v_a_4_2)*
349 dpt*(1-dpt2/v_a_2_2*r_E_prim);
350
351 double ea_2_0 = A*ea[2][0] + B*ea[4][0];
352 double ea_2_1 = A*ea[2][1] + B*ea[4][1];
353 double ea_2_2 = A*A*ea[2][2] + 2*A*B*ea[2][4] + B*B*ea[4][4];
354 double ea_2_3 = A*ea[2][3] + B*ea[4][3];
355 double ea_2_4 = A*ea[2][4] + B*ea[4][4];
356
357 v_a[2] = dpt;
358 a(v_a);
359
360 ea[2][0] = ea_2_0;
361 //std::cout<<"ea[2][0] is "<<ea[2][0]<<" ea(3,1) is "<<ea(3,1)<<std::endl;
362 ea[2][1] = ea_2_1;
363 //std::cout<<"ea[2][2] is "<<ea[2][2]<<" ea(3,3) is "<<ea(3,3)<<std::endl;
364 ea[2][2] = ea_2_2;
365 ea[2][3] = ea_2_3;
366 ea[2][4] = ea_2_4;
367
368 Ea(ea);
369 //cout<<"eloss():dE: "<<dE<<endl;
370 //cout<<"eloss():A: "<<A<<" B: "<<B<<endl;
371 //cout<<"eloss():ea after: "<<Ea()<<endl;
372 r0_ = fabs(center().perp() - fabs(radius()));
373}
KalFitTrack::tofall_
static int tofall_
Definition: KalFitTrack.h:314
KalFitTrack::factor_strag_
static double factor_strag_
factor of energy loss straggling for electron
Definition: KalFitTrack.h:310
KalFitTrack::strag_
static int strag_
Flag to take account of energy loss straggling :
Definition: KalFitTrack.h:308

Referenced by KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), KalFitAlg::smoother_anal(), KalFitAlg::smoother_calib(), KalFitElement::updateTrack(), KalFitElement::updateTrack_alreadyfound(), and KalFitElement::updateTrack_rphi().

◆ filter()

double KalFitTrack::filter ( double  v_m,
const CLHEP::HepVector &  m_H,
double  v_d,
double  m_V 
)

Definition at line 1254 of file KalFitTrack.cxx.

1256{
1257 HepMatrix m_HT = m_H.T();
1258 HepPoint3D x0 = x(0);
1259 HepVector v_x(3);
1260 v_x[0] = x0.x();
1261 v_x[1] = x0.y();
1262 v_x[2] = x0.z();
1263 HepMatrix m_X = delXDelA(0);
1264 HepMatrix m_XT = m_X.T();
1265 HepMatrix m_C = m_X * Ea() * m_XT;
1266 //int err;
1267 HepVector m_K = 1 / (m_V + (m_HT * m_C * m_H)[0]) * m_H;
1268 HepVector v_a = a();
1269 HepSymMatrix ea = Ea();
1270 HepMatrix mXTmK = m_XT * m_K;
1271 double v_r = v_m - v_d - (m_HT * v_x)[0];
1272 v_a += ea * mXTmK * v_r;
1273 a(v_a);
1274 ea.assign(ea - ea * mXTmK * m_HT * m_X * ea);
1275 Ea(ea);
1276 // Record last hit included parameters :
1277 update_last();
1278 HepMatrix mCmK = m_C * m_K;
1279 v_x += mCmK * v_r;
1280 m_C -= mCmK * m_HT * m_C;
1281 v_r = v_m - v_d - (m_HT * v_x)[0];
1282 double m_R = m_V - (m_HT * m_C * m_H)[0];
1283 double chiSq = v_r * v_r / m_R;
1284 chiSq_ += chiSq;
1285 return chiSq;
1286}
m_C
********INTEGER modcns REAL m_C
Definition: PseuMar.h:13
KalFitTrack::chiSq
double chiSq(void) const
Definition: KalFitTrack.h:183

◆ fiTerm()

void KalFitTrack::fiTerm ( double  fi)

Definition at line 1306 of file KalFitTrack.cxx.

1306 {
1307 fiTerm_ = fi;
1308}

Referenced by KalFitAlg::smoother_anal().

◆ getDigi()

double KalFitTrack::getDigi ( ) const

◆ getDriftDist()

double KalFitTrack::getDriftDist ( KalFitHitMdc &  hitmdc,
double  drifttime,
int  lr 
) const

Definition at line 194 of file KalFitTrack2.cxx.

195{
196 int layerid = hitmdc.wire().layer().layerId();
197 int cellid = MdcID::wire(hitmdc.rechitptr()->getMdcId());
198 if(debug_ == 4 ){
199 std::cout<<"the cellid is .."<<cellid<<std::endl;
200 }
201 double entrangle = hitmdc.rechitptr()->getEntra();
202
203 //std::cout<<" entrangle: "<<entrangle<<std::endl;
204
205 return CalibFunSvc_->driftTimeToDist(drifttime, layerid, cellid, lr, entrangle);
206}
KalFitTrack::debug_
static int debug_
for debug
Definition: KalFitTrack.h:280
MdcCalibFunSvc::driftTimeToDist
double driftTimeToDist(double drifttime, int layid, int cellid, int lr, double entrance=0.0) const
Definition: MdcCalibFunSvc.cxx:149
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
RecMdcHit::getMdcId
const Identifier getMdcId(void) const
Definition: RecMdcHit.h:49

Referenced by chi2_next(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ getDriftTime()

double KalFitTrack::getDriftTime ( KalFitHitMdc &  hitmdc,
double  toftime 
) const

Definition at line 74 of file KalFitTrack2.cxx.

75{
76 const double vinner = 220.0e8; // cm/s
77 const double vouter = 240.0e8; // cm/s
78
79 int layerid = hitmdc.wire().layer().layerId();
80 double zhit = (hitmdc.rechitptr())->getZhit();
81 const KalFitWire* wire = &(hitmdc.wire());
82 HepPoint3D fPoint = wire->fwd();
83 HepPoint3D bPoint = wire->bck();
84
85 // unit is centimeter
86 double wireLen = (fPoint-bPoint).x()*(fPoint-bPoint).x()+(fPoint-bPoint).y()*(fPoint-bPoint).y()+(fPoint-bPoint).z()*(fPoint-bPoint).z();
87 wireLen = sqrt(wireLen);
88 double wireZLen = fabs(fPoint.z() - bPoint.z());
89 double tp = 0.;
90 double vp = 0.;
91 // west readout
92 if(0==layerid%2){
93 // inner chamber
94 if(layerid<8){
95 vp = vinner;
96 }
97 else {
98 vp = vouter;
99 }
100 tp = wireLen*fabs(zhit - bPoint.z())/wireZLen/vp;
101 }
102
103 // east readout
104 if(1==layerid%2){
105 // inner chamber
106 if(layerid<8){
107 vp = vinner;
108 }
109 else {
110 vp = vouter;
111 }
112 tp = wireLen*fabs(zhit - fPoint.z())/wireZLen/vp;
113 }
114
115 // s to ns
116 tp = 1.0e9*tp;
117
118 if(0==tprop_) tp = 0.;
119
120 //std::cout<<"propogation time: "<<tp<<std::endl;
121
122 int wireid = hitmdc.wire().geoID();
123 double drifttime1(0.);
124 double drifttime2(0.);
125 double drifttime3(0.);
126
127 MdcDigiCol::iterator iter = mdcDigiCol_->begin();
128 for(; iter != mdcDigiCol_->end(); iter++ ) {
129 if((*iter)->identify() == (hitmdc.rechitptr())->getMdcId()) break;
130 }
131 //double t0 = get_T0();
132 // see the code of wulh in /Mdc/MdcCalibFunSvc/MdcCalibFunSvc/MdcCalibFunSvc.h,
133 // double getT0(int wireid) const { return m_t0[wireid]; }
134
135 //double getTimeWalk(int layid, double Q) const ;
136 double Q = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
137 double timewalk = CalibFunSvc_->getTimeWalk(layerid, Q);
138
139 if(debug_ == 4) {
140 std::cout<<"CalibFunSvc_->getTimeWalk, timewalk ="<<timewalk<<std::endl;
141 }
142
143 double timeoffset = CalibFunSvc_->getT0(wireid);
144 if(debug_ == 4) {
145 std::cout<<"CalibFunSvc_->getT0, timeoffset ="<<timeoffset<<std::endl;
146 }
147
148 double eventt0 = getT0();
149
150 if(debug_ == 4) {
151 std::cout<<"the Event T0 we get in the function getDriftTime(...) is "<<eventt0<<std::endl;
152 }
153
154 // this tdc value come from MdcRecHit assigned by zhangyao
155 double tdctime1 = hitmdc.tdc();
156 double tdctime2(0.);
157 double tdctime3(0.);
158
159 if(debug_ == 4) {
160 std::cout<<"tdctime1 be here is .."<<tdctime1<<std::endl;
161 }
162
163 // this tdc value come from MdcDigiCol time channel
164 // attention, if we use the iter like this: for(MdcDigiCol::iterator iter = mdcDigiCol_->begin(); iter != mdcDigiCol_->end(); iter++) it cannot pass the gmake , throw an error !
165 if(debug_ == 4) {
166 // std::cout<<"the size of the mdcDigiCol_ is "<<mdcDigiCol_.size()<<std::endl;
167 }
168 // MdcDigiCol::iterator iter = mdcDigiCol_->begin();
169 // for(; iter != mdcDigiCol_->end(); iter++ ) {
170 // if((*iter)->identify() == (hitmdc.rechitptr())->getMdcId()) break;
171 // }
172 if(debug_ == 4) {
173 std::cout<<"the time channel be here is .."<<(*iter)->getTimeChannel()<<std::endl;
174 }
175 tdctime2 = RawDataUtil::MdcTime((*iter)->getTimeChannel());
176 tdctime3 = hitmdc.rechitptr()->getTdc();
177 drifttime1 = tdctime1 - eventt0 - toftime - timewalk -timeoffset - tp;
178 drifttime2 = tdctime2 - eventt0 - toftime - timewalk -timeoffset - tp;
179 drifttime3 = tdctime3 - eventt0 - toftime - timewalk -timeoffset - tp;
180 if(debug_ == 4 ) {
181 std::cout<<"we now compare the three kind of tdc , the tdc get from timeChannel() is "<<tdctime2<<" the tdc get from KalFitHitMdc is "<<tdctime1<<" the tdc from MdcRecHit is "<<tdctime3<<" the drifttime1 is ..."<<drifttime1<<" the drifttime 2 is ..."<<drifttime2<<" the drifttime3 is ..."<<drifttime3<<std::endl;
182 }
183 //return drifttime3;
184 //return drifttime1;
185 if(drifttime_choice_ == 0)
186 return drifttime2;
187 if(drifttime_choice_ == 1)
188 // use the driftT caluculated by track-finding
189 return hitmdc.rechitptr()->getDriftT();
190}
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
KalFitHitMdc::tdc
double tdc(void) const
Definition: KalFitHitMdc.h:32
KalFitTrack::getT0
double getT0(void) const
Definition: KalFitTrack2.cxx:57
KalFitTrack::drifttime_choice_
static int drifttime_choice_
the drifttime choice
Definition: KalFitTrack.h:321
KalFitTrack::tprop_
static int tprop_
for signal propagation correction
Definition: KalFitTrack.h:277
MdcCalibFunSvc::getT0
double getT0(int layid, int cellid) const
Definition: MdcCalibFunSvc.cxx:507
MdcCalibFunSvc::getTimeWalk
double getTimeWalk(int layid, double Q) const
Definition: MdcCalibFunSvc.cxx:514
RawDataUtil::MdcTime
static double MdcTime(int timeChannel)
Definition: RawDataUtil.h:8
RawDataUtil::MdcCharge
static double MdcCharge(int chargeChannel)
Definition: RawDataUtil.h:11
RecMdcHit::getTdc
const double getTdc(void) const
Definition: RecMdcHit.h:50
RecMdcHit::getDriftT
const double getDriftT(void) const
Definition: RecMdcHit.h:52
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8

Referenced by chi2_next(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ getFiTerm()

double KalFitTrack::getFiTerm ( void  )
inline

Definition at line 164 of file KalFitTrack.h.

164{ return fiTerm_;}

Referenced by KalFitAlg::fillTds_back().

◆ getInitMatrix()

HepSymMatrix KalFitTrack::getInitMatrix ( void  ) const

Definition at line 42 of file KalFitTrack2.cxx.

43{
44 return initMatrix_ ;
45}

Referenced by KalFitAlg::smoother_calib().

◆ getPathSM()

double KalFitTrack::getPathSM ( void  )
inline

Definition at line 158 of file KalFitTrack.h.

158{ return pathSM_;}

Referenced by KalFitAlg::fillTds_back().

◆ getSigma() [1/2]

double KalFitTrack::getSigma ( int  layerId,
double  driftDist 
) const

Definition at line 3096 of file KalFitTrack.cxx.

3096 {
3097 double sigma1,sigma2,f;
3098 driftDist *= 10;//mm
3099 if(layerId<8){
3100 if(driftDist<0.5){
3101 sigma1=0.112784; sigma2=0.229274; f=0.666;
3102 }else if(driftDist<1.){
3103 sigma1=0.103123; sigma2=0.269797; f=0.934;
3104 }else if(driftDist<1.5){
3105 sigma1=0.08276; sigma2=0.17493; f=0.89;
3106 }else if(driftDist<2.){
3107 sigma1=0.070109; sigma2=0.149859; f=0.89;
3108 }else if(driftDist<2.5){
3109 sigma1=0.064453; sigma2=0.130149; f=0.886;
3110 }else if(driftDist<3.){
3111 sigma1=0.062383; sigma2=0.138806; f=0.942;
3112 }else if(driftDist<3.5){
3113 sigma1=0.061873; sigma2=0.145696; f=0.946;
3114 }else if(driftDist<4.){
3115 sigma1=0.061236; sigma2=0.119584; f=0.891;
3116 }else if(driftDist<4.5){
3117 sigma1=0.066292; sigma2=0.148426; f=0.917;
3118 }else if(driftDist<5.){
3119 sigma1=0.078074; sigma2=0.188148; f=0.911;
3120 }else if(driftDist<5.5){
3121 sigma1=0.088657; sigma2=0.27548; f=0.838;
3122 }else{
3123 sigma1=0.093089; sigma2=0.115556; f=0.367;
3124 }
3125 }else{
3126 if(driftDist<0.5){
3127 sigma1=0.112433; sigma2=0.327548; f=0.645;
3128 }else if(driftDist<1.){
3129 sigma1=0.096703; sigma2=0.305206; f=0.897;
3130 }else if(driftDist<1.5){
3131 sigma1=0.082518; sigma2=0.248913; f= 0.934;
3132 }else if(driftDist<2.){
3133 sigma1=0.072501; sigma2=0.153868; f= 0.899;
3134 }else if(driftDist<2.5){
3135 sigma1= 0.065535; sigma2=0.14246; f=0.914;
3136 }else if(driftDist<3.){
3137 sigma1=0.060497; sigma2=0.126489; f=0.918;
3138 }else if(driftDist<3.5){
3139 sigma1=0.057643; sigma2= 0.112927; f=0.892;
3140 }else if(driftDist<4.){
3141 sigma1=0.055266; sigma2=0.094833; f=0.887;
3142 }else if(driftDist<4.5){
3143 sigma1=0.056263; sigma2=0.124419; f= 0.932;
3144 }else if(driftDist<5.){
3145 sigma1=0.056599; sigma2=0.124248; f=0.923;
3146 }else if(driftDist<5.5){
3147 sigma1= 0.061377; sigma2=0.146147; f=0.964;
3148 }else if(driftDist<6.){
3149 sigma1=0.063978; sigma2=0.150591; f=0.942;
3150 }else if(driftDist<6.5){
3151 sigma1=0.072951; sigma2=0.15685; f=0.913;
3152 }else if(driftDist<7.){
3153 sigma1=0.085438; sigma2=0.255109; f=0.931;
3154 }else if(driftDist<7.5){
3155 sigma1=0.101635; sigma2=0.315529; f=0.878;
3156 }else{
3157 sigma1=0.149529; sigma2=0.374697; f=0.89;
3158 }
3159 }
3160 double sigmax = sqrt(f*sigma1*sigma1+(1 - f)*sigma2*sigma2)*0.1;
3161 return sigmax;//cm
3162}
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")

Referenced by chi2_next(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ getSigma() [2/2]

double KalFitTrack::getSigma ( KalFitHitMdc &  hitmdc,
double  tanlam,
int  lr,
double  dist 
) const

Definition at line 209 of file KalFitTrack2.cxx.

210{
211 int layerid = hitmdc.wire().layer().layerId();
212 double entrangle = hitmdc.rechitptr()->getEntra();
213 // double tanlam = hitmdc.rechitptr()->getTanl();
214 double z = hitmdc.rechitptr()->getZhit();
215 double Q = hitmdc.rechitptr()->getAdc();
216 //std::cout<<" the driftdist before getsigma is "<<dist<<" the layer is"<<layerid<<std::endl;
217 //cout<<"layerid, lr, dist, entrangle, tanlam, z , Q = "<<layerid<<", "<<lr<<", "<<dist<<", "<<entrangle<<", "<<tanlam<<", "<<z<<", "<<Q<<endl;//wangll
218 double temp = CalibFunSvc_->getSigma(layerid, lr, dist, entrangle, tanlam, z , Q );
219 //std::cout<<" the sigma is "<<temp<<std::endl;
220 return temp;
221}
MdcCalibFunSvc::getSigma
double getSigma(int layid, int lr, double dist, double entrance=0.0, double tanlam=0.0, double z=0.0, double Q=1000.0) const
Definition: MdcCalibFunSvc.cxx:278
RecMdcHit::getZhit
const double getZhit(void) const
Definition: RecMdcHit.h:55
RecMdcHit::getAdc
const double getAdc(void) const
Definition: RecMdcHit.h:51

◆ getT0()

double KalFitTrack::getT0 ( void  ) const

Definition at line 57 of file KalFitTrack2.cxx.

58{
59 //------------------get event start time-----------
60
61 if(debug_ == 4) {
62 std::cout<<"in function KalFitTrack::getT0 ( ), EventT0_ = "<<EventT0_<<std::endl;
63 }
64 return EventT0_ ;
65}

Referenced by getDriftTime().

◆ getTofSM()

double KalFitTrack::getTofSM ( void  )
inline

Definition at line 161 of file KalFitTrack.h.

161{ return tof2_;}

Referenced by KalFitAlg::fillTds_back().

◆ HelixSeg()

KalFitHelixSeg & KalFitTrack::HelixSeg ( int  i)
inline

Definition at line 235 of file KalFitTrack.h.

235{ return HelixSegs_[i];}

Referenced by KalFitAlg::smoother_calib(), and update_hits_csmalign().

◆ HelixSegs() [1/2]

void KalFitTrack::HelixSegs ( vector< KalFitHelixSeg > &  vs)

Referenced by KalFitAlg::complete_track(), KalFitAlg::fillTds_back(), KalFitAlg::filter_fwd_calib(), KalFitAlg::smoother_anal(), and KalFitAlg::smoother_calib().

◆ HelixSegs() [2/2]

vector< KalFitHelixSeg > & KalFitTrack::HelixSegs ( void  )
inline

Definition at line 234 of file KalFitTrack.h.

234{ return HelixSegs_;}

◆ HitMdc()

KalFitHitMdc & KalFitTrack::HitMdc ( int  i)
inline

Definition at line 230 of file KalFitTrack.h.

230{ return HitsMdc_[i];}

Referenced by chi2_next(), KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), KalFitAlg::smoother_anal(), smoother_Mdc_csmalign(), KalFitAlg::start_seed(), and update_hits_csmalign().

◆ HitsMdc() [1/2]

void KalFitTrack::HitsMdc ( vector< KalFitHitMdc > &  lh)

Definition at line 1836 of file KalFitTrack.cxx.

1836{ HitsMdc_ = lh;}

Referenced by KalFitAlg::complete_track(), KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew(), KalFitAlg::smoother_anal(), and KalFitAlg::start_seed().

◆ HitsMdc() [2/2]

vector< KalFitHitMdc > & KalFitTrack::HitsMdc ( void  )
inline

Definition at line 229 of file KalFitTrack.h.

229{ return HitsMdc_;}

Referenced by order_hits().

◆ insist() [1/2]

void KalFitTrack::insist ( int  t)
inline

Definition at line 209 of file KalFitTrack.h.

209{ insist_ = t;}

◆ insist() [2/2]

int KalFitTrack::insist ( void  ) const
inline

Extractor :

Definition at line 178 of file KalFitTrack.h.

178{ return insist_; }

◆ intersect_cylinder()

double KalFitTrack::intersect_cylinder ( double  r) const

Intersection with different geometry.

Definition at line 123 of file KalFitTrack.cxx.

124{
125 double m_rad = radius();
126 double l = center().perp();
127
128 double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
129
130 if(cosPhi < -1 || cosPhi > 1) return 0;
131
132 double dPhi = center().phi() - acos(cosPhi) - phi0();
133
134 if(dPhi < -M_PI) dPhi += 2 * M_PI;
135
136 return dPhi;
137}

Referenced by KalFitCylinder::intersect(), intersect_yz_plane(), and intersect_zx_plane().

◆ intersect_xy_plane()

double KalFitTrack::intersect_xy_plane ( double  z) const

Definition at line 175 of file KalFitTrack.cxx.

176{
177 if (tanl() != 0 && radius() != 0)
178 return (pivot().z() + dz() - z) / (radius() * tanl());
179 else return 0;
180}

Referenced by KalFitCylinder::intersect().

◆ intersect_yz_plane()

double KalFitTrack::intersect_yz_plane ( const HepTransform3D &  plane,
double  x 
) const

Definition at line 157 of file KalFitTrack.cxx.

159{
160 HepPoint3D xc = plane * center();
161 double r = radius();
162 double d = r * r - (x - xc.x()) * (x - xc.x());
163 if(d < 0) return 0;
164
165 double yy = xc.y();
166 if(yy > 0) yy -= sqrt(d);
167 else yy += sqrt(d);
168
169 double l = (plane.inverse() *
170 HepPoint3D(x, yy, 0)).perp();
171
172 return intersect_cylinder(l);
173}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
KalFitTrack::intersect_cylinder
double intersect_cylinder(double r) const
Intersection with different geometry.
Definition: KalFitTrack.cxx:123

◆ intersect_zx_plane()

double KalFitTrack::intersect_zx_plane ( const HepTransform3D &  plane,
double  y 
) const

Definition at line 139 of file KalFitTrack.cxx.

141{
142 HepPoint3D xc = plane * center();
143 double r = radius();
144 double d = r * r - (y - xc.y()) * (y - xc.y());
145 if(d < 0) return 0;
146
147 double xx = xc.x();
148 if(xx > 0) xx -= sqrt(d);
149 else xx += sqrt(d);
150
151 double l = (plane.inverse() *
152 HepPoint3D(xx, y, 0)).perp();
153
154 return intersect_cylinder(l);
155}

◆ lead() [1/2]

void KalFitTrack::lead ( int  i)
static

Definition at line 1827 of file KalFitTrack.cxx.

1827{ lead_ = i;}

◆ lead() [2/2]

int KalFitTrack::lead ( void  )
static

Magnetic field map.

Definition at line 1828 of file KalFitTrack.cxx.

1828{ return lead_;}

Referenced by KalFitAlg::initialize().

◆ LR()

void KalFitTrack::LR ( int  x)
static

Definition at line 1835 of file KalFitTrack.cxx.

1835{ LR_ = x;}

Referenced by KalFitAlg::initialize().

◆ mass() [1/2]

double KalFitTrack::mass ( int  i)
static

Definition at line 1821 of file KalFitTrack.cxx.

1821{ return MASS[i];}

◆ mass() [2/2]

double KalFitTrack::mass ( void  ) const
inline

Definition at line 182 of file KalFitTrack.h.

182{ return mass_; }

Referenced by KalFitAlg::fillTds(), KalFitAlg::fillTds_back(), KalFitAlg::fillTds_lead(), and KalFitAlg::smoother_anal().

◆ mom()

CLHEP::Hep3Vector * KalFitTrack::mom ( void  )
inline

Definition at line 190 of file KalFitTrack.h.

190{ return mom_; }

◆ ms()

void KalFitTrack::ms ( double  path,
const KalFitMaterial &  m,
int  index 
)

Definition at line 182 of file KalFitTrack.cxx.

184{
185 HepSymMatrix ea = Ea();
186 //cout<<"ms():path "<<path<<endl;
187 //cout<<"ms():ea before: "<<ea<<endl;
188 double k = kappa();
189 double t = tanl();
190 double t2 = t * t;
191 double pt2 = 1 + t2;
192 double k2 = k * k;
193
194 double pmag = 1 / fabs(k) * sqrt(pt2);
195 double dth = material.mcs_angle(mass_, path, pmag);
196 double dth2 = dth * dth;
197 double pt2dth2 = pt2 * dth2;
198
199 ea[1][1] += pt2dth2;
200 ea[2][2] += k2 * t2 * dth2;
201 ea[2][4] += k * t * pt2dth2;
202 ea[4][4] += pt2 * pt2dth2;
203
204 ea[3][3] += dth2 * path * path /3 / (1 + t2);
205 ea[3][4] += dth2 * path/2 * sqrt(1 + t2);
206 ea[3][2] += dth2 * t / sqrt(1 + t2) * k * path/2;
207 ea[0][0] += dth2 * path * path/3;
208 ea[0][1] += dth2 * sqrt(1 + t2) * path/2;
209
210 Ea(ea);
211 //cout<<"ms():ms angle in this: "<<dth<<endl;
212 //cout<<"ms():ea after: "<<Ea()<<endl;
213 if (index < 0) {
214 double x0 = material.X0();
215 if (x0) path_rd_ += path/x0;
216 }
217}

Referenced by KalFitElement::updateTrack(), KalFitElement::updateTrack_alreadyfound(), and KalFitElement::updateTrack_rphi().

◆ msgasmdc()

void KalFitTrack::msgasmdc ( double  path,
int  index 
)

Calculate multiple scattering angle.

Definition at line 219 of file KalFitTrack.cxx.

220{
221 double k = kappa();
222 double t = tanl();
223 double t2 = t * t;
224 double k2 = k * k;
225
226 double pmag = ( 1 / fabs(k) ) * sqrt(1 + t2);
227 double psq = pmag*pmag;
228 /*
229 double Zprims = 3/2*0.076 + 0.580/9.79*4.99*(4.99+1) +
230 0.041/183.85*74*(74+1) + 0.302/26.98 * 13 * (13+1);
231 double chicc = 0.00039612 * sqrt(Zprims * 0.001168);
232 double dth = 2.557 * chicc * sqrt(path * (mass_*mass_ + psq)) / psq;
233 */
234
235 //std::cout<<" mdcGasRadlen: "<<mdcGasRadlen_<<std::endl;
236 double pathx = path/mdcGasRadlen_;
237 double dth = 0.0136* sqrt(pathx * (mass_*mass_ + psq))/psq
238 *(1 + 0.038 * log(pathx));;
239 HepSymMatrix ea = Ea();
240#ifdef YDEBUG
241 cout<<"msgasmdc():path "<<path<<" pathx "<<pathx<<endl;
242 cout<<"msgasmdc():dth "<<dth<<endl;
243 cout<<"msgasmdc():ea before: "<<ea<<endl;
244#endif
245 double dth2 = dth * dth;
246
247 ea[1][1] += (1 + t2) * dth2;
248 ea[2][2] += k2 * t2 * dth2;
249 ea[2][4] += k * t * (1 + t2) * dth2;
250 ea[4][4] += (1 + t2) * (1 + t2) * dth2;
251
252 // additionnal terms (terms proportional to l and l^2)
253
254 ea[3][3] += dth2 * path * path /3 / (1 + t2);
255 ea[3][4] += dth2 * path/2 * sqrt(1 + t2);
256 ea[3][2] += dth2 * t / sqrt(1 + t2) * k * path/2;
257 ea[0][0] += dth2 * path * path/3;
258 ea[0][1] += dth2 * sqrt(1 + t2) * path/2;
259
260 Ea(ea);
261#ifdef YDEBUG
262 cout<<"msgasmdc():ea after: "<<Ea()<<endl;
263#endif
264 if (index < 0) {
265 pathip_ += path;
266 // RMK : put by hand, take care !!
267 double x0 = mdcGasRadlen_; // for the Mdc gas
268 path_rd_ += path/x0;
269 tof(path);
270#ifdef YDEBUG
271 cout<<"ms...pathip..."<<pathip_<<endl;
272#endif
273 }
274}
KalFitTrack::tof
double tof(void) const
Definition: KalFitTrack.h:191
KalFitTrack::mdcGasRadlen_
static double mdcGasRadlen_
Definition: KalFitTrack.h:274

Referenced by KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), KalFitAlg::smoother_anal(), and KalFitAlg::smoother_calib().

◆ ncath()

unsigned int KalFitTrack::ncath ( void  ) const
inline

Definition at line 200 of file KalFitTrack.h.

200{ return ncath_; }

◆ nchits() [1/2]

void KalFitTrack::nchits ( int  n)
inline

Definition at line 217 of file KalFitTrack.h.

217{ nchits_ = n; }

◆ nchits() [2/2]

unsigned int KalFitTrack::nchits ( void  ) const
inline

Definition at line 198 of file KalFitTrack.h.

198{ return nchits_; }

Referenced by KalFitAlg::fillTds(), KalFitAlg::fillTds_ip(), KalFitAlg::fillTds_lead(), and KalFitAlg::start_seed().

◆ ndf_back() [1/2]

void KalFitTrack::ndf_back ( int  n)
inline

Definition at line 216 of file KalFitTrack.h.

216{ ndf_back_ = n; }

◆ ndf_back() [2/2]

int KalFitTrack::ndf_back ( void  ) const
inline

Definition at line 185 of file KalFitTrack.h.

185{ return ndf_back_; }

Referenced by KalFitAlg::fillTds_back(), and KalFitAlg::start_seed().

◆ nhit_r()

int KalFitTrack::nhit_r ( void  ) const
inline

Definition at line 203 of file KalFitTrack.h.

203{ return nhit_r_; }

◆ nhit_z()

int KalFitTrack::nhit_z ( void  ) const
inline

Definition at line 204 of file KalFitTrack.h.

204{ return nhit_z_; }

◆ nLayerUsed()

int KalFitTrack::nLayerUsed ( )
inline

Definition at line 336 of file KalFitTrack.h.

337 {
338 int n=0;
339 for(int i=0; i<43; i++) n+=myLayerUsed[i];
340 return n;
341 }

Referenced by KalFitAlg::fillTds_ip().

◆ nmass()

int KalFitTrack::nmass ( void  )
static

Definition at line 1820 of file KalFitTrack.cxx.

1820{ return NMASS;}

Referenced by KalFitAlg::complete_track().

◆ nster() [1/2]

void KalFitTrack::nster ( int  n)
inline

Definition at line 218 of file KalFitTrack.h.

218{ nster_ = n; }

◆ nster() [2/2]

unsigned int KalFitTrack::nster ( void  ) const
inline

Definition at line 199 of file KalFitTrack.h.

199{ return nster_; }

Referenced by KalFitAlg::fillTds(), KalFitAlg::fillTds_ip(), KalFitAlg::fillTds_lead(), and KalFitAlg::start_seed().

◆ number_wirhit()

void KalFitTrack::number_wirhit ( void  )

Definition at line 452 of file KalFitTrack.cxx.

453{
454 unsigned int nhit = HitsMdc_.size();
455 int Num[50] = {0};
456 for( unsigned i=0 ; i < nhit; i++ )
457 Num[HitsMdc_[i].wire().layer().layerId()]++;
458 for( unsigned i=0 ; i < nhit; i++ )
459 if (Num[HitsMdc_[i].wire().layer().layerId()]>2)
460 HitsMdc_[i].chi2(-2);
461}

◆ numf() [1/2]

void KalFitTrack::numf ( int  i)
static

Definition at line 1833 of file KalFitTrack.cxx.

1833{ numf_ = i;}
KalFitTrack::numf_
static int numf_
Flag for treatment of non-uniform mag field.
Definition: KalFitTrack.h:284

◆ numf() [2/2]

int KalFitTrack::numf ( void  )
static

Definition at line 1834 of file KalFitTrack.cxx.

1834{ return numf_;}

Referenced by KalFitAlg::initialize().

◆ order_hits()

void KalFitTrack::order_hits ( void  )

Definition at line 437 of file KalFitTrack.cxx.

437 {
438 for(int it=0; it<HitsMdc().size()-1; it++){
439 if(HitsMdc_[it].wire().layer().layerId() == HitsMdc_[it+1].wire().layer().layerId())
440 {
441 if((kappa()<0)&&(HitsMdc_[it].wire().localId() > HitsMdc_[it+1].wire().localId())){
442 std::swap(HitsMdc_[it], HitsMdc_[it+1]);
443 }
444 if((kappa()>0)&&(HitsMdc_[it].wire().localId() < HitsMdc_[it+1].wire().localId())){
445 std::swap(HitsMdc_[it], HitsMdc_[it+1]);
446 }
447 }
448 }
449}
KalFitTrack::HitsMdc
vector< KalFitHitMdc > & HitsMdc(void)
Definition: KalFitTrack.h:229

◆ order_wirhit()

void KalFitTrack::order_wirhit ( int  index)

Modifier Order the wire hits for mdc track

Definition at line 375 of file KalFitTrack.cxx.

376{
377 unsigned int nhit = HitsMdc_.size();
378 Helix tracktest = *(Helix*)this;
379 int ind = 0;
380 double* Rad = new double[nhit];
381 double* Ypos = new double[nhit];
382 for( unsigned i=0 ; i < nhit; i++ ){
383
384 HepPoint3D fwd(HitsMdc_[i].wire().fwd());
385 HepPoint3D bck(HitsMdc_[i].wire().bck());
386 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
387
388 // Modification for stereo wires :
389 Helix work = tracktest;
390 work.ignoreErrorMatrix();
391 work.pivot((fwd + bck) * .5);
392 HepPoint3D x0 = (work.x(0).z() - bck.z())
393 / wire.z() * wire + bck;
394
395 tracktest.pivot(x0);
396 Rad[ind] = tracktest.x(0).perp();
397 Ypos[ind] = x0.y();
398 ind++;
399 //cout<<"Ypos: "<<Ypos[ind-1]<<endl;
400 }
401
402 // Reorder...
403 if (index < 0)
404 for(int j, k = nhit - 1; k >= 0; k = j){
405 j = -1;
406 for(int i = 1; i <= k; i++)
407 if(Rad[i - 1] > Rad[i]){
408 j = i - 1;
409 std::swap(Rad[i], Rad[j]);
410 std::swap(HitsMdc_[i], HitsMdc_[j]);
411 }
412 }
413 if (index > 0)
414 for(int j, k = nhit - 1; k >= 0; k = j){
415 j = -1;
416 for(int i = 1; i <= k; i++)
417 if(Rad[i - 1] < Rad[i]){
418 j = i - 1;
419 std::swap(Rad[i], Rad[j]);
420 std::swap(HitsMdc_[i], HitsMdc_[j]);
421 }
422 }
423 if (index == 0)
424 for(int j, k = nhit - 1; k >= 0; k = j){
425 j = -1;
426 for(int i = 1; i <= k; i++)
427 if(Ypos[i - 1] > Ypos[i]){
428 j = i - 1;
429 std::swap(Ypos[i], Ypos[j]);
430 std::swap(HitsMdc_[i], HitsMdc_[j]);
431 }
432 }
433 delete [] Rad;
434 delete [] Ypos;
435}

Referenced by KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), and KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew().

◆ p_kaon() [1/2]

void KalFitTrack::p_kaon ( double  pl)
inline

Definition at line 212 of file KalFitTrack.h.

212{ p_kaon_ = pl;}

◆ p_kaon() [2/2]

double KalFitTrack::p_kaon ( void  ) const
inline

Definition at line 194 of file KalFitTrack.h.

194{ return p_kaon_; }

Referenced by KalFitAlg::complete_track().

◆ p_proton() [1/2]

void KalFitTrack::p_proton ( double  pl)
inline

Definition at line 213 of file KalFitTrack.h.

213{ p_proton_ = pl;}

◆ p_proton() [2/2]

double KalFitTrack::p_proton ( void  ) const
inline

Definition at line 195 of file KalFitTrack.h.

195{ return p_proton_; }

Referenced by KalFitAlg::complete_track().

◆ pat1()

int KalFitTrack::pat1 ( void  ) const
inline

Definition at line 201 of file KalFitTrack.h.

201{ return pat1_; }

◆ pat2()

int KalFitTrack::pat2 ( void  ) const
inline

Definition at line 202 of file KalFitTrack.h.

202{ return pat2_; }

◆ path_ab()

double KalFitTrack::path_ab ( void  ) const
inline

Definition at line 188 of file KalFitTrack.h.

188{ return path_ab_; }

◆ path_add()

void KalFitTrack::path_add ( double  path)

Update the path length estimation.

Definition at line 1289 of file KalFitTrack.cxx.

1290{
1291 pathip_ += path;
1292 tof(path);
1293}

◆ path_rd()

double KalFitTrack::path_rd ( void  ) const
inline

Definition at line 187 of file KalFitTrack.h.

187{ return path_rd_; }

◆ pathip() [1/2]

void KalFitTrack::pathip ( double  pl)
inline

Definition at line 211 of file KalFitTrack.h.

211{ pathip_ = pl;}

◆ pathip() [2/2]

double KalFitTrack::pathip ( void  ) const
inline

Definition at line 186 of file KalFitTrack.h.

186{ return pathip_; }

Referenced by KalFitAlg::fillTds_back().

◆ PathL()

double KalFitTrack::PathL ( int  layer)

Function to calculate the path length in the layer.

◆ pathl()

double * KalFitTrack::pathl ( void  )
inline

Definition at line 189 of file KalFitTrack.h.

189{ return PathL_; }

Referenced by KalFitAlg::fillTds_back(), and pivot_numf().

◆ pivot_forMdc()

const HepPoint3D & KalFitTrack::pivot_forMdc ( void  ) const
inline

Definition at line 150 of file KalFitTrack.h.

150{ return pivot_forMdc_;}

◆ pivot_last()

const HepPoint3D & KalFitTrack::pivot_last ( void  ) const
inline

returns helix parameters

Definition at line 146 of file KalFitTrack.h.

146{ return pivot_last_; }

◆ pivot_numf() [1/2]

const HepPoint3D & KalFitTrack::pivot_numf ( const HepPoint3D &  newPivot)

Sets pivot position in a given mag field.

Definition at line 1374 of file KalFitTrack.cxx.

1374 {
1375
1376 int nstep(1);
1377 HepPoint3D delta_x((newPivot-pivot()).x()/double(inner_steps_),
1378 (newPivot-pivot()).y()/double(inner_steps_),
1379 (newPivot-pivot()).z()/double(inner_steps_));
1380 int i = 1;
1381
1382 while (i <= inner_steps_) {
1383 HepPoint3D nextPivot(pivot()+delta_x);
1384 double xnp(nextPivot.x()), ynp(nextPivot.y()), znp(nextPivot.z());
1385 HepSymMatrix Ea_now = Ea();
1386 HepPoint3D piv(pivot());
1387 double xp(piv.x()), yp(piv.y()), zp(piv.z());
1388 double dr = a()[0];
1389 double phi0 = a()[1];
1390 double kappa = a()[2];
1391 double dz = a()[3];
1392 double tanl = a()[4];
1393 double m_rad(0);
1394 if (numfcor_ == 1)
1395 m_rad = radius_numf();
1396 else
1397 m_rad = radius();
1398 double rdr = dr + m_rad;
1399 double phi = fmod(phi0 + M_PI4, M_PI2);
1400 double csf0 = cos(phi);
1401 double snf0 = (1. - csf0) * (1. + csf0);
1402 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1403 if(phi > M_PI) snf0 = - snf0;
1404
1405 double xc = xp + rdr * csf0;
1406 double yc = yp + rdr * snf0;
1407 double csf = (xc - xnp) / m_rad;
1408 double snf = (yc - ynp) / m_rad;
1409 double anrm = sqrt(csf * csf + snf * snf);
1410 csf /= anrm;
1411 snf /= anrm;
1412 phi = atan2(snf, csf);
1413 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
1414 if(phid > M_PI) phid = phid - M_PI2;
1415 double drp = (xp + dr * csf0 + m_rad * (csf0 - csf) - xnp)
1416 * csf
1417 + (yp + dr * snf0 + m_rad * (snf0 - snf) - ynp) * snf;
1418 double dzp = zp + dz - m_rad * tanl * phid - znp;
1419
1420 HepVector ap(5);
1421 ap[0] = drp;
1422 ap[1] = fmod(phi + M_PI4, M_PI2);
1423 ap[2] = kappa;
1424 ap[3] = dzp;
1425 ap[4] = tanl;
1426
1427 // Modification due to non uniform magnetic field :
1428 if (numf_ > 10) {
1429
1430 Hep3Vector x1(xp + dr*csf0, yp + dr*snf0, zp + dz);
1431 double csf0p = cos(ap[1]);
1432 double snf0p = (1. - csf0p) * (1. + csf0p);
1433 snf0p = sqrt((snf0p > 0.) ? snf0p : 0.);
1434 if(ap[1] > M_PI) snf0p = - snf0p;
1435
1436 Hep3Vector x2(xnp + drp*csf0p,
1437 ynp + drp*snf0p,
1438 znp + dzp);
1439 Hep3Vector dist((x1 - x2).x()/100.0,
1440 (x1 - x2).y()/100.0,
1441 (x1 - x2).z()/100.0);
1442 HepPoint3D middlebis((x1.x() + x2.x())/2,
1443 (x1.y() + x2.y())/2,
1444 (x1.z() + x2.z())/2);
1445 HepVector3D field;
1446 MFSvc_->fieldVector(10.*middlebis,field);
1447 field = 1000.*field;
1448 Hep3Vector dB(field.x(),
1449 field.y(),
1450 (field.z()-0.1*Bznom_));
1451 if (field.z()) {
1452 double akappa(fabs(kappa));
1453 double sign = kappa/akappa;
1454 HepVector dp(3);
1455 dp = 0.299792458 * sign * dB.cross(dist);
1456 HepVector dhp(3);
1457 dhp[0] = -akappa*(dp[0]*csf0p+dp[1]*snf0p);
1458 dhp[1] = kappa*akappa*(dp[0]*snf0p-dp[1]*csf0p);
1459 dhp[2] = dp[2]*akappa+dhp[1]*tanl/kappa;
1460 if (numfcor_ == 0){
1461 ap[1] += dhp[0];
1462 }
1463 ap[2] += dhp[1];
1464 ap[4] += dhp[2];
1465 }
1466 }
1467 HepMatrix m_del = delApDelA(ap);
1468 Ea_now.assign(m_del * Ea_now * m_del.T());
1469 pivot(nextPivot);
1470 a(ap);
1471 Ea(Ea_now);
1472 i++;
1473 }
1474 return newPivot;
1475}
IMagneticFieldSvc::fieldVector
virtual StatusCode fieldVector(const HepGeom::Point3D< double > &xyz, HepGeom::Vector3D< double > &fvec) const =0
KalFitTrack::inner_steps_
static int inner_steps_
Definition: KalFitTrack.h:285
KalFitTrack::numfcor_
static int numfcor_
NUMF treatment improved.
Definition: KalFitTrack.h:297
KalFitTrack::radius_numf
double radius_numf(void) const
Estimation of the radius in a given mag field.
Definition: KalFitTrack.cxx:1336
KalmanFit::Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/helix/Helix.h:218

Referenced by KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), KalFitAlg::smoother_anal(), KalFitAlg::smoother_calib(), KalFitElement::updateTrack(), KalFitElement::updateTrack_alreadyfound(), and KalFitElement::updateTrack_rphi().

◆ pivot_numf() [2/2]

const HepPoint3D & KalFitTrack::pivot_numf ( const HepPoint3D &  newPivot,
double &  pathl 
)

Definition at line 1478 of file KalFitTrack.cxx.

1478 {
1479
1480 Helix tracktest = *(Helix*)this;
1481 tracktest.ignoreErrorMatrix();
1482 double tl = a()[4];
1483 double th = 90.0 - 180.0*M_1_PI*atan( tl );
1484 /*
1485 int nstep(1);
1486 if (steplev_ == 1)
1487 nstep = 3;
1488 else if (steplev_ == 2 && (th > 140 || th <45))
1489 if ((pivot()-newPivot).mag()<3.)
1490 nstep = 3;
1491 else
1492 nstep = 6;
1493 */
1494 Hep3Vector delta_x((newPivot-pivot()).x()/double(outer_steps_),
1495 (newPivot-pivot()).y()/double(outer_steps_),
1496 (newPivot-pivot()).z()/double(outer_steps_));
1497 int i = 1;
1498
1499 while (i <= outer_steps_) {
1500 HepPoint3D nextPivot(pivot()+delta_x);
1501 double xnp(nextPivot.x()), ynp(nextPivot.y()), znp(nextPivot.z());
1502
1503 HepSymMatrix Ea_now = Ea();
1504 HepPoint3D piv(pivot());
1505 double xp(piv.x()), yp(piv.y()), zp(piv.z());
1506
1507 double dr = a()[0];
1508 double phi0 = a()[1];
1509 double kappa = a()[2];
1510 double dz = a()[3];
1511 double tanl = a()[4];
1512
1513 double m_rad(0);
1514 m_rad = radius();
1515
1516 double rdr = dr + m_rad;
1517 double phi = fmod(phi0 + M_PI4, M_PI2);
1518 double csf0 = cos(phi);
1519 double snf0 = (1. - csf0) * (1. + csf0);
1520 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1521 if(phi > M_PI) snf0 = - snf0;
1522
1523 double xc = xp + rdr * csf0;
1524 double yc = yp + rdr * snf0;
1525 double csf = (xc - xnp) / m_rad;
1526 double snf = (yc - ynp) / m_rad;
1527 double anrm = sqrt(csf * csf + snf * snf);
1528 csf /= anrm;
1529 snf /= anrm;
1530 phi = atan2(snf, csf);
1531 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
1532 if(phid > M_PI) phid = phid - M_PI2;
1533 double drp = (xp + dr * csf0 + m_rad * (csf0 - csf) - xnp)
1534 * csf
1535 + (yp + dr * snf0 + m_rad * (snf0 - snf) - ynp) * snf;
1536 double dzp = zp + dz - m_rad * tanl * phid - znp;
1537
1538 HepVector ap(5);
1539 ap[0] = drp;
1540 ap[1] = fmod(phi + M_PI4, M_PI2);
1541 ap[2] = kappa;
1542 ap[3] = dzp;
1543 ap[4] = tanl;
1544
1545 //std::cout<<" numf_: "<<numf_<<std::endl;
1546
1547 // Modification due to non uniform magnetic field :
1548 if (numf_ > 10) {
1549
1550 Hep3Vector x1(xp + dr*csf0, yp + dr*snf0, zp + dz);
1551
1552 double csf0p = cos(ap[1]);
1553 double snf0p = (1. - csf0p) * (1. + csf0p);
1554 snf0p = sqrt((snf0p > 0.) ? snf0p : 0.);
1555 if(ap[1] > M_PI) snf0p = - snf0p;
1556
1557 Hep3Vector x2(xnp + drp*csf0p,
1558 ynp + drp*snf0p,
1559 znp + dzp);
1560
1561 Hep3Vector dist((x1 - x2).x()/100.0,
1562 (x1 - x2).y()/100.0,
1563 (x1 - x2).z()/100.0);
1564
1565 HepPoint3D middlebis((x1.x() + x2.x())/2,
1566 (x1.y() + x2.y())/2,
1567 (x1.z() + x2.z())/2);
1568
1569 HepVector3D field;
1570 MFSvc_->fieldVector(10.*middlebis,field);
1571 field = 1000.*field;
1572
1573 //std::cout<<"B field: "<<field<<std::endl;
1574 Hep3Vector dB(field.x(),
1575 field.y(),
1576 (field.z()-0.1*Bznom_));
1577
1578
1579 //std::cout<<" dB: "<<dB<<std::endl;
1580
1581
1582 if (field.z()) {
1583 double akappa(fabs(kappa));
1584 double sign = kappa/akappa;
1585 HepVector dp(3);
1586 dp = 0.299792458 * sign * dB.cross(dist);
1587
1588 //std::cout<<"dp: "<<dp<<std::endl;
1589
1590 HepVector dhp(3);
1591 dhp[0] = -akappa*(dp[0]*csf0p+dp[1]*snf0p);
1592 dhp[1] = kappa*akappa*(dp[0]*snf0p-dp[1]*csf0p);
1593 dhp[2] = dp[2]*akappa+dhp[1]*tanl/kappa;
1594
1595
1596 //std::cout<<"dhp: "<<dhp<<std::endl;
1597
1598
1599 ap[1] += dhp[0];
1600 ap[2] += dhp[1];
1601 ap[4] += dhp[2];
1602 }
1603 }
1604 HepMatrix m_del = delApDelA(ap);
1605 Ea_now.assign(m_del * Ea_now * m_del.T());
1606 pivot(nextPivot);
1607 a(ap);
1608 Ea(Ea_now);
1609 i++;
1610
1611 //std::cout<<" a: "<<a()<<std::endl;
1612 }
1613
1614 // Estimation of the path length :
1615 double tanl_0(tracktest.a()[4]);
1616 double phi0_0(tracktest.a()[1]);
1617 double radius_0(tracktest.radius());
1618 tracktest.pivot(newPivot);
1619
1620 double phi0_1 = tracktest.a()[1];
1621 if (fabs(phi0_1 - phi0_0) > M_PI) {
1622 if (phi0_1 > phi0_0) phi0_1 -= 2 * M_PI;
1623 else phi0_0 -= 2 * M_PI;
1624 }
1625 if(phi0_1 == phi0_0) phi0_1 = phi0_0 + 1.e-10;
1626 pathl = fabs(radius_0 * (phi0_1 - phi0_0)
1627 * sqrt(1 + tanl_0 * tanl_0));
1628 return newPivot;
1629}
KalFitTrack::pathl
double * pathl(void)
Definition: KalFitTrack.h:189
KalFitTrack::outer_steps_
static int outer_steps_
Definition: KalFitTrack.h:286

◆ point_last() [1/2]

void KalFitTrack::point_last ( const HepPoint3D &  point)
inline

set and give out the last point of the track

Definition at line 141 of file KalFitTrack.h.

141{ point_last_ = point;}

Referenced by KalFitAlg::fillTds_back(), and KalFitAlg::smoother_anal().

◆ point_last() [2/2]

const HepPoint3D & KalFitTrack::point_last ( void  )
inline

Definition at line 142 of file KalFitTrack.h.

142{ return point_last_;}

◆ r0()

double KalFitTrack::r0 ( void  ) const
inline

Definition at line 181 of file KalFitTrack.h.

181{ return r0_; }

◆ r_max()

double KalFitTrack::r_max ( void  ) const
inline

Definition at line 197 of file KalFitTrack.h.

197{ return r_max_; }

◆ radius_numf()

double KalFitTrack::radius_numf ( void  ) const

Estimation of the radius in a given mag field.

Definition at line 1336 of file KalFitTrack.cxx.

1336 {
1337
1338 double Bz(Bznom_);
1339 //std::cout<<"Bz: "<<Bz<<std::endl;
1340 if (numf_ > 10){
1341 double dr = a()[0];
1342 double phi0 = a()[1];
1343 double dz = a()[3];
1344 double phi = fmod(phi0 + M_PI4, M_PI2);
1345 double csf0 = cos(phi);
1346 double snf0 = (1. - csf0) * (1. + csf0);
1347 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1348 if(phi > M_PI) snf0 = - snf0;
1349 //XYZPoint
1350 HepPoint3D x0((pivot().x() + dr*csf0),
1351 (pivot().y() + dr*snf0),
1352 (pivot().z() + dz));
1353
1354 //XYZVector b;
1355 HepVector3D b;
1356
1357 //std::cout<<"b: "<<b<<std::endl;
1358
1359 MFSvc_->fieldVector(10.*x0, b);
1360
1361 //std::cout<<"b: "<<b<<std::endl;
1362
1363
1364 b = 10000.*b;
1365 Bz = b.z();
1366 }
1367 if (Bz == 0)
1368 Bz = Bznom_;
1369 double ALPHA_loc = 10000./2.99792458/Bz;
1370 return ALPHA_loc / a()[2];
1371}
b
const double b
Definition: slope.cxx:9

Referenced by KalFitAlg::fillTds_back(), and pivot_numf().

◆ resetLayerUsed()

void KalFitTrack::resetLayerUsed ( )
inline

Definition at line 343 of file KalFitTrack.h.

343 {
344 for(int i=0; i<43; i++) myLayerUsed[i]=0;
345 }

◆ resol() [1/2]

void KalFitTrack::resol ( int  i)
static

Definition at line 1831 of file KalFitTrack.cxx.

1831{ resolflag_ = i;}

◆ resol() [2/2]

int KalFitTrack::resol ( void  )
static

Definition at line 1832 of file KalFitTrack.cxx.

1832{ return resolflag_;}

Referenced by KalFitAlg::initialize().

◆ setIMdcGeomSvc()

void KalFitTrack::setIMdcGeomSvc ( IMdcGeomSvc *  igeomsvc)
static

Definition at line 239 of file KalFitTrack2.cxx.

240{
241 iGeomSvc_ = igeomsvc;
242}

Referenced by KalFitAlg::setGeomSvc_init().

◆ setInitMatrix()

void KalFitTrack::setInitMatrix ( HepSymMatrix  m)
static

Definition at line 37 of file KalFitTrack2.cxx.

38{
39 initMatrix_ = m;
40}

Referenced by KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), and KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew().

◆ setMagneticFieldSvc()

void KalFitTrack::setMagneticFieldSvc ( IMagneticFieldSvc *  mf)
static

Definition at line 228 of file KalFitTrack2.cxx.

229{
230 /*ISvcLocator* svcLocator = Gaudi::svcLocator();
231 StatusCode sc = svcLocator->service("MagneticFieldSvc",MFSvc_);
232 if (sc.isFailure()){
233 std::cout << "Could not load MagneticFieldSvc!" << std::endl;
234 }*/
235 MFSvc_ = mf;
236 if(MFSvc_==0) cout<<"KalFitTrack2:: Could not load MagneticFieldSvc!"<<endl;
237}

Referenced by KalFitAlg::initialize().

◆ setMdcCalibFunSvc()

void KalFitTrack::setMdcCalibFunSvc ( const MdcCalibFunSvc *  calibsvc)
static

Definition at line 223 of file KalFitTrack2.cxx.

224{
225 CalibFunSvc_ = calibsvc;
226}

Referenced by KalFitAlg::setCalibSvc_init().

◆ setMdcDigiCol()

void KalFitTrack::setMdcDigiCol ( MdcDigiCol *  digicol)
static

Definition at line 244 of file KalFitTrack2.cxx.

245{
246 mdcDigiCol_ = digicol;
247}

Referenced by KalFitAlg::execute().

◆ setT0()

void KalFitTrack::setT0 ( double  t0)
static

Definition at line 47 of file KalFitTrack2.cxx.

48{
49 //------------------set event start time-----------
50
51 EventT0_ = eventstart;
52 if(debug_ == 4) {
53 std::cout<<"in function KalFitTrack::setT0(...), EventT0_ = "<<EventT0_<<std::endl;
54 }
55}

Referenced by KalFitAlg::execute().

◆ smoother_Mdc() [1/2]

double KalFitTrack::smoother_Mdc ( KalFitHelixSeg &  seg,
CLHEP::Hep3Vector &  meas,
int &  flg,
int  csmflag 
)

Kalman smoother for Mdc.

Referenced by KalFitAlg::smoother_anal(), and KalFitAlg::smoother_calib().

◆ smoother_Mdc() [2/2]

double KalFitTrack::smoother_Mdc ( KalFitHitMdc &  HitMdc,
CLHEP::Hep3Vector &  meas,
KalFitHelixSeg &  seg,
double &  dchi2,
int  csmflag 
)

◆ smoother_Mdc_csmalign()

double KalFitTrack::smoother_Mdc_csmalign ( KalFitHelixSeg &  seg,
CLHEP::Hep3Vector &  meas,
int &  flg,
int  csmflag 
)

move the pivot of the helixseg to IP(0,0,0)

Definition at line 656 of file KalFitTrack.cxx.

657{
658
659
660 HepPoint3D ip(0., 0., 0.);
661 // attention! we should not to decide the left&right in the smoother process,
662 // because we choose the left&right of hits from the filter process.
663
664 KalFitHitMdc* HitMdc = seg.HitMdc();
665 double lr = HitMdc->LR();
666
667 // For taking the propagation delay into account :
668 int layerid = (*HitMdc).wire().layer().layerId();
669 int wire_ID = HitMdc->wire().geoID();
670 double entrangle = HitMdc->rechitptr()->getEntra();
671
672 if (wire_ID<0 || wire_ID>6796){ //bes
673 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID << std::endl;
674 return 0;
675 }
676
677 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
678 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
679 double dd(0.);
680 double tofest(0);
681 double phi = fmod(phi0() + M_PI4, M_PI2);
682 double csf0 = cos(phi);
683 double snf0 = (1. - csf0) * (1. + csf0);
684 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
685 if(phi > M_PI) snf0 = - snf0;
686
687 if (Tof_correc_) {
688 Hep3Vector ip(0, 0, 0);
689 Helix work = *(Helix*)this;
690 work.ignoreErrorMatrix();
691 work.pivot(ip);
692 double phi_ip = work.phi0();
693 if (fabs(phi - phi_ip) > M_PI) {
694 if (phi > phi_ip) phi -= 2 * M_PI;
695 else phi_ip -= 2 * M_PI;
696 }
697 double t = tanl();
698 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
699 double pmag( sqrt( 1.0 + t*t ) / kappa());
700 double mass_over_p( mass_ / pmag );
701 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
702 tofest = l / ( 29.9792458 * beta );
703 if(csmflag==1 && (*HitMdc).wire().y()>0.) tofest= -1. * tofest;
704 }
705
706 const HepSymMatrix& ea = Ea();
707 const HepVector& v_a = a();
708
709
710 HepPoint3D pivot_work = pivot();
711
712 /*
713 HepVector a_work = a();
714 HepSymMatrix ea_work = Ea();
715
716 KalFitTrack helix_work(pivot_work, a_work, ea_work, 0, 0, 0);
717 helix_work.pivot(ip);
718
719 HepVector a_temp = helix_work.a();
720 HepSymMatrix ea_temp = helix_work.Ea();
721
722 seg.Ea_pre_bwd(ea_temp);
723 seg.a_pre_bwd(a_temp);
724
725 */
726
727 seg.a_pre_bwd(a());
728 seg.Ea_pre_bwd(Ea());
729
730 double dchi2R(99999999), dchi2L(99999999);
731 HepVector v_H(5, 0);
732 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
733 v_H[3] = -meas.z();
734 HepMatrix v_HT = v_H.T();
735
736 double estim = (v_HT * v_a)[0];
737 HepVector ea_v_H = ea * v_H;
738 HepMatrix ea_v_HT = (ea_v_H).T();
739 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
740 HepSymMatrix eaNew(5);
741 HepVector aNew(5);
742 double time = (*HitMdc).tdc();
743
744 //seg.dt(time);
745 // if (Tof_correc_)
746 // {
747 // time -= tofest;
748 // seg.dt(time);
749 // }
750 // double dd = 1.0e-4 * 40.0 * time;
751 // seg.dd(dd);
752
753 double drifttime = getDriftTime(*HitMdc , tofest);
754 seg.dt(drifttime);
755 double ddl = getDriftDist(*HitMdc, drifttime, 0 );
756 double ddr = getDriftDist(*HitMdc, drifttime, 1 );
757 double erddl = getSigma( *HitMdc, a()[4], 0, ddl);
758 double erddr = getSigma( *HitMdc, a()[4], 1, ddr);
759
760 double dmeas2[2] = { 0.1*ddl, 0.1*ddr }; // millimeter to centimeter
761 double er_dmeas2[2] = {0., 0.};
762 if(resolflag_ == 1) {
763 er_dmeas2[0] = 0.1*erddl;
764 er_dmeas2[1] = 0.1*erddr;
765 }else if(resolflag_ == 0)
766 {
767 }
768
769 // Left hypothesis :
770 if (lr == -1) {
771 double er_dmeasL, dmeasL;
772 if(Tof_correc_) {
773 dmeasL = (-1.0)*dmeas2[0];
774 er_dmeasL = er_dmeas2[0];
775 } else {
776 dmeasL = (-1.0)*fabs((*HitMdc).dist()[0]);
777 er_dmeasL = (*HitMdc).erdist()[0];
778 }
779
780
781 //if(layerid < 4) er_dmeasL*=2.0;
782
783 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
784 eaNew.assign(ea - ea_v_H * AL * ea_v_HT);
785 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
786 if(0. == RkL) RkL = 1.e-4;
787
788 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
789 aNew = v_a + diffL;
790 double sigL = (dmeasL - (v_HT * aNew)[0]);
791 dchi2L = (sigL*sigL) / (RkL * er_dmeasL*er_dmeasL);
792 } else if (lr == 1) {
793 // Right hypothesis :
794
795 double er_dmeasR, dmeasR;
796 if(Tof_correc_) {
797 dmeasR = dmeas2[1];
798 er_dmeasR = er_dmeas2[1];
799 } else {
800 dmeasR = fabs((*HitMdc).dist()[1]);
801 er_dmeasR = (*HitMdc).erdist()[1];
802 }
803
804
805 //if(layerid < 4) er_dmeasR*=2.0;
806
807
808 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
809 eaNew.assign(ea - ea_v_H * AR * ea_v_HT);
810 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
811 if(0. == RkR) RkR = 1.e-4;
812
813 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
814 aNew = v_a + diffR;
815 double sigR = (dmeasR- (v_HT * aNew)[0]);
816 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
817 }
818
819 // Update Kalman result :
820#ifdef YDEBUG
821 cout<<"in smoother_Mdc: lr= "<<lr<<" a: "<<v_a<<" a_NEW: "<<aNew<<endl;
822#endif
823 double dchi2_loc(0);
824 if ((dchi2R < dchi2cuts_calib[layerid] && lr == 1) ||
825 (dchi2L < dchi2cuts_calib[layerid] && lr == -1)) {
826
827 ndf_back_++;
828 flg = 1;
829 int chge(1);
830 if (!(aNew[2] && fabs(aNew[2]-a()[2]) < 1.0)) chge = 0;
831 if (lr == 1) {
832 chiSq_back_ += dchi2R;
833 dchi2_loc = dchi2R;
834 dd = 0.1*ddr;
835 // if(debug_ ==4) std::cout<<"in smoother "<<std::endl;
836
837 } else if (lr == -1) {
838 chiSq_back_ += dchi2L;
839 dchi2_loc = dchi2L;
840 dd = -0.1*ddl;
841
842 }
843 if (chge){
844 a(aNew);
845 Ea(eaNew);
846 }
847
848 seg.a_filt_bwd(aNew);
849 seg.Ea_filt_bwd(eaNew);
850
851 HepVector a_pre_fwd_temp=seg.a_pre_fwd();
852 if( (a_pre_fwd_temp[1]-seg.a_pre_bwd()[1]) > 3. * M_PI /2.) a_pre_fwd_temp[1]-= M_PI2;
853 if( (a_pre_fwd_temp[1]-seg.a_pre_bwd()[1]) < -3. * M_PI /2. ) a_pre_fwd_temp[1]+= M_PI2;
854
855 seg.a_pre_fwd(a_pre_fwd_temp);
856
857 HepVector a_pre_filt_temp=seg.a_filt_fwd();
858 if( (a_pre_filt_temp[1]-seg.a_pre_bwd()[1]) > 3. * M_PI /2. ) a_pre_filt_temp[1] -= M_PI2;
859 if( (a_pre_filt_temp[1]-seg.a_pre_bwd()[1]) < -3. * M_PI /2.) a_pre_filt_temp[1] += M_PI2;
860 seg.a_filt_fwd(a_pre_filt_temp);
861
862 /*
863 KalFitTrack helixNew(pivot_work, aNew, eaNew, 0, 0, 0);
864 helixNew.pivot(ip);
865 a_temp = helixNew.a();
866 ea_temp = helixNew.Ea();
867 seg.Ea_filt_bwd(ea_temp);
868 seg.a_filt_bwd(a_temp);
869 */
870
871 int inv;
872
873 if(debug_ == 4){
874 std::cout<<"seg.Ea_pre_bwd().inverse(inv) ..."<<seg.Ea_pre_bwd().inverse(inv)<<std::endl;
875 std::cout<<"seg.Ea_pre_fwd().inverse(inv) ..."<<seg.Ea_pre_fwd().inverse(inv)<<std::endl;
876 }
877
878 HepSymMatrix Ea_pre_comb = (seg.Ea_pre_bwd().inverse(inv)+seg.Ea_pre_fwd().inverse(inv)).inverse(inv);
879 seg.Ea_exclude(Ea_pre_comb);
880
881
882 if(debug_ == 4) {
883 std::cout<<"Ea_pre_comb_ ... "<<Ea_pre_comb<<std::endl;
884 std::cout<<"seg.a_pre_bwd()..."<<seg.a_pre_bwd()<<std::endl;
885 std::cout<<"seg.a_pre_fwd()..."<<seg.a_pre_fwd()<<std::endl;
886 }
887 HepVector a_pre_comb = Ea_pre_comb*((seg.Ea_pre_bwd().inverse(inv))*seg.a_pre_bwd()+(seg.Ea_pre_fwd().inverse(inv))*seg.a_pre_fwd());
888 seg.a_exclude(a_pre_comb);
889
890 if(debug_ == 4) {
891 std::cout<<"seg.Ea_filt_fwd()..."<<seg.Ea_filt_fwd()<<std::endl;
892 std::cout<<"seg.Ea_filt_fwd().inverse(inv)..."<<seg.Ea_filt_fwd().inverse(inv)<<std::endl;
893 }
894 seg.Ea((seg.Ea_filt_fwd().inverse(inv)+seg.Ea_pre_bwd().inverse(inv)).inverse(inv));
895 seg.Ea_include((seg.Ea_filt_fwd().inverse(inv)+seg.Ea_pre_bwd().inverse(inv)).inverse(inv));
896
897 if(debug_ == 4) {
898 std::cout<<"seg.Ea() is ..."<<seg.Ea()<<std::endl;
899 std::cout<<"seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd() ..."<<seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()<<std::endl;
900 std::cout<<"seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd() ... "<<seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()<<std::endl;
901 }
902
903 seg.a(seg.Ea()*(seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()+seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()));
904 seg.a_include(seg.Ea()*(seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()+seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()));
905
906 if(inv != 0) {
907 //std::cout<<"ERROR OCCUR WHEN INVERSE MATRIX !"<<std::endl;
908 }
909
910 seg.residual_exclude(dd - (v_HT*a_pre_comb)[0]);
911 seg.residual_include(dd - (v_HT*seg.a())[0]);
912 seg.doca_exclude(fabs((v_HT*a_pre_comb)[0]));
913 seg.doca_include(fabs((v_HT*seg.a())[0]));
914
915 if(debug_ == 4){
916 std::cout<<"the dd in smoother_Mdc is "<<dd
917 <<" the (v_HT*a_pre_comb)[0] is "<<(v_HT*a_pre_comb)[0]<<std::endl;
918 }
919
920 //compare the two method to calculate the include doca value :
921 if(debug_ == 4){
922 std::cout<<"method 1 ..."<<sqrt(seg.a()[0]*seg.a()[0]+seg.a()[3]*seg.a()[3])<<std::endl;
923 std::cout<<"method 2 ..."<<fabs((v_HT*seg.a())[0])<<std::endl;
924 }
925
926
927 /// move the pivot of the helixseg to IP(0,0,0)
928 KalFitTrack helixNew1(pivot_work, seg.a(), seg.Ea(), 0, 0, 0);
929 helixNew1.pivot(ip);
930 HepVector a_temp1 = helixNew1.a();
931 HepSymMatrix ea_temp1 = helixNew1.Ea();
932 seg.Ea(ea_temp1);
933 seg.a(a_temp1);
934 seg.Ea_include(ea_temp1);
935 seg.a_include(a_temp1);
936
937 KalFitTrack helixNew2(pivot_work, seg.a_exclude(), seg.Ea_exclude(), 0, 0, 0);
938 helixNew2.pivot(ip);
939 HepVector a_temp2 = helixNew2.a();
940 HepSymMatrix ea_temp2 = helixNew2.Ea();
941 seg.Ea_exclude(ea_temp2);
942 seg.a_exclude(a_temp2);
943
944 seg.tof(tofest);
945 seg.dd(dd);
946
947 }
948 return chiSq_back_;
949}
KalFitHelixSeg::a_pre_fwd
CLHEP::HepVector a_pre_fwd(void)
Definition: KalFitHelixSeg.h:45
KalFitHelixSeg::Ea_filt_fwd
CLHEP::HepSymMatrix Ea_filt_fwd(void)
Definition: KalFitHelixSeg.h:54
KalFitHelixSeg::a_filt_bwd
CLHEP::HepVector a_filt_bwd(void)
Definition: KalFitHelixSeg.h:48
KalFitHelixSeg::a_include
CLHEP::HepVector a_include(void)
Definition: KalFitHelixSeg.h:49
KalFitHelixSeg::tof
double tof(void)
Definition: KalFitHelixSeg.h:31
KalFitHelixSeg::Ea_pre_bwd
CLHEP::HepSymMatrix Ea_pre_bwd(void)
Definition: KalFitHelixSeg.h:55
KalFitHelixSeg::doca_include
double doca_include(void)
Definition: KalFitHelixSeg.h:30
KalFitHelixSeg::a_filt_fwd
CLHEP::HepVector a_filt_fwd(void)
Definition: KalFitHelixSeg.h:47
KalFitHelixSeg::dt
double dt(void)
Definition: KalFitHelixSeg.h:32
KalFitHelixSeg::doca_exclude
double doca_exclude(void)
Definition: KalFitHelixSeg.h:29
KalFitHelixSeg::a_pre_bwd
CLHEP::HepVector a_pre_bwd(void)
Definition: KalFitHelixSeg.h:46
KalFitHelixSeg::Ea_pre_fwd
CLHEP::HepSymMatrix & Ea_pre_fwd(void)
Definition: KalFitHelixSeg.h:53
KalFitHelixSeg::dd
double dd(void)
Definition: KalFitHelixSeg.h:33
KalFitHelixSeg::a_exclude
CLHEP::HepVector a_exclude(void)
Definition: KalFitHelixSeg.h:50
KalFitHelixSeg::Ea_filt_bwd
CLHEP::HepSymMatrix Ea_filt_bwd(void)
Definition: KalFitHelixSeg.h:56
KalFitHelixSeg::residual_exclude
double residual_exclude(void)
Definition: KalFitHelixSeg.h:27
KalFitHelixSeg::Ea_include
CLHEP::HepSymMatrix Ea_include(void)
Definition: KalFitHelixSeg.h:57
KalFitHelixSeg::Ea_exclude
CLHEP::HepSymMatrix Ea_exclude(void)
Definition: KalFitHelixSeg.h:58
KalFitHelixSeg::HitMdc
KalFitHitMdc * HitMdc(void)
Definition: KalFitHelixSeg.h:26
KalFitHelixSeg::residual_include
double residual_include(void)
Definition: KalFitHelixSeg.h:28
KalFitHitMdc
Description of a Hit in Mdc.
Definition: KalFitHitMdc.h:17
KalFitTrack
Description of a track class (<- Helix.cc)
Definition: KalFitTrack.h:35
KalFitTrack::dchi2cuts_calib
static double dchi2cuts_calib[43]
Definition: KalFitTrack.h:294

Referenced by KalFitAlg::smoother_calib().

◆ tof() [1/2]

void KalFitTrack::tof ( double  path)

Update the tof estimation.

Definition at line 1311 of file KalFitTrack.cxx.

1312{
1313 double light_speed( 29.9792458 ); // light speed in cm/nsec
1314 double t = tanl();
1315 double pmag( sqrt( 1.0 + t*t ) / kappa());
1316 if (pmag !=0) {
1317 double mass_over_p( mass_ / pmag );
1318 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
1319 tof_ += path / ( light_speed * beta );
1320 }
1321
1322 if (tofall_) {
1323 if (p_kaon_ > 0){
1324 double massk_over_p( MASS[3] / p_kaon_ );
1325 double beta_kaon( 1.0 / sqrt( 1.0 + massk_over_p * massk_over_p ) );
1326 tof_kaon_ += path / (light_speed * beta_kaon);
1327 }
1328 if (p_proton_ > 0){
1329 double massp_over_p( MASS[4] / p_proton_ );
1330 double beta_proton( 1.0 / sqrt( 1.0 + massp_over_p * massp_over_p ) );
1331 tof_proton_ += path / (light_speed * beta_proton);
1332 }
1333 }
1334}

Referenced by KalFitAlg::fillTds_back().

◆ tof() [2/2]

double KalFitTrack::tof ( void  ) const
inline

Definition at line 191 of file KalFitTrack.h.

191{ return tof_; }

Referenced by msgasmdc(), and path_add().

◆ tof_kaon()

double KalFitTrack::tof_kaon ( void  ) const
inline

Definition at line 192 of file KalFitTrack.h.

192{ return tof_kaon_; }

Referenced by KalFitAlg::fillTds_back().

◆ tof_proton()

double KalFitTrack::tof_proton ( void  ) const
inline

Definition at line 193 of file KalFitTrack.h.

193{ return tof_proton_; }

Referenced by KalFitAlg::fillTds_back().

◆ trasan_id() [1/2]

void KalFitTrack::trasan_id ( int  t)
inline

Definition at line 208 of file KalFitTrack.h.

208{ trasan_id_ = t;}

◆ trasan_id() [2/2]

int KalFitTrack::trasan_id ( void  ) const
inline

Definition at line 180 of file KalFitTrack.h.

180{ return trasan_id_; }

◆ type() [1/2]

void KalFitTrack::type ( int  t)
inline

Reinitialize (modificator)

Definition at line 207 of file KalFitTrack.h.

207{ type_ = t;}

◆ type() [2/2]

int KalFitTrack::type ( void  ) const
inline

Definition at line 179 of file KalFitTrack.h.

179{ return type_; }

Referenced by KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), and KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew().

◆ update_bit()

void KalFitTrack::update_bit ( int  i)

Definition at line 1807 of file KalFitTrack.cxx.

1807 {
1808 int j(0);
1809 if (i < 31){
1810 j = (int) pow(2.,i);
1811 if (!(pat1_ & j))
1812 pat1_ = pat1_ | j;
1813 } else if (i < 50) {
1814 j = (int) pow(2.,(i-31));
1815 if (!(pat2_ & j))
1816 pat2_ = pat2_ | j;
1817 }
1818}

Referenced by update_hits_csmalign().

◆ update_forMdc()

void KalFitTrack::update_forMdc ( void  )

Definition at line 116 of file KalFitTrack.cxx.

117{
118 pivot_forMdc_ = pivot();
119 a_forMdc_ = a();
120 Ea_forMdc_ = Ea();
121}

Referenced by KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), and KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew().

◆ update_hits() [1/2]

double KalFitTrack::update_hits ( KalFitHelixSeg &  HelixSeg,
int  inext,
CLHEP::Hep3Vector &  meas,
int  way,
double &  dchi2,
int  csmflag 
)

◆ update_hits() [2/2]

double KalFitTrack::update_hits ( KalFitHitMdc &  HitMdc,
int  inext,
CLHEP::Hep3Vector &  meas,
int  way,
double &  dchi2,
double &  dtrack,
double &  dtracknew,
double &  dtdc,
int  csmflag 
)

Include the Mdc wire hits.

Referenced by KalFitAlg::filter_fwd_anal(), and KalFitAlg::filter_fwd_calib().

◆ update_hits_csmalign()

double KalFitTrack::update_hits_csmalign ( KalFitHelixSeg &  HelixSeg,
int  inext,
CLHEP::Hep3Vector &  meas,
int  way,
double &  dchi2,
int  csmflag 
)

Definition at line 2104 of file KalFitTrack.cxx.

2104 {
2105
2106
2107 HepPoint3D ip(0.,0.,0.);
2108
2109 KalFitHitMdc* HitMdc = HelixSeg.HitMdc();
2110 double lr = HitMdc->LR();
2111 int layerid = (*HitMdc).wire().layer().layerId();
2112 // cout<<"layerid: "<<layerid<<endl;
2113 double entrangle = HitMdc->rechitptr()->getEntra();
2114
2115 if(debug_ == 4) {
2116 std::cout<<"in update_hits(kalFitHelixSeg,...) : the layerid ="<<layerid<<std::endl;
2117 std::cout<<" update_hits: lr= " << lr <<std::endl;
2118 }
2119 // For taking the propagation delay into account :
2120 const KalFitWire& Wire = HitMdc->wire();
2121 int wire_ID = Wire.geoID();
2122
2123 if (wire_ID<0 || wire_ID>6796){ //bes
2124 std::cout << " KalFitTrack: wire_ID problem : " << wire_ID
2125 << std::endl;
2126 return 0;
2127 }
2128 int flag_ster = Wire.stereo();
2129 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
2130 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
2131 double tofest(0);
2132 double phi = fmod(phi0() + M_PI4, M_PI2);
2133 double csf0 = cos(phi);
2134 double snf0 = (1. - csf0) * (1. + csf0);
2135 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
2136 if(phi > M_PI) snf0 = - snf0;
2137 if (Tof_correc_){
2138 double t = tanl();
2139 double dphi(0);
2140
2141 Helix work = *(Helix*)this;
2142 work.ignoreErrorMatrix();
2143 double phi_ip(0);
2144 Hep3Vector ip(0, 0, 0);
2145 work.pivot(ip);
2146 phi_ip = work.phi0();
2147 if (fabs(phi - phi_ip) > M_PI) {
2148 if (phi > phi_ip) phi -= 2 * M_PI;
2149 else phi_ip -= 2 * M_PI;
2150 }
2151 dphi = phi - phi_ip;
2152
2153 double l = fabs(radius() * dphi * sqrt(1 + t * t));
2154 double pmag( sqrt( 1.0 + t*t ) / kappa());
2155 double mass_over_p( mass_ / pmag );
2156 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2157 tofest = l / ( 29.9792458 * beta );
2158 if(csmflag==1 && (*HitMdc).wire().y()>0.) tofest= -1. * tofest;
2159 }
2160
2161 const HepSymMatrix& ea = Ea();
2162 HelixSeg.Ea_pre_fwd(ea);
2163 HelixSeg.Ea_filt_fwd(ea);
2164
2165
2166 const HepVector& v_a = a();
2167 HelixSeg.a_pre_fwd(v_a);
2168 HelixSeg.a_filt_fwd(v_a);
2169
2170 /*
2171
2172 HepPoint3D pivot_work = pivot();
2173 HepVector a_work = a();
2174 HepSymMatrix ea_work = Ea();
2175 KalFitTrack helix_work(pivot_work, a_work, ea_work, 0, 0, 0);
2176 helix_work.pivot(ip);
2177
2178 HepVector a_temp = helix_work.a();
2179 HepSymMatrix ea_temp = helix_work.Ea();
2180
2181 HelixSeg.Ea_pre_fwd(ea_temp);
2182 HelixSeg.a_pre_fwd(a_temp);
2183
2184 // My God, for the protection purpose
2185 HelixSeg.Ea_filt_fwd(ea_temp);
2186 HelixSeg.a_filt_fwd(a_temp);
2187
2188 */
2189
2190 double dchi2R(99999999), dchi2L(99999999);
2191
2192 HepVector v_H(5, 0);
2193 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
2194 v_H[3] = -meas.z();
2195 HepMatrix v_HT = v_H.T();
2196
2197 double estim = (v_HT * v_a)[0];
2198 HepVector ea_v_H = ea * v_H;
2199 HepMatrix ea_v_HT = (ea_v_H).T();
2200 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
2201
2202 HepSymMatrix eaNewL(5), eaNewR(5);
2203 HepVector aNewL(5), aNewR(5);
2204#ifdef YDEBUG
2205 cout<<"phi "<<phi<<" snf0 "<<snf0<<" csf0 "<<csf0<<endl
2206 <<" meas: "<<meas <<endl;
2207 cout<<"the related matrix:"<<endl;
2208 cout<<"v_H "<<v_H<<endl;
2209 cout<<"v_a "<<v_a<<endl;
2210 cout<<"ea "<<ea<<endl;
2211 cout<<"v_H_T_ea_v_H "<<v_H_T_ea_v_H<<endl;
2212 cout<<"LR_= "<< LR_ << "lr= " << lr <<endl;
2213#endif
2214
2215 double time = (*HitMdc).tdc();
2216 // if (Tof_correc_)
2217 // time = time - way*tofest;
2218
2219 // double dd = 1.0e-4 * 40.0 * time;
2220 //the following 3 line : tempory
2221
2222 double drifttime = getDriftTime(*HitMdc , tofest);
2223 double ddl = getDriftDist(*HitMdc, drifttime, 0 );
2224 double ddr = getDriftDist(*HitMdc, drifttime, 1 );
2225 double erddl = getSigma( *HitMdc, a()[4], 0, ddl);
2226 double erddr = getSigma( *HitMdc, a()[4], 1, ddr);
2227
2228 if(debug_ == 4){
2229 std::cout<<"the pure drift time is "<<drifttime<<std::endl;
2230 std::cout<<"the drift dist left hypothesis is "<<ddl<<std::endl;
2231 std::cout<<"the drift dist right hypothesis is "<<ddr<<std::endl;
2232 std::cout<<"the error sigma left hypothesis is "<<erddl<<std::endl;
2233 std::cout<<"the error sigma right hypothesis is "<<erddr<<std::endl;
2234 }
2235 double dmeas2[2] = { 0.1*ddl , 0.1*ddr }; //millimeter to centimeter
2236 double er_dmeas2[2] = {0. , 0.} ;
2237
2238 if(resolflag_ == 1) {
2239 er_dmeas2[0] = 0.1*erddl;
2240 er_dmeas2[1] = 0.1*erddr;
2241 }
2242 else if(resolflag_ == 0) {
2243 // int layid = (*HitMdc).wire().layer().layerId();
2244 // double sigma = getSigma(layid, dd);
2245 // er_dmeas2[0] = er_dmeas2[1] = sigma;
2246 }
2247
2248 // Left hypothesis :
2249 if ((LR_==0 && lr != 1.0) ||
2250 (LR_==1 && lr == -1.0)) {
2251
2252 double er_dmeasL, dmeasL;
2253 if(Tof_correc_) {
2254 dmeasL = (-1.0)*fabs(dmeas2[0]);
2255 er_dmeasL = er_dmeas2[0];
2256 } else {
2257 dmeasL = (-1.0)*fabs((*HitMdc).dist()[0]);
2258 er_dmeasL = (*HitMdc).erdist()[0];
2259 }
2260
2261 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
2262 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
2263 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
2264 if(0. == RkL) RkL = 1.e-4;
2265
2266 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
2267
2268 aNewL = v_a + diffL;
2269 double sigL = dmeasL -(v_HT * aNewL)[0];
2270 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
2271 }
2272
2273 if ((LR_==0 && lr != -1.0) ||
2274 (LR_==1 && lr == 1.0)) {
2275
2276 double er_dmeasR, dmeasR;
2277 if(Tof_correc_) {
2278 dmeasR = dmeas2[1];
2279 er_dmeasR = er_dmeas2[1];
2280 } else {
2281 dmeasR = fabs((*HitMdc).dist()[1]);
2282 er_dmeasR = (*HitMdc).erdist()[1];
2283 }
2284
2285 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
2286 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
2287 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
2288 if(0. == RkR) RkR = 1.e-4;
2289
2290 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
2291
2292 aNewR = v_a + diffR;
2293 double sigR = dmeasR -(v_HT * aNewR)[0];
2294 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
2295 }
2296
2297 // Update Kalman result :
2298 double dchi2_loc(0);
2299
2300#ifdef YDEBUG
2301 cout<<" track::update_hits......"<<endl;
2302 std::cout << " dchi2R = " << dchi2R << ", dchi2L = " << dchi2L
2303 << " estimate = "<<estim<< std::endl;
2304 std::cout << " dmeasR = " << dmeas2[1]
2305 << ", dmeasL = " << (-1.0)*fabs(dmeas2[0]) << ", check HitMdc.ddl = "
2306 << (*HitMdc).dist()[0] << std::endl;
2307 std::cout<<"dchi2L : "<<dchi2L<<" ,dchi2R: "<<dchi2R<<std::endl;
2308#endif
2309
2310
2311 if (fabs(dchi2R-dchi2L)<10. && inext>0) {
2312
2313 KalFitHitMdc & HitMdc_next = HitsMdc_[inext];
2314
2315 Helix H_fromR(pivot(), aNewR, eaNewR);
2316 double chi2_fromR(chi2_next(H_fromR, HitMdc_next,csmflag));
2317 //double chi2_fromR(chi2_next(H_fromR, HitMdc_next));
2318
2319 Helix H_fromL(pivot(), aNewL, eaNewL);
2320 double chi2_fromL(chi2_next(H_fromL, HitMdc_next,csmflag));
2321 //double chi2_fromL(chi2_next(H_fromL, HitMdc_next));
2322#ifdef YDEBUG
2323 std::cout << " chi2_fromL = " << chi2_fromL
2324 << ", chi2_fromR = " << chi2_fromR << std::endl;
2325#endif
2326 if (fabs(chi2_fromR-chi2_fromL)<10.){
2327 int inext2(-1);
2328 if (inext+1<HitsMdc_.size())
2329 for( int k=inext+1 ; k < HitsMdc_.size(); k++ )
2330 if (!(HitsMdc_[k].chi2()<0)) {
2331 inext2 = k;
2332 break;
2333 }
2334
2335 if (inext2 != -1){
2336 KalFitHitMdc & HitMdc_next2 = HitsMdc_[inext2];
2337 // double chi2_fromR2(chi2_next(H_fromR, HitMdc_next2));
2338 // double chi2_fromL2(chi2_next(H_fromL, HitMdc_next2));
2339 double chi2_fromR2(chi2_next(H_fromR, HitMdc_next2, csmflag));
2340 double chi2_fromL2(chi2_next(H_fromL, HitMdc_next2, csmflag));
2341#ifdef YDEBUG
2342 std::cout << " chi2_fromL2 = " << chi2_fromL2
2343 << ", chi2_fromR2 = " << chi2_fromR2 << std::endl;
2344#endif
2345 if (fabs(dchi2R+chi2_fromR+chi2_fromR2-(dchi2L+chi2_fromL+chi2_fromL2))<2) {
2346 if (chi2_fromR2<chi2_fromL2)
2347 dchi2L = DBL_MAX;
2348 else
2349 dchi2R = DBL_MAX;
2350 }
2351 }
2352 }
2353
2354 if (!(dchi2L == DBL_MAX && dchi2R == DBL_MAX)) {
2355 if (dchi2R+chi2_fromR < dchi2L+chi2_fromL){
2356 dchi2L = DBL_MAX;
2357#ifdef YDEBUG
2358 std::cout << " choose right..." << std::endl;
2359#endif
2360 } else {
2361 dchi2R = DBL_MAX;
2362#ifdef YDEBUG
2363 std::cout << " choose left..." << std::endl;
2364#endif
2365 }
2366 }
2367 }
2368
2369 if ((0 < dchi2R && dchi2R < dchi2cutf_calib[layerid]) ||
2370 (0 < dchi2L && dchi2L < dchi2cutf_calib[layerid])) {
2371
2372 if (((LR_==0 && dchi2R < dchi2L && lr != -1) ||
2373 (LR_==1 && lr == 1)) &&
2374 fabs(aNewR[2]-a()[2]) < 1000. && aNewR[2]) {
2375 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2R < dchi2cutf_calib[layerid]){
2376 nchits_++;
2377 if (flag_ster) nster_++;
2378 if(layerid>19){
2379 Ea(eaNewR);
2380 HelixSeg.Ea_filt_fwd(eaNewR);
2381 a(aNewR);
2382 HelixSeg.a_filt_fwd(aNewR);
2383 }
2384 /*
2385 Ea(eaNewR);
2386 a(aNewR);
2387
2388 KalFitTrack helixR(pivot_work, aNewR, eaNewR, 0, 0, 0);
2389 helixR.pivot(ip);
2390
2391 a_temp = helixR.a();
2392 ea_temp = helixR.Ea();
2393
2394 HelixSeg.Ea_filt_fwd(ea_temp);
2395 HelixSeg.a_filt_fwd(a_temp);
2396 //HelixSeg.filt_include(1);
2397
2398 */
2399
2400 chiSq_ += dchi2R;
2401 dchi2_loc = dchi2R;
2402 if (l_mass_ == lead_){
2403 (*HitMdc).LR(1);
2404 HelixSeg.LR(1);
2405 }
2406 update_bit((*HitMdc).wire().layer().layerId());
2407 }
2408 } else if (((LR_==0 && dchi2L <= dchi2R && lr != 1) ||
2409 (LR_==1 && lr == -1)) &&
2410 fabs(aNewL[2]-a()[2]) < 1000. && aNewL[2]){
2411 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2L < dchi2cutf_calib[layerid]){
2412 nchits_++;
2413 if (flag_ster) nster_++;
2414 if(layerid>19){
2415 Ea(eaNewL);
2416 HelixSeg.Ea_filt_fwd(eaNewL);
2417 a(aNewL);
2418 HelixSeg.a_filt_fwd(aNewL);
2419 }
2420
2421 /*
2422 Ea(eaNewL);
2423 a(aNewL);
2424
2425 KalFitTrack helixL(pivot_work, aNewL, eaNewL, 0, 0, 0);
2426 helixL.pivot(ip);
2427 a_temp = helixL.a();
2428 ea_temp = helixL.Ea();
2429 HelixSeg.Ea_filt_fwd(ea_temp);
2430 HelixSeg.a_filt_fwd(a_temp);
2431 //HelixSeg.filt_include(1);
2432
2433 */
2434
2435 chiSq_ += dchi2L;
2436 dchi2_loc = dchi2L;
2437 if (l_mass_ == lead_){
2438 (*HitMdc).LR(-1);
2439 HelixSeg.LR(-1);
2440 }
2441 update_bit((*HitMdc).wire().layer().layerId());
2442 }
2443 }
2444 }
2445
2446 if (dchi2_loc > dchi2_max_) {
2447 dchi2_max_ = dchi2_loc ;
2448 r_max_ = pivot().perp();
2449 }
2450 dchi2out = dchi2_loc;
2451 // if(dchi2out == 0) {
2452 // dchi2out = ( (dchi2L < dchi2R ) ? dchi2L : dchi2R ) ;
2453 // }
2454 return chiSq_;
2455}
KalFitTrack::chi2_next
double chi2_next(Helix &H, KalFitHitMdc &HitMdc, int csmflag)
Definition: KalFitTrack.cxx:2950
KalFitTrack::update_bit
void update_bit(int i)
Definition: KalFitTrack.cxx:1807
KalFitTrack::nmdc_hit2_
static int nmdc_hit2_
Cut chi2 for each hit.
Definition: KalFitTrack.h:312
KalFitTrack::HelixSeg
KalFitHelixSeg & HelixSeg(int i)
Definition: KalFitTrack.h:235
KalFitTrack::dchi2cutf_calib
static double dchi2cutf_calib[43]
Definition: KalFitTrack.h:292
KalFitWire::stereo
unsigned int stereo(void) const
Definition: KalFitWire.h:75

Referenced by KalFitAlg::filter_fwd_calib().

◆ update_last()

void KalFitTrack::update_last ( void  )

Record the current parameters as ..._last information :

Definition at line 109 of file KalFitTrack.cxx.

110{
111 pivot_last_ = pivot();
112 a_last_ = a();
113 Ea_last_ = Ea();
114}

Referenced by filter(), and KalFitTrack().

◆ useLayer()

void KalFitTrack::useLayer ( int  iLay)
inline

Definition at line 347 of file KalFitTrack.h.

347 {
348 if(iLay>=0 && iLay<=43) myLayerUsed[iLay]=1;
349 }

Member Data Documentation

◆ Bznom_

double KalFitTrack::Bznom_ = 10
static

Definition at line 299 of file KalFitTrack.h.

Referenced by KalFitAlg::beginRun(), KalFitAlg::execute(), KalFitAlg::initialize(), KalFitTrack(), KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew(), pivot_numf(), and radius_numf().

◆ chi2_hitf_

double KalFitTrack::chi2_hitf_ = 1000
static

Cut chi2 for each hit.

Definition at line 282 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize().

◆ chi2_hits_

double KalFitTrack::chi2_hits_ = 1000
static

Definition at line 282 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize().

◆ chi2mdc_hit2_

double KalFitTrack::chi2mdc_hit2_
static

Definition at line 313 of file KalFitTrack.h.

◆ dchi2cutf_anal

double KalFitTrack::dchi2cutf_anal = {0.}
static

Definition at line 288 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut().

◆ dchi2cutf_calib

double KalFitTrack::dchi2cutf_calib = {0.}
static

Definition at line 292 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut(), and update_hits_csmalign().

◆ dchi2cuts_anal

double KalFitTrack::dchi2cuts_anal = {0.}
static

Definition at line 290 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut().

◆ dchi2cuts_calib

double KalFitTrack::dchi2cuts_calib = {0.}
static

Definition at line 294 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut(), and smoother_Mdc_csmalign().

◆ debug_

int KalFitTrack::debug_ = 0
static

for debug

Definition at line 280 of file KalFitTrack.h.

Referenced by getDriftDist(), getDriftTime(), getT0(), KalFitAlg::initialize(), setT0(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ drifttime_choice_

int KalFitTrack::drifttime_choice_ = 0
static

the drifttime choice

Definition at line 321 of file KalFitTrack.h.

Referenced by getDriftTime(), and KalFitAlg::initialize().

◆ factor_strag_

double KalFitTrack::factor_strag_ = 0.4
static

factor of energy loss straggling for electron

Definition at line 310 of file KalFitTrack.h.

Referenced by eloss(), and KalFitAlg::initialize().

◆ inner_steps_

int KalFitTrack::inner_steps_ = 3
static

Definition at line 285 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize(), and pivot_numf().

◆ LR_

int KalFitTrack::LR_ = 1
static

Use L/R decision from MdcRecHit information :

Definition at line 319 of file KalFitTrack.h.

Referenced by chi2_next(), KalFitAlg::initialize(), KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew(), LR(), and update_hits_csmalign().

◆ mdcGasRadlen_

double KalFitTrack::mdcGasRadlen_ = 0.
static

Definition at line 274 of file KalFitTrack.h.

Referenced by msgasmdc(), and KalFitAlg::setBesFromGdml().

◆ nmdc_hit2_

int KalFitTrack::nmdc_hit2_ = 500
static

Cut chi2 for each hit.

Definition at line 312 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize(), and update_hits_csmalign().

◆ numf_

int KalFitTrack::numf_ = 0
static

Flag for treatment of non-uniform mag field.

Definition at line 284 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize(), numf(), pivot_numf(), and radius_numf().

◆ numfcor_

int KalFitTrack::numfcor_ = 1
static

NUMF treatment improved.

Definition at line 297 of file KalFitTrack.h.

Referenced by KalFitAlg::beginRun(), KalFitAlg::execute(), KalFitAlg::initialize(), and pivot_numf().

◆ outer_steps_

int KalFitTrack::outer_steps_ = 3
static

Definition at line 286 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize(), and pivot_numf().

◆ resolflag_

int KalFitTrack::resolflag_ = 0
static

wire resoltion flag

Definition at line 316 of file KalFitTrack.h.

Referenced by chi2_next(), resol(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ steplev_

int KalFitTrack::steplev_ = 0
static

Level of precision (= 1 : 5steps for all tracks; = 2: 5 steps for very non uniform part)

Definition at line 302 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize().

◆ strag_

int KalFitTrack::strag_ = 1
static

Flag to take account of energy loss straggling :

Definition at line 308 of file KalFitTrack.h.

Referenced by eloss(), and KalFitAlg::initialize().

◆ Tof_correc_

int KalFitTrack::Tof_correc_ = 1
static

Flag for TOF correction.

Definition at line 305 of file KalFitTrack.h.

Referenced by chi2_next(), KalFitAlg::initialize(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ tofall_

int KalFitTrack::tofall_ = 1
static

Definition at line 314 of file KalFitTrack.h.

Referenced by KalFitAlg::complete_track(), eloss(), KalFitAlg::fillTds_back(), KalFitAlg::initialize(), and tof().

◆ tprop_

int KalFitTrack::tprop_ = 1
static

for signal propagation correction

Definition at line 277 of file KalFitTrack.h.

Referenced by getDriftTime(), and KalFitAlg::initialize().

The documentation for this class was generated from the following files: