#include <KalFitTrack.h>
Inheritance diagram for KalFitTrack:
Public Member Functions | |
| const HepVector & | a (const HepVector &newA) |
| sets helix parameters. | |
| const HepVector & | a (void) const |
| returns helix parameters. | |
| const HepVector & | a (const HepVector &newA) |
| sets helix parameters. | |
| const HepVector & | a (void) const |
| returns helix parameters. | |
| const CLHEP::HepVector & | a_forMdc (void) const |
| const CLHEP::HepVector & | a_forMdc (void) const |
| const CLHEP::HepVector & | a_last (void) const |
| const CLHEP::HepVector & | a_last (void) const |
| void | add_nhit_r (void) |
| void | add_nhit_r (void) |
| void | add_nhit_z (void) |
| void | add_nhit_z (void) |
| void | addPathSM (double path) |
| void | addPathSM (double path) |
| void | addTofSM (double time) |
| void | addTofSM (double time) |
| double | alpha (void) const |
| double | alpha (void) const |
| void | appendHelixSegs (KalFitHelixSeg s) |
| void | appendHelixSegs (KalFitHelixSeg s) |
| void | appendHitsMdc (KalFitHitMdc h) |
| Functions for Mdc hits list. | |
| void | appendHitsMdc (KalFitHitMdc h) |
| Functions for Mdc hits list. | |
| double | approach (HepPoint3D pfwd, HepPoint3D pbwd, bool doSagCorrection) const |
| double | approach (KalFitHitMdc &hit, bool doSagCorrection) const |
| double | approach (HepPoint3D pfwd, HepPoint3D pbwd, bool doSagCorrection) const |
| double | approach (KalFitHitMdc &hit, bool doSagCorrection) const |
| double | bFieldZ (void) const |
| double | bFieldZ (double) |
| sets/returns z componet of the magnetic field. | |
| double | bFieldZ (void) const |
| double | bFieldZ (double) |
| sets/returns z componet of the magnetic field. | |
| const HepPoint3D & | center (void) const |
| returns position of helix center(z = 0.); | |
| const HepPoint3D & | center (void) const |
| returns position of helix center(z = 0.); | |
| void | chgmass (int i) |
| void | chgmass (int i) |
| double | chi2_next (Helix &H, KalFitHitMdc &HitMdc) |
| double | chi2_next (Helix &H, KalFitHitMdc &HitMdc, int csmflag) |
| double | chi2_next (Helix &H, KalFitHitMdc &HitMdc) |
| double | chi2_next (Helix &H, KalFitHitMdc &HitMdc, int csmflag) |
| void | chiSq (double c) |
| double | chiSq (void) const |
| void | chiSq (double c) |
| double | chiSq (void) const |
| void | chiSq_back (double c) |
| double | chiSq_back (void) const |
| void | chiSq_back (double c) |
| double | chiSq_back (void) const |
| double | cor_tanldep (double *p, double er) |
| Correct the error according the current tanl value :. | |
| double | cor_tanldep (double *p, double er) |
| Correct the error according the current tanl value :. | |
| double | cosPhi0 (void) const |
| double | cosPhi0 (void) const |
| double | curv (void) const |
| double | curv (void) const |
| double | dchi2_max (void) const |
| double | dchi2_max (void) const |
| HepMatrix | del4MDelA (double phi, double mass) const |
| HepMatrix | del4MDelA (double phi, double mass) const |
| HepMatrix | del4MXDelA (double phi, double mass) const |
| HepMatrix | del4MXDelA (double phi, double mass) const |
| HepMatrix | delApDelA (const HepVector &ap) const |
| HepMatrix | delApDelA (const HepVector &ap) const |
| HepMatrix | delMDelA (double phi) const |
| HepMatrix | delMDelA (double phi) const |
| HepMatrix | delXDelA (double phi) const |
| HepMatrix | delXDelA (double phi) const |
| Hep3Vector | direction (double dPhi=0.) const |
| returns direction vector after rotating angle dPhi in phi direction. | |
| Hep3Vector | direction (double dPhi=0.) const |
| returns direction vector after rotating angle dPhi in phi direction. | |
| double | dr (void) const |
| returns an element of parameters. | |
| double | dr (void) const |
| returns an element of parameters. | |
| double | dz (void) const |
| double | dz (void) const |
| const HepSymMatrix & | Ea (const HepSymMatrix &newdA) |
| sets helix paramters and error matrix. | |
| const HepSymMatrix & | Ea (void) const |
| returns error matrix. | |
| const HepSymMatrix & | Ea (const HepSymMatrix &newdA) |
| sets helix paramters and error matrix. | |
| const HepSymMatrix & | Ea (void) const |
| returns error matrix. | |
| const CLHEP::HepSymMatrix & | Ea_forMdc (void) const |
| const CLHEP::HepSymMatrix & | Ea_forMdc (void) const |
| const CLHEP::HepSymMatrix & | Ea_last (void) const |
| const CLHEP::HepSymMatrix & | Ea_last (void) const |
| void | eloss (double path, const KalFitMaterial &m, int index) |
| Calculate total energy lost in material. | |
| void | eloss (double path, const KalFitMaterial &m, int index) |
| Calculate total energy lost in material. | |
| double | filter (double v_m, const CLHEP::HepVector &m_H, double v_d, double m_V) |
| double | filter (double v_m, const CLHEP::HepVector &m_H, double v_d, double m_V) |
| void | fiTerm (double fi) |
| void | fiTerm (double fi) |
| double | getDigi () const |
| double | getDigi () const |
| double | getDriftDist (KalFitHitMdc &hitmdc, double drifttime, int lr) const |
| double | getDriftDist (KalFitHitMdc &hitmdc, double drifttime, int lr) const |
| double | getDriftTime (KalFitHitMdc &hitmdc, double toftime) const |
| double | getDriftTime (KalFitHitMdc &hitmdc, double toftime) const |
| double | getFiTerm (void) |
| double | getFiTerm (void) |
| HepSymMatrix | getInitMatrix (void) const |
| HepSymMatrix | getInitMatrix (void) const |
| double | getPathSM (void) |
| double | getPathSM (void) |
| double | getSigma (KalFitHitMdc &hitmdc, double tanlam, int lr, double dist) const |
| double | getSigma (int layerId, double driftDist) const |
| double | getSigma (KalFitHitMdc &hitmdc, double tanlam, int lr, double dist) const |
| double | getSigma (int layerId, double driftDist) const |
| double | getT0 (void) const |
| double | getT0 (void) const |
| double | getTofSM (void) |
| double | getTofSM (void) |
| KalFitHelixSeg & | HelixSeg (int i) |
| KalFitHelixSeg & | HelixSeg (int i) |
| vector< KalFitHelixSeg > & | HelixSegs (void) |
| void | HelixSegs (vector< KalFitHelixSeg > &vs) |
| vector< KalFitHelixSeg > & | HelixSegs (void) |
| void | HelixSegs (vector< KalFitHelixSeg > &vs) |
| KalFitHitMdc & | HitMdc (int i) |
| KalFitHitMdc & | HitMdc (int i) |
| vector< KalFitHitMdc > & | HitsMdc (void) |
| void | HitsMdc (vector< KalFitHitMdc > &lh) |
| vector< KalFitHitMdc > & | HitsMdc (void) |
| void | HitsMdc (vector< KalFitHitMdc > &lh) |
| 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. | |
| 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. | |
| void | insist (int t) |
| int | insist (void) const |
| Extractor :. | |
| void | insist (int t) |
| int | insist (void) const |
| Extractor :. | |
| double | intersect_cylinder (double r) const |
| Intersection with different geometry. | |
| double | intersect_cylinder (double r) const |
| Intersection with different geometry. | |
| double | intersect_xy_plane (double z) const |
| double | intersect_xy_plane (double z) const |
| double | intersect_yz_plane (const HepTransform3D &plane, double x) const |
| double | intersect_yz_plane (const HepTransform3D &plane, double x) const |
| double | intersect_zx_plane (const HepTransform3D &plane, double y) const |
| double | intersect_zx_plane (const HepTransform3D &plane, double y) const |
| KalFitTrack (const HepPoint3D &pivot, const CLHEP::HepVector &a, const CLHEP::HepSymMatrix &Ea, unsigned int m, double chiSq, unsigned int nhits) | |
| constructor | |
| KalFitTrack (const HepPoint3D &pivot, const CLHEP::HepVector &a, const CLHEP::HepSymMatrix &Ea, unsigned int m, double chiSq, unsigned int nhits) | |
| constructor | |
| double | kappa (void) const |
| double | kappa (void) const |
| double | mass (void) const |
| double | mass (void) const |
| CLHEP::Hep3Vector * | mom (void) |
| CLHEP::Hep3Vector * | mom (void) |
| HepLorentzVector | momentum (double dPhi, double mass, HepPoint3D &x, HepSymMatrix &Emx) 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) const |
| returns 4momentum 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. | |
| Hep3Vector | momentum (double dPhi=0.) const |
| returns momentum 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. | |
| 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) const |
| returns 4momentum 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. | |
| Hep3Vector | momentum (double dPhi=0.) const |
| returns momentum vector after rotating angle dPhi in phi direction. | |
| void | ms (double path, const KalFitMaterial &m, int index) |
| void | ms (double path, const KalFitMaterial &m, int index) |
| void | msgasmdc (double path, int index) |
| Calculate multiple scattering angle. | |
| void | msgasmdc (double path, int index) |
| Calculate multiple scattering angle. | |
| unsigned int | ncath (void) const |
| unsigned int | ncath (void) const |
| void | nchits (int n) |
| unsigned int | nchits (void) const |
| void | nchits (int n) |
| unsigned int | nchits (void) const |
| void | ndf_back (int n) |
| int | ndf_back (void) const |
| void | ndf_back (int n) |
| int | ndf_back (void) const |
| int | nhit_r (void) const |
| int | nhit_r (void) const |
| int | nhit_z (void) const |
| int | nhit_z (void) const |
| void | nster (int n) |
| unsigned int | nster (void) const |
| void | nster (int n) |
| unsigned int | nster (void) const |
| void | number_wirhit (void) |
| void | number_wirhit (void) |
| void | order_hits (void) |
| void | order_hits (void) |
| void | order_wirhit (int index) |
| void | order_wirhit (int index) |
| void | p_kaon (double pl) |
| double | p_kaon (void) const |
| void | p_kaon (double pl) |
| double | p_kaon (void) const |
| void | p_proton (double pl) |
| double | p_proton (void) const |
| void | p_proton (double pl) |
| double | p_proton (void) const |
| int | pat1 (void) const |
| int | pat1 (void) const |
| int | pat2 (void) const |
| int | pat2 (void) const |
| double | path_ab (void) const |
| double | path_ab (void) const |
| void | path_add (double path) |
| Update the path length estimation. | |
| void | path_add (double path) |
| Update the path length estimation. | |
| double | path_rd (void) const |
| double | path_rd (void) const |
| void | pathip (double pl) |
| double | pathip (void) const |
| void | pathip (double pl) |
| double | pathip (void) const |
| double | PathL (int layer) |
| Function to calculate the path length in the layer. | |
| double * | pathl (void) |
| double | PathL (int layer) |
| Function to calculate the path length in the layer. | |
| double * | pathl (void) |
| double | phi0 (void) const |
| double | phi0 (void) const |
| const HepPoint3D & | pivot (const HepPoint3D &newPivot) |
| sets pivot position. | |
| const HepPoint3D & | pivot (void) const |
| returns pivot position. | |
| const HepPoint3D & | pivot (const HepPoint3D &newPivot) |
| sets pivot position. | |
| const HepPoint3D & | pivot (void) const |
| returns pivot position. | |
| const HepPoint3D & | pivot_forMdc (void) const |
| const HepPoint3D & | pivot_forMdc (void) const |
| const HepPoint3D & | pivot_last (void) const |
| returns helix parameters | |
| const HepPoint3D & | pivot_last (void) const |
| returns helix parameters | |
| const HepPoint3D & | pivot_numf (const HepPoint3D &newPivot, double &pathl) |
| const HepPoint3D & | pivot_numf (const HepPoint3D &newPivot) |
| Sets pivot position in a given mag field. | |
| const HepPoint3D & | pivot_numf (const HepPoint3D &newPivot, double &pathl) |
| const HepPoint3D & | pivot_numf (const HepPoint3D &newPivot) |
| Sets pivot position in a given mag field. | |
| const HepPoint3D & | point_last (void) |
| void | point_last (const HepPoint3D &point) |
| set and give out the last point of the track | |
| const HepPoint3D & | point_last (void) |
| void | point_last (const HepPoint3D &point) |
| set and give out the last point of the track | |
| double | r0 (void) const |
| double | r0 (void) const |
| double | r_max (void) const |
| double | r_max (void) const |
| double | radius (void) const |
| returns radious of helix. | |
| double | radius (void) const |
| returns radious of helix. | |
| double | radius_numf (void) const |
| Estimation of the radius in a given mag field. | |
| double | radius_numf (void) const |
| Estimation of the radius in a given mag field. | |
| void | set (const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea) |
| sets helix pivot position, parameters, and error matrix. | |
| void | set (const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea) |
| sets helix pivot position, parameters, and error matrix. | |
| double | sinPhi0 (void) const |
| double | sinPhi0 (void) const |
| double | smoother_Mdc (KalFitHitMdc &HitMdc, CLHEP::Hep3Vector &meas, KalFitHelixSeg &seg, double &dchi2, int csmflag) |
| double | smoother_Mdc (KalFitHelixSeg &seg, CLHEP::Hep3Vector &meas, int &flg, int csmflag) |
| Kalman smoother for Mdc. | |
| double | smoother_Mdc (KalFitHitMdc &HitMdc, CLHEP::Hep3Vector &meas, KalFitHelixSeg &seg, double &dchi2, int csmflag) |
| 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_csmalign (KalFitHelixSeg &seg, CLHEP::Hep3Vector &meas, int &flg, int csmflag) |
| double | tanl (void) const |
| double | tanl (void) const |
| double | tof (void) const |
| void | tof (double path) |
| Update the tof estimation. | |
| double | tof (void) const |
| void | tof (double path) |
| Update the tof estimation. | |
| double | tof_kaon (void) const |
| double | tof_kaon (void) const |
| double | tof_proton (void) const |
| double | tof_proton (void) const |
| void | trasan_id (int t) |
| int | trasan_id (void) const |
| void | trasan_id (int t) |
| int | trasan_id (void) const |
| void | type (int t) |
| Reinitialize (modificator). | |
| int | type (void) const |
| void | type (int t) |
| Reinitialize (modificator). | |
| int | type (void) const |
| void | update_bit (int i) |
| void | update_bit (int i) |
| void | update_forMdc (void) |
| void | update_forMdc (void) |
| double | update_hits (KalFitHelixSeg &HelixSeg, int inext, CLHEP::Hep3Vector &meas, int way, 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 (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_csmalign (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) |
| void | update_last (void) |
| Record the current parameters as ..._last information :. | |
| void | update_last (void) |
| Record the current parameters as ..._last information :. | |
| HepPoint3D | x (double dPhi, HepSymMatrix &Ex) const |
| returns position and convariance matrix(Ex) after rotation. | |
| double * | x (double dPhi, double p[3]) const |
| HepPoint3D | x (double dPhi=0.) const |
| returns position after rotating angle dPhi in phi direction. | |
| HepPoint3D | x (double dPhi, HepSymMatrix &Ex) const |
| returns position and convariance matrix(Ex) after rotation. | |
| double * | x (double dPhi, double p[3]) const |
| HepPoint3D | x (double dPhi=0.) const |
| returns position after rotating angle dPhi in phi direction. | |
| ~KalFitTrack (void) | |
| destructor | |
| ~KalFitTrack (void) | |
| destructor | |
Static Public Member Functions | |
| void | back (int i) |
| int | back (void) |
| void | back (int i) |
| int | back (void) |
| void | lead (int i) |
| int | lead (void) |
| Magnetic field map. | |
| void | lead (int i) |
| int | lead (void) |
| Magnetic field map. | |
| void | LR (int x) |
| void | LR (int x) |
| double | mass (int i) |
| double | mass (int i) |
| int | nmass (void) |
| int | nmass (void) |
| void | numf (int i) |
| int | numf (void) |
| void | numf (int i) |
| int | numf (void) |
| void | resol (int i) |
| int | resol (void) |
| void | resol (int i) |
| int | resol (void) |
| void | setIMdcGeomSvc (IMdcGeomSvc *igeomsvc) |
| void | setIMdcGeomSvc (IMdcGeomSvc *igeomsvc) |
| void | setInitMatrix (HepSymMatrix m) |
| void | setInitMatrix (HepSymMatrix m) |
| void | setMagneticFieldSvc (void) |
| void | setMagneticFieldSvc (void) |
| void | setMdcCalibFunSvc (const MdcCalibFunSvc *calibsvc) |
| void | setMdcCalibFunSvc (const MdcCalibFunSvc *calibsvc) |
| void | setMdcDigiCol (MdcDigiCol *digicol) |
| void | setMdcDigiCol (MdcDigiCol *digicol) |
| void | setT0 (double t0) |
| void | setT0 (double t0) |
Static Public Attributes | |
| double | Bznom_ = 10 |
| double | chi2_hitf_ = 1000 |
| Cut chi2 for each hit. | |
| double | chi2_hits_ = 1000 |
| Cut chi2 for each hit. | |
| double | chi2mdc_hit2_ |
| const double | ConstantAlpha = 333.564095 |
| Constant alpha for uniform field. | |
| double | dchi2cutf_anal [43] = {0.} |
| double | dchi2cutf_calib [43] = {0.} |
| double | dchi2cuts_anal [43] = {0.} |
| double | dchi2cuts_calib [43] = {0.} |
| int | debug_ = 0 |
| for debug | |
| int | drifttime_choice_ = 0 |
| the drifttime choice | |
| double | factor_strag_ = 0.4 |
| factor of energy loss straggling for electron | |
| int | inner_steps_ = 3 |
| int | LR_ = 1 |
| Use L/R decision from MdcRecHit information :. | |
| double | mdcGasRadlen_ = 0. |
| int | nmdc_hit2_ = 500 |
| Cut chi2 for each hit. | |
| int | numf_ = 0 |
| Flag for treatment of non-uniform mag field. | |
| int | numfcor_ = 1 |
| NUMF treatment improved. | |
| int | outer_steps_ = 3 |
| int | resolflag_ = 0 |
| wire resoltion flag | |
| int | steplev_ = 0 |
| int | strag_ = 1 |
| Flag to take account of energy loss straggling :. | |
| int | Tof_correc_ = 1 |
| Flag for TOF correction. | |
| int | tofall_ = 1 |
| int | tprop_ = 1 |
| for signal propagation correction | |
Private Types | |
| enum | { NMASS = 5 } |
| enum | { NMASS = 5 } |
Private Attributes | |
| CLHEP::HepVector | a_forMdc_ |
| CLHEP::HepVector | a_last_ |
| double | chiSq_ |
| double | chiSq_back_ |
| double | dchi2_max_ |
| CLHEP::HepSymMatrix | Ea_forMdc_ |
| CLHEP::HepSymMatrix | Ea_last_ |
| double | fiTerm_ |
| vector< KalFitHelixSeg > | HelixSegs_ |
| vector< KalFitHelixSeg > | HelixSegs_ |
| vector< KalFitHitMdc > | HitsMdc_ |
| vector< KalFitHitMdc > | HitsMdc_ |
| int | insist_ |
| int | l_mass_ |
| int | layer_prec_ |
| double | mass_ |
| CLHEP::Hep3Vector | mom_ [43] |
| unsigned int | ncath_ |
| unsigned int | nchits_ |
| unsigned int | ndf_back_ |
| int | nhit_r_ |
| int | nhit_z_ |
| unsigned int | nster_ |
| double | p_kaon_ |
| double | p_proton_ |
| int | pat1_ |
| int | pat2_ |
| double | path_ab_ |
| double | path_rd_ |
| double | pathip_ |
| double | PathL_ [43] |
| double | pathSM_ |
| HepPoint3D | pivot_forMdc_ |
| HepPoint3D | pivot_last_ |
| HepPoint3D | point_last_ |
| double | r0_ |
| double | r_max_ |
| double | tof2_ |
| double | tof_ |
| double | tof_kaon_ |
| double | tof_proton_ |
| int | trasan_id_ |
| int | type_ |
Static Private Attributes | |
| int | back_ = 1 |
| Flags. | |
| const MdcCalibFunSvc * | CalibFunSvc_ |
| const MdcCalibFunSvc * | CalibFunSvc_ = 0 |
| double | EventT0_ = 0. |
| IMdcGeomSvc * | iGeomSvc_ |
| IMdcGeomSvc * | iGeomSvc_ = 0 |
| HepSymMatrix | initMatrix_ |
| int | lead_ = 2 |
| Flags. | |
| const double | MASS [NMASS] |
| MdcDigiCol * | mdcDigiCol_ |
| MdcDigiCol * | mdcDigiCol_ = 0 |
| const IMagneticFieldSvc * | MFSvc_ |
| const IMagneticFieldSvc * | MFSvc_ = 0 |
|
|
00073 { NMASS = 5 };
|
|
|
00073 { NMASS = 5 };
|
|
||||||||||||||||||||||||||||
|
constructor
00085 : Helix(pivot, a, Ea), type_(0), insist_(0), chiSq_(0), 00086 nchits_(0), nster_(0), ncath_(0), 00087 ndf_back_(0), chiSq_back_(0), pathip_(0), 00088 path_rd_(0), path_ab_(0), tof_(0), dchi2_max_(0), r_max_(0), 00089 tof_kaon_(0), tof_proton_(0), pat1_(0), pat2_(0), layer_prec_(0), 00090 trasan_id_(0), nhit_r_(0), nhit_z_(0),pathSM_(0),tof2_(0) 00091 { 00092 memset(PathL_, 0, sizeof(PathL_)); 00093 l_mass_ = m; 00094 if (m < NMASS) mass_ = MASS[m]; 00095 else mass_ = MASS[2]; 00096 r0_ = fabs(center().perp() - fabs(radius())); 00097 bFieldZ(Bznom_); 00098 update_last(); 00099 }
|
|
|
destructor
00103 {
00104 // delete all objects
00105
00106 }
|
|
||||||||||||||||||||||||||||
|
constructor
|
|
|
destructor
|
|
|
sets helix parameters.
|
|
|
returns helix parameters.
|
|
|
sets helix parameters.
00266 {
00267 m_a = i;
00268 updateCache();
00269 return m_a;
00270 }
|
|
|
returns helix parameters.
00254 {
00255 return m_a;
00256 }
|
|
|
00151 { return a_forMdc_; }
|
|
|
00151 { return a_forMdc_; }
|
|
|
00147 { return a_last_; }
|
|
|
00147 { return a_last_; }
|
|
|
00220 { nhit_r_++;}
|
|
|
00220 { nhit_r_++;}
|
|
|
00221 { nhit_z_++;}
|
|
|
00221 { nhit_z_++;}
|
|
|
|
|
|
01284 {
01285 pathSM_ += path;
01286 }
|
|
|
|
|
|
01289 {
01290 tof2_ += time;
01291 }
|
|
|
|
|
|
00292 {
00293
00294 return m_alpha;
00295 }
|
|
|
|
|
|
|
|
|
Functions for Mdc hits list.
|
|
|
Functions for Mdc hits list.
01825 { HitsMdc_.push_back(h);}
|
|
||||||||||||||||
|
|
|
||||||||||||
|
|
|
||||||||||||||||
|
00209 {
00210 // ...Cal. dPhi to rotate...
00211 // const TMDCWire & w = * link.wire();
00212 HepPoint3D positionOnWire, positionOnTrack;
00213 HepPoint3D pv = pivot();
00214 HepVector Va = a();
00215 HepSymMatrix Ma = Ea();
00216
00217 Helix _helix(pv, Va ,Ma);
00218 Hep3Vector Wire;
00219 Wire[0] = (pfwd - pbwd).x();
00220 Wire[1] = (pfwd - pbwd).y();
00221 Wire[2] = (pfwd - pbwd).z();
00222 // xyPosition(), returns middle position of a wire. z componet is 0.
00223 // w.xyPosition(wp);
00224 double wp[3];
00225 wp[0] = 0.5*(pfwd + pbwd).x();
00226 wp[1] = 0.5*(pfwd + pbwd).y();
00227 wp[2] = 0.;
00228 double wb[3];
00229 // w.backwardPosition(wb);
00230 wb[0] = pbwd.x();
00231 wb[1] = pbwd.y();
00232 wb[2] = pbwd.z();
00233 double v[3];
00234 v[0] = Wire.unit().x();
00235 v[1] = Wire.unit().y();
00236 v[2] = Wire.unit().z();
00237 // std::cout<<"Wire.unit() is "<<Wire.unit()<<std::endl;
00238
00239 // ...Sag correction...
00240 /* if (doSagCorrection) {
00241 HepVector3D dir = w.direction();
00242 HepPoint3D xw(wp[0], wp[1], wp[2]);
00243 HepPoint3D wireBackwardPosition(wb[0], wb[1], wb[2]);
00244 w.wirePosition(link.positionOnTrack().z(),
00245 xw,
00246 wireBackwardPosition,
00247 dir);
00248 v[0] = dir.x();
00249 v[1] = dir.y();
00250 v[2] = dir.z();
00251 wp[0] = xw.x();
00252 wp[1] = xw.y();
00253 wp[2] = xw.z();
00254 wb[0] = wireBackwardPosition.x();
00255 wb[1] = wireBackwardPosition.y();
00256 wb[2] = wireBackwardPosition.z();
00257 }
00258 */
00259 // ...Cal. dPhi to rotate...
00260 const HepPoint3D & xc = _helix.center();
00261 double xt[3];
00262 _helix.x(0., xt);
00263 double x0 = - xc.x();
00264 double y0 = - xc.y();
00265 double x1 = wp[0] + x0;
00266 double y1 = wp[1] + y0;
00267 x0 += xt[0];
00268 y0 += xt[1];
00269 //std::cout<<" x0 is: "<<x0<<" y0 is: "<<y0<<std::endl;
00270 //std::cout<<" x1 is: "<<x1<<" y1 is: "<<y1<<std::endl;
00271 //std::cout<<" xt[0] is: "<<xt[0]<<" xt[1] is: "<<xt[1]<<std::endl;
00272
00273 double dPhi = atan2(x0 * y1 - y0 * x1, x0 * x1 + y0 * y1);
00274 //std::cout<<" x0 * y1 - y0 * x1 is: "<<(x0 * y1 - y0 * x1)<<std::endl;
00275 //std::cout<<" x0 * x1 + y0 * y1 is: "<<(x0 * x1 + y0 * y1)<<std::endl;
00276 //std::cout<<" before loop dPhi is "<<dPhi<<std::endl;
00277 //...Setup...
00278 double kappa = _helix.kappa();
00279 double phi0 = _helix.phi0();
00280
00281 //...Axial case...
00282 /* if (!w.stereo()) {
00283 positionOnTrack = _helix.x(dPhi);
00284 HepPoint3D x(wp[0], wp[1], wp[2]);
00285 x.setZ(positionOnTrack.z());
00286 positionOnWire = x;
00287 //link.dPhi(dPhi);
00288 std::cout<<" in axial wire : positionOnTrack is "<<positionOnTrack
00289 <<" positionOnWire is "<<positionOnWire<<std::endl;
00290 return (positionOnTrack - positionOnWire).mag();
00291 }
00292 */
00293 double firstdfdphi = 0.;
00294 static bool first = true;
00295 if (first) {
00296 // extern BelleTupleManager * BASF_Histogram;
00297 // BelleTupleManager * m = BASF_Histogram;
00298 // h_nTrial = m->histogram("TTrack::approach nTrial", 100, 0., 100.);
00299 }
00300 //#endif
00301
00302 //...Stereo case...
00303 double rho = Helix::ConstantAlpha / kappa;
00304 double tanLambda = _helix.tanl();
00305 static HepPoint3D x;
00306 double t_x[3];
00307 double t_dXdPhi[3];
00308 const double convergence = 1.0e-5;
00309 double l;
00310 unsigned nTrial = 0;
00311 while (nTrial < 100) {
00312
00313 // x = link.positionOnTrack(_helix->x(dPhi));
00314 positionOnTrack = _helix.x(dPhi);
00315 x = _helix.x(dPhi);
00316 t_x[0] = x[0];
00317 t_x[1] = x[1];
00318 t_x[2] = x[2];
00319
00320 l = v[0] * t_x[0] + v[1] * t_x[1] + v[2] * t_x[2]
00321 - v[0] * wb[0] - v[1] * wb[1] - v[2] * wb[2];
00322
00323 double rcosPhi = rho * cos(phi0 + dPhi);
00324 double rsinPhi = rho * sin(phi0 + dPhi);
00325 t_dXdPhi[0] = rsinPhi;
00326 t_dXdPhi[1] = - rcosPhi;
00327 t_dXdPhi[2] = - rho * tanLambda;
00328
00329 //...f = d(Distance) / d phi...
00330 double t_d2Xd2Phi[2];
00331 t_d2Xd2Phi[0] = rcosPhi;
00332 t_d2Xd2Phi[1] = rsinPhi;
00333
00334 //...iw new...
00335 double n[3];
00336 n[0] = t_x[0] - wb[0];
00337 n[1] = t_x[1] - wb[1];
00338 n[2] = t_x[2] - wb[2];
00339
00340 double a[3];
00341 a[0] = n[0] - l * v[0];
00342 a[1] = n[1] - l * v[1];
00343 a[2] = n[2] - l * v[2];
00344 double dfdphi = a[0] * t_dXdPhi[0]
00345 + a[1] * t_dXdPhi[1]
00346 + a[2] * t_dXdPhi[2];
00347
00348 if (nTrial == 0) {
00349 // break;
00350 firstdfdphi = dfdphi;
00351 }
00352
00353 //...Check bad case...
00354 if (nTrial > 3) {
00355 // std::cout<<" BAD CASE!!, calculate approach ntrial = "<<nTrial<< endl;
00356 }
00357 //...Is it converged?...
00358 if (fabs(dfdphi) < convergence)
00359 break;
00360
00361 double dv = v[0] * t_dXdPhi[0]
00362 + v[1] * t_dXdPhi[1]
00363 + v[2] * t_dXdPhi[2];
00364 double t0 = t_dXdPhi[0] * t_dXdPhi[0]
00365 + t_dXdPhi[1] * t_dXdPhi[1]
00366 + t_dXdPhi[2] * t_dXdPhi[2];
00367 double d2fd2phi = t0 - dv * dv
00368 + a[0] * t_d2Xd2Phi[0]
00369 + a[1] * t_d2Xd2Phi[1];
00370 // + a[2] * t_d2Xd2Phi[2];
00371
00372 dPhi -= dfdphi / d2fd2phi;
00373 ++nTrial;
00374 }
00375 //std::cout<<" dPhi is: "<<dPhi<<std::endl;
00376 //...Cal. positions...
00377 positionOnWire[0] = wb[0] + l * v[0];
00378 positionOnWire[1] = wb[1] + l * v[1];
00379 positionOnWire[2] = wb[2] + l * v[2];
00380
00381 //std::cout<<"wb[0] is: "<<wb[0]<<" l is: "<<l<<" v[0] is: "<<v[0]<<std::endl;
00382 //std::cout<<"wb[1] is: "<<wb[1]<<" v[1] is: "<<v[1]<<std::endl;
00383 //std::cout<<"wb[2] is: "<<wb[2]<<" v[2] is: "<<v[2]<<std::endl;
00384
00385 //std::cout<<" positionOnTrack is "<<positionOnTrack
00386 // <<" positionOnWire is "<<positionOnWire<<std::endl;
00387
00388 return (positionOnTrack - positionOnWire).mag();
00389 // link.dPhi(dPhi);
00390 // return nTrial;
00391 }
|
|
||||||||||||
|
00016 {
00017 //...Cal. dPhi to rotate...
00018 //const TMDCWire & w = * link.wire();
00019 HepPoint3D positionOnWire, positionOnTrack;
00020 HepPoint3D pv = pivot();
00021 //std::cout<<"the track pivot in approach is "<<pv<<std::endl;
00022 HepVector Va = a();
00023 //std::cout<<"the track parameters is "<<Va<<std::endl;
00024 HepSymMatrix Ma = Ea();
00025 //std::cout<<"the error matrix is "<<Ma<<std::endl;
00026
00027 Helix _helix(pv, Va ,Ma);
00028 //_helix.pivot(IP);
00029 const KalFitWire& w = hit.wire();
00030 Hep3Vector Wire = w.fwd() - w.bck();
00031 //xyPosition(), returns middle position of a wire. z componet is 0.
00032 //w.xyPosition(wp);
00033 double wp[3];
00034 wp[0] = 0.5*(w.fwd() + w.bck()).x();
00035 wp[1] = 0.5*(w.fwd() + w.bck()).y();
00036 wp[2] = 0.;
00037 double wb[3];
00038 //w.backwardPosition(wb);
00039 wb[0] = w.bck().x();
00040 wb[1] = w.bck().y();
00041 wb[2] = w.bck().z();
00042 double v[3];
00043 v[0] = Wire.unit().x();
00044 v[1] = Wire.unit().y();
00045 v[2] = Wire.unit().z();
00046 //std::cout<<"Wire.unit() is "<<Wire.unit()<<std::endl;
00047
00048 //...Sag correction...
00049 /* if (doSagCorrection) {
00050 HepVector3D dir = w.direction();
00051 HepPoint3D xw(wp[0], wp[1], wp[2]);
00052 HepPoint3D wireBackwardPosition(wb[0], wb[1], wb[2]);
00053 w.wirePosition(link.positionOnTrack().z(),
00054 xw,
00055 wireBackwardPosition,
00056 dir);
00057 v[0] = dir.x();
00058 v[1] = dir.y();
00059 v[2] = dir.z();
00060 wp[0] = xw.x();
00061 wp[1] = xw.y();
00062 wp[2] = xw.z();
00063 wb[0] = wireBackwardPosition.x();
00064 wb[1] = wireBackwardPosition.y();
00065 wb[2] = wireBackwardPosition.z();
00066 }
00067 */
00068 //...Cal. dPhi to rotate...
00069 const HepPoint3D & xc = _helix.center();
00070
00071 //std::cout<<" helix center: "<<xc<<std::endl;
00072
00073 double xt[3]; _helix.x(0., xt);
00074 double x0 = - xc.x();
00075 double y0 = - xc.y();
00076 double x1 = wp[0] + x0;
00077 double y1 = wp[1] + y0;
00078 x0 += xt[0];
00079 y0 += xt[1];
00080 double dPhi = atan2(x0 * y1 - y0 * x1, x0 * x1 + y0 * y1);
00081
00082 //std::cout<<"dPhi is "<<dPhi<<std::endl;
00083
00084 //...Setup...
00085 double kappa = _helix.kappa();
00086 double phi0 = _helix.phi0();
00087
00088 //...Axial case...
00089 /* if (!w.stereo()) {
00090 positionOnTrack = _helix.x(dPhi);
00091 HepPoint3D x(wp[0], wp[1], wp[2]);
00092 x.setZ(positionOnTrack.z());
00093 positionOnWire = x;
00094 //link.dPhi(dPhi);
00095 std::cout<<" in axial wire : positionOnTrack is "<<positionOnTrack
00096 <<" positionOnWire is "<<positionOnWire<<std::endl;
00097 return (positionOnTrack - positionOnWire).mag();
00098 }
00099 */
00100 double firstdfdphi = 0.;
00101 static bool first = true;
00102 if (first) {
00103 // extern BelleTupleManager * BASF_Histogram;
00104 // BelleTupleManager * m = BASF_Histogram;
00105 // h_nTrial = m->histogram("TTrack::approach nTrial", 100, 0., 100.);
00106 }
00107 //#endif
00108
00109 //...Stereo case...
00110 double rho = Helix::ConstantAlpha / kappa;
00111 double tanLambda = _helix.tanl();
00112 static HepPoint3D x;
00113 double t_x[3];
00114 double t_dXdPhi[3];
00115 const double convergence = 1.0e-5;
00116 double l;
00117 unsigned nTrial = 0;
00118 while (nTrial < 100) {
00119
00120 // x = link.positionOnTrack(_helix->x(dPhi));
00121 positionOnTrack = _helix.x(dPhi);
00122 x = _helix.x(dPhi);
00123 t_x[0] = x[0];
00124 t_x[1] = x[1];
00125 t_x[2] = x[2];
00126
00127 l = v[0] * t_x[0] + v[1] * t_x[1] + v[2] * t_x[2]
00128 - v[0] * wb[0] - v[1] * wb[1] - v[2] * wb[2];
00129
00130 double rcosPhi = rho * cos(phi0 + dPhi);
00131 double rsinPhi = rho * sin(phi0 + dPhi);
00132 t_dXdPhi[0] = rsinPhi;
00133 t_dXdPhi[1] = - rcosPhi;
00134 t_dXdPhi[2] = - rho * tanLambda;
00135
00136 //...f = d(Distance) / d phi...
00137 double t_d2Xd2Phi[2];
00138 t_d2Xd2Phi[0] = rcosPhi;
00139 t_d2Xd2Phi[1] = rsinPhi;
00140
00141 //...iw new...
00142 double n[3];
00143 n[0] = t_x[0] - wb[0];
00144 n[1] = t_x[1] - wb[1];
00145 n[2] = t_x[2] - wb[2];
00146
00147 double a[3];
00148 a[0] = n[0] - l * v[0];
00149 a[1] = n[1] - l * v[1];
00150 a[2] = n[2] - l * v[2];
00151 double dfdphi = a[0] * t_dXdPhi[0]
00152 + a[1] * t_dXdPhi[1]
00153 + a[2] * t_dXdPhi[2];
00154
00155 //#ifdef TRKRECO_DEBUG
00156 if (nTrial == 0) {
00157 // break;
00158 firstdfdphi = dfdphi;
00159 }
00160
00161 //...Check bad case...
00162 if (nTrial > 3) {
00163 std::cout<<" bad case, calculate approach ntrial = "<<nTrial<< endl;
00164 }
00165 //#endif
00166
00167 //...Is it converged?...
00168 if (fabs(dfdphi) < convergence)
00169 break;
00170
00171 double dv = v[0] * t_dXdPhi[0]
00172 + v[1] * t_dXdPhi[1]
00173 + v[2] * t_dXdPhi[2];
00174 double t0 = t_dXdPhi[0] * t_dXdPhi[0]
00175 + t_dXdPhi[1] * t_dXdPhi[1]
00176 + t_dXdPhi[2] * t_dXdPhi[2];
00177 double d2fd2phi = t0 - dv * dv
00178 + a[0] * t_d2Xd2Phi[0]
00179 + a[1] * t_d2Xd2Phi[1];
00180 // + a[2] * t_d2Xd2Phi[2];
00181
00182 dPhi -= dfdphi / d2fd2phi;
00183
00184 // cout << "nTrial=" << nTrial << endl;
00185 // cout << "iw f,df,dphi=" << dfdphi << "," << d2fd2phi << "," << dPhi << endl;
00186
00187 ++nTrial;
00188 }
00189 //...Cal. positions...
00190 positionOnWire[0] = wb[0] + l * v[0];
00191 positionOnWire[1] = wb[1] + l * v[1];
00192 positionOnWire[2] = wb[2] + l * v[2];
00193
00194 //std::cout<<" positionOnTrack is "<<positionOnTrack
00195 // <<" positionOnWire is "<<positionOnWire<<std::endl;
00196 return (positionOnTrack - positionOnWire).mag();
00197
00198 //link.dPhi(dPhi);
00199
00200 // #ifdef TRKRECO_DEBUG
00201 // h_nTrial->accumulate((float) nTrial + .5);
00202 // #endif
00203 // return nTrial;
00204 }
|
|
|
|
|
|
|
|
|
01817 { back_ = i;}
|
|
|
01818 { return back_;}
|
|
|
|
|
|
sets/returns z componet of the magnetic field.
|
|
|
00300 {
00301 return m_bField;
00302 }
|
|
|
sets/returns z componet of the magnetic field.
00280 {
00281 m_bField = a;
00282 m_alpha = 10000. / 2.99792458 / m_bField;
00283
00284 //std::cout<<"m_alpha: "<<m_alpha<<std::endl;
00285 updateCache();
00286 return m_bField;
00287 }
|
|
|
returns position of helix center(z = 0.);
|
|
|
returns position of helix center(z = 0.);
00194 {
00195 return m_center;
00196 }
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
||||||||||||||||
|
|
|
||||||||||||
|
02793 {
02794
02795 double lr = HitMdc.LR();
02796 const KalFitWire& Wire = HitMdc.wire();
02797 int wire_ID = Wire.geoID();
02798 int layerid = HitMdc.wire().layer().layerId();
02799 double entrangle = HitMdc.rechitptr()->getEntra();
02800
02801 HepPoint3D fwd(Wire.fwd());
02802 HepPoint3D bck(Wire.bck());
02803 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
02804 Helix work = H;
02805 work.ignoreErrorMatrix();
02806 work.pivot((fwd + bck) * .5);
02807 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
02808 H.pivot(x0kal);
02809
02810 Hep3Vector meas = H.momentum(0).cross(wire).unit();
02811
02812 if (wire_ID<0 || wire_ID>6796){ //bes
02813 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
02814 << std::endl;
02815 return DBL_MAX;
02816 }
02817
02818 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
02819 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
02820 double tofest(0);
02821 double phi = fmod(phi0() + M_PI4, M_PI2);
02822 double csf0 = cos(phi);
02823 double snf0 = (1. - csf0) * (1. + csf0);
02824 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
02825 if(phi > M_PI) snf0 = - snf0;
02826
02827 if (Tof_correc_) {
02828 Hep3Vector ip(0, 0, 0);
02829 Helix work = *(Helix*)this;
02830 work.ignoreErrorMatrix();
02831 work.pivot(ip);
02832 double phi_ip = work.phi0();
02833 if (fabs(phi - phi_ip) > M_PI) {
02834 if (phi > phi_ip) phi -= 2 * M_PI;
02835 else phi_ip -= 2 * M_PI;
02836 }
02837 double t = tanl();
02838 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
02839 double pmag( sqrt( 1.0 + t*t ) / kappa());
02840 double mass_over_p( mass_ / pmag );
02841 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
02842 tofest = l / ( 29.9792458 * beta );
02843 // if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
02844 }
02845
02846 const HepSymMatrix& ea = H.Ea();
02847 const HepVector& v_a = H.a();
02848 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
02849
02850 HepVector v_H(5, 0);
02851 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
02852 v_H[3] = -meas.z();
02853 HepMatrix v_HT = v_H.T();
02854
02855 double estim = (v_HT * v_a)[0];
02856 HepVector ea_v_H = ea * v_H;
02857 HepMatrix ea_v_HT = (ea_v_H).T();
02858 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
02859
02860 HepSymMatrix eaNewL(5), eaNewR(5);
02861 HepVector aNewL(5), aNewR(5);
02862
02863 //double time = HitMdc.tdc();
02864 //if (Tof_correc_)
02865 // time = time - tofest;
02866 double drifttime = getDriftTime(HitMdc , tofest);
02867 double ddl = getDriftDist(HitMdc, drifttime , 0 );
02868 double ddr = getDriftDist(HitMdc, drifttime , 1 );
02869 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
02870 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
02871
02872 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
02873 double er_dmeas2[2] = {0. , 0.};
02874 if(resolflag_ == 1) {
02875 er_dmeas2[0] = 0.1*erddl;
02876 er_dmeas2[1] = 0.1*erddr;
02877 }
02878 else if(resolflag_ == 0) {
02879 // int layid = HitMdc.wire().layer().layerId();
02880 // double sigma = getSigma(layid, dd);
02881 // er_dmeas2[0] = er_dmeas2[1] = sigma;
02882 }
02883
02884 if ((LR_==0 && lr != 1.0) ||
02885 (LR_==1 && lr == -1.0)) {
02886
02887 double er_dmeasL, dmeasL;
02888 if(Tof_correc_) {
02889 dmeasL = (-1.0)*fabs(dmeas2[0]);
02890 er_dmeasL = er_dmeas2[0];
02891 } else {
02892 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
02893 er_dmeasL = HitMdc.erdist()[0];
02894 }
02895
02896 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
02897 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
02898 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
02899 if(0. == RkL) RkL = 1.e-4;
02900
02901 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
02902 aNewL = v_a + diffL;
02903 double sigL = dmeasL -(v_HT * aNewL)[0];
02904 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
02905 }
02906
02907 if ((LR_==0 && lr != -1.0) ||
02908 (LR_==1 && lr == 1.0)) {
02909
02910 double er_dmeasR, dmeasR;
02911 if(Tof_correc_) {
02912 dmeasR = dmeas2[1];
02913 er_dmeasR = er_dmeas2[1];
02914 } else {
02915 dmeasR = fabs(HitMdc.dist()[1]);
02916 er_dmeasR = HitMdc.erdist()[1];
02917 }
02918
02919 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
02920 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
02921 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
02922 if(0. == RkR) RkR = 1.e-4;
02923
02924 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
02925 aNewR = v_a + diffR;
02926 double sigR = dmeasR -(v_HT * aNewR)[0];
02927 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
02928 }
02929
02930 if (dchi2R < dchi2L){
02931 H.a(aNewR); H.Ea(eaNewR);
02932 } else {
02933 H.a(aNewL); H.Ea(eaNewL);
02934 }
02935 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
02936 }
|
|
||||||||||||||||
|
02938 {
02939
02940 double lr = HitMdc.LR();
02941 const KalFitWire& Wire = HitMdc.wire();
02942 int wire_ID = Wire.geoID();
02943 int layerid = HitMdc.wire().layer().layerId();
02944 double entrangle = HitMdc.rechitptr()->getEntra();
02945
02946 HepPoint3D fwd(Wire.fwd());
02947 HepPoint3D bck(Wire.bck());
02948 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
02949 Helix work = H;
02950 work.ignoreErrorMatrix();
02951 work.pivot((fwd + bck) * .5);
02952 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
02953 H.pivot(x0kal);
02954
02955 Hep3Vector meas = H.momentum(0).cross(wire).unit();
02956
02957 if (wire_ID<0 || wire_ID>6796){ //bes
02958 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
02959 << std::endl;
02960 return DBL_MAX;
02961 }
02962
02963 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
02964 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
02965 double tofest(0);
02966 double phi = fmod(phi0() + M_PI4, M_PI2);
02967 double csf0 = cos(phi);
02968 double snf0 = (1. - csf0) * (1. + csf0);
02969 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
02970 if(phi > M_PI) snf0 = - snf0;
02971
02972 if (Tof_correc_) {
02973 Hep3Vector ip(0, 0, 0);
02974 Helix work = *(Helix*)this;
02975 work.ignoreErrorMatrix();
02976 work.pivot(ip);
02977 double phi_ip = work.phi0();
02978 if (fabs(phi - phi_ip) > M_PI) {
02979 if (phi > phi_ip) phi -= 2 * M_PI;
02980 else phi_ip -= 2 * M_PI;
02981 }
02982 double t = tanl();
02983 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
02984 double pmag( sqrt( 1.0 + t*t ) / kappa());
02985 double mass_over_p( mass_ / pmag );
02986 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
02987 tofest = l / ( 29.9792458 * beta );
02988 if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
02989 }
02990
02991 const HepSymMatrix& ea = H.Ea();
02992 const HepVector& v_a = H.a();
02993 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
02994
02995 HepVector v_H(5, 0);
02996 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
02997 v_H[3] = -meas.z();
02998 HepMatrix v_HT = v_H.T();
02999
03000 double estim = (v_HT * v_a)[0];
03001 HepVector ea_v_H = ea * v_H;
03002 HepMatrix ea_v_HT = (ea_v_H).T();
03003 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
03004
03005 HepSymMatrix eaNewL(5), eaNewR(5);
03006 HepVector aNewL(5), aNewR(5);
03007
03008 //double time = HitMdc.tdc();
03009 //if (Tof_correc_)
03010 // time = time - tofest;
03011 double drifttime = getDriftTime(HitMdc , tofest);
03012 double ddl = getDriftDist(HitMdc, drifttime , 0 );
03013 double ddr = getDriftDist(HitMdc, drifttime , 1 );
03014 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
03015 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
03016
03017 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
03018 double er_dmeas2[2] = {0. , 0.};
03019 if(resolflag_ == 1) {
03020 er_dmeas2[0] = 0.1*erddl;
03021 er_dmeas2[1] = 0.1*erddr;
03022 }
03023 else if(resolflag_ == 0) {
03024 // int layid = HitMdc.wire().layer().layerId();
03025 // double sigma = getSigma(layid, dd);
03026 // er_dmeas2[0] = er_dmeas2[1] = sigma;
03027 }
03028
03029 if ((LR_==0 && lr != 1.0) ||
03030 (LR_==1 && lr == -1.0)) {
03031
03032 double er_dmeasL, dmeasL;
03033 if(Tof_correc_) {
03034 dmeasL = (-1.0)*fabs(dmeas2[0]);
03035 er_dmeasL = er_dmeas2[0];
03036 } else {
03037 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
03038 er_dmeasL = HitMdc.erdist()[0];
03039 }
03040
03041 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
03042 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
03043 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
03044 if(0. == RkL) RkL = 1.e-4;
03045
03046 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
03047 aNewL = v_a + diffL;
03048 double sigL = dmeasL -(v_HT * aNewL)[0];
03049 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
03050 }
03051
03052 if ((LR_==0 && lr != -1.0) ||
03053 (LR_==1 && lr == 1.0)) {
03054
03055 double er_dmeasR, dmeasR;
03056 if(Tof_correc_) {
03057 dmeasR = dmeas2[1];
03058 er_dmeasR = er_dmeas2[1];
03059 } else {
03060 dmeasR = fabs(HitMdc.dist()[1]);
03061 er_dmeasR = HitMdc.erdist()[1];
03062 }
03063
03064 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
03065 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
03066 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
03067 if(0. == RkR) RkR = 1.e-4;
03068
03069 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
03070 aNewR = v_a + diffR;
03071 double sigR = dmeasR -(v_HT * aNewR)[0];
03072 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
03073 }
03074
03075 if (dchi2R < dchi2L){
03076 H.a(aNewR); H.Ea(eaNewR);
03077 } else {
03078 H.a(aNewL); H.Ea(eaNewL);
03079 }
03080 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
03081 }
|
|
|
00214 { chiSq_ = c; }
|
|
|
00183 { return chiSq_; }
|
|
|
00214 { chiSq_ = c; }
|
|
|
00183 { return chiSq_; }
|
|
|
00215 { chiSq_back_ = c; }
|
|
|
00184 { return chiSq_back_; }
|
|
|
00215 { chiSq_back_ = c; }
|
|
|
00184 { return chiSq_back_; }
|
|
||||||||||||
|
Correct the error according the current tanl value :.
|
|
||||||||||||
|
Correct the error according the current tanl value :.
|
|
|
|
|
|
00312 {
00313 return m_cp;
00314 }
|
|
|
|
|
|
00248 {
00249 return m_r;
00250 }
|
|
|
00196 { return dchi2_max_; }
|
|
|
00196 { return dchi2_max_; }
|
|
||||||||||||
|
|
|
||||||||||||
|
00605 {
00606 //
00607 // Calculate Jacobian (@4m/@a)
00608 // Vector a is helix parameters and phi is internal parameter.
00609 // Vector 4m is 4 momentum.
00610 //
00611
00612 HepMatrix d4MDA(4,5,0);
00613
00614 double phi0 = m_ac[1];
00615 double cpa = m_ac[2];
00616 double tnl = m_ac[4];
00617
00618 double cosf0phi = cos(phi0+phi);
00619 double sinf0phi = sin(phi0+phi);
00620
00621 double rho;
00622 if(cpa != 0.)rho = 1./cpa;
00623 else rho = (DBL_MAX);
00624
00625 double charge = 1.;
00626 if(cpa < 0.)charge = -1.;
00627
00628 double E = sqrt(rho*rho*(1.+tnl*tnl)+mass*mass);
00629
00630 d4MDA[0][1] = -fabs(rho)*cosf0phi;
00631 d4MDA[0][2] = charge*rho*rho*sinf0phi;
00632
00633 d4MDA[1][1] = -fabs(rho)*sinf0phi;
00634 d4MDA[1][2] = -charge*rho*rho*cosf0phi;
00635
00636 d4MDA[2][2] = -charge*rho*rho*tnl;
00637 d4MDA[2][4] = fabs(rho);
00638
00639 if (cpa != 0.0 && E != 0.0) {
00640 d4MDA[3][2] = (-1.-tnl*tnl)/(cpa*cpa*cpa*E);
00641 d4MDA[3][4] = tnl/(cpa*cpa*E);
00642 } else {
00643 d4MDA[3][2] = (DBL_MAX);
00644 d4MDA[3][4] = (DBL_MAX);
00645 }
00646 return d4MDA;
00647 }
|
|
||||||||||||
|
|
|
||||||||||||
|
00651 {
00652 //
00653 // Calculate Jacobian (@4mx/@a)
00654 // Vector a is helix parameters and phi is internal parameter.
00655 // Vector 4xm is 4 momentum and position.
00656 //
00657
00658 HepMatrix d4MXDA(7,5,0);
00659
00660 const double & dr = m_ac[0];
00661 const double & phi0 = m_ac[1];
00662 const double & cpa = m_ac[2];
00663 const double & dz = m_ac[3];
00664 const double & tnl = m_ac[4];
00665
00666 double cosf0phi = cos(phi0+phi);
00667 double sinf0phi = sin(phi0+phi);
00668
00669 double rho;
00670 if(cpa != 0.)rho = 1./cpa;
00671 else rho = (DBL_MAX);
00672
00673 double charge = 1.;
00674 if(cpa < 0.)charge = -1.;
00675
00676 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
00677
00678 d4MXDA[0][1] = - fabs(rho) * cosf0phi;
00679 d4MXDA[0][2] = charge * rho * rho * sinf0phi;
00680
00681 d4MXDA[1][1] = - fabs(rho) * sinf0phi;
00682 d4MXDA[1][2] = - charge * rho * rho * cosf0phi;
00683
00684 d4MXDA[2][2] = - charge * rho * rho * tnl;
00685 d4MXDA[2][4] = fabs(rho);
00686
00687 if (cpa != 0.0 && E != 0.0) {
00688 d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
00689 d4MXDA[3][4] = tnl / (cpa * cpa * E);
00690 } else {
00691 d4MXDA[3][2] = (DBL_MAX);
00692 d4MXDA[3][4] = (DBL_MAX);
00693 }
00694
00695 d4MXDA[4][0] = m_cp;
00696 d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
00697 if (cpa != 0.0) {
00698 d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
00699 } else {
00700 d4MXDA[4][2] = (DBL_MAX);
00701 }
00702
00703 d4MXDA[5][0] = m_sp;
00704 d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
00705 if (cpa != 0.0) {
00706 d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);
00707
00708 d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
00709 } else {
00710 d4MXDA[5][2] = (DBL_MAX);
00711
00712 d4MXDA[6][2] = (DBL_MAX);
00713 }
00714
00715 d4MXDA[6][3] = 1.;
00716 d4MXDA[6][4] = - m_r * phi;
00717
00718 return d4MXDA;
00719 }
|
|
|
|
|
|
00446 {
00447 //
00448 // Calculate Jacobian (@ap/@a)
00449 // Vector ap is new helix parameters and a is old helix parameters.
00450 //
00451
00452 HepMatrix dApDA(5,5,0);
00453
00454 const double & dr = m_ac[0];
00455 const double & phi0 = m_ac[1];
00456 const double & cpa = m_ac[2];
00457 const double & dz = m_ac[3];
00458 const double & tnl = m_ac[4];
00459
00460 double drp = ap[0];
00461 double phi0p = ap[1];
00462 double cpap = ap[2];
00463 double dzp = ap[3];
00464 double tnlp = ap[4];
00465
00466 double rdr = m_r + dr;
00467 double rdrpr;
00468 if ((m_r + drp) != 0.0) {
00469 rdrpr = 1. / (m_r + drp);
00470 } else {
00471 rdrpr = (DBL_MAX);
00472 }
00473 // double csfd = cos(phi0)*cos(phi0p) + sin(phi0)*sin(phi0p);
00474 // double snfd = cos(phi0)*sin(phi0p) - sin(phi0)*cos(phi0p);
00475 double csfd = cos(phi0p - phi0);
00476 double snfd = sin(phi0p - phi0);
00477 double phid = fmod(phi0p - phi0 + M_PI8, M_PI2);
00478 if (phid > M_PI) phid = phid - M_PI2;
00479
00480 dApDA[0][0] = csfd;
00481 dApDA[0][1] = rdr*snfd;
00482 if(cpa!=0.0) {
00483 dApDA[0][2] = (m_r/cpa)*( 1.0 - csfd );
00484 } else {
00485 dApDA[0][2] = (DBL_MAX);
00486 }
00487
00488 dApDA[1][0] = - rdrpr*snfd;
00489 dApDA[1][1] = rdr*rdrpr*csfd;
00490 if(cpa!=0.0) {
00491 dApDA[1][2] = (m_r/cpa)*rdrpr*snfd;
00492 } else {
00493 dApDA[1][2] = (DBL_MAX);
00494 }
00495
00496 dApDA[2][2] = 1.0;
00497
00498 dApDA[3][0] = m_r*rdrpr*tnl*snfd;
00499 dApDA[3][1] = m_r*tnl*(1.0 - rdr*rdrpr*csfd);
00500 if(cpa!=0.0) {
00501 dApDA[3][2] = (m_r/cpa)*tnl*(phid - m_r*rdrpr*snfd);
00502 } else {
00503 dApDA[3][2] = (DBL_MAX);
00504 }
00505 dApDA[3][3] = 1.0;
00506 dApDA[3][4] = - m_r*phid;
00507
00508 dApDA[4][4] = 1.0;
00509
00510 return dApDA;
00511 }
|
|
|
|
|
|
00568 {
00569 //
00570 // Calculate Jacobian (@m/@a)
00571 // Vector a is helix parameters and phi is internal parameter.
00572 // Vector m is momentum.
00573 //
00574
00575 HepMatrix dMDA(3,5,0);
00576
00577 const double & phi0 = m_ac[1];
00578 const double & cpa = m_ac[2];
00579 const double & tnl = m_ac[4];
00580
00581 double cosf0phi = cos(phi0+phi);
00582 double sinf0phi = sin(phi0+phi);
00583
00584 double rho;
00585 if(cpa != 0.)rho = 1./cpa;
00586 else rho = (DBL_MAX);
00587
00588 double charge = 1.;
00589 if(cpa < 0.)charge = -1.;
00590
00591 dMDA[0][1] = -fabs(rho)*cosf0phi;
00592 dMDA[0][2] = charge*rho*rho*sinf0phi;
00593
00594 dMDA[1][1] = -fabs(rho)*sinf0phi;
00595 dMDA[1][2] = -charge*rho*rho*cosf0phi;
00596
00597 dMDA[2][2] = -charge*rho*rho*tnl;
00598 dMDA[2][4] = fabs(rho);
00599
00600 return dMDA;
00601 }
|
|
|
|
|
|
00514 {
00515 //
00516 // Calculate Jacobian (@x/@a)
00517 // Vector a is helix parameters and phi is internal parameter
00518 // which specifys the point to be calculated for Ex(phi).
00519 //
00520
00521 HepMatrix dXDA(3,5,0);
00522
00523 const double & dr = m_ac[0];
00524 const double & phi0 = m_ac[1];
00525 const double & cpa = m_ac[2];
00526 const double & dz = m_ac[3];
00527 const double & tnl = m_ac[4];
00528
00529 double cosf0phi = cos(phi0 + phi);
00530 double sinf0phi = sin(phi0 + phi);
00531
00532 dXDA[0][0] = m_cp;
00533 dXDA[0][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
00534 if(cpa!=0.0) {
00535 dXDA[0][2] = - (m_r / cpa) * (m_cp - cosf0phi);
00536 } else {
00537 dXDA[0][2] = (DBL_MAX);
00538 }
00539 // dXDA[0][3] = 0.0;
00540 // dXDA[0][4] = 0.0;
00541
00542 dXDA[1][0] = m_sp;
00543 dXDA[1][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
00544 if(cpa!=0.0) {
00545 dXDA[1][2] = - (m_r / cpa) * (m_sp - sinf0phi);
00546 } else {
00547 dXDA[1][2] = (DBL_MAX);
00548 }
00549 // dXDA[1][3] = 0.0;
00550 // dXDA[1][4] = 0.0;
00551
00552 // dXDA[2][0] = 0.0;
00553 // dXDA[2][1] = 0.0;
00554 if(cpa!=0.0) {
00555 dXDA[2][2] = (m_r / cpa) * tnl * phi;
00556 } else {
00557 dXDA[2][2] = (DBL_MAX);
00558 }
00559 dXDA[2][3] = 1.0;
00560 dXDA[2][4] = - m_r * phi;
00561
00562 return dXDA;
00563 }
|
|
|
returns direction vector after rotating angle dPhi in phi direction.
|
|
|
returns direction vector after rotating angle dPhi in phi direction.
00212 {
00213 return momentum(phi).unit();
00214 }
|
|
|
returns an element of parameters.
|
|
|
returns an element of parameters.
00218 {
00219 return m_ac[0];
00220 }
|
|
|
|
|
|
00236 {
00237 return m_ac[3];
00238 }
|
|
|
sets helix paramters and error matrix.
|
|
|
returns error matrix.
|
|
|
sets helix paramters and error matrix.
00274 {
00275 return m_Ea = i;
00276 }
|
|
|
returns error matrix.
00260 {
00261 return m_Ea;
00262 }
|
|
|
00152 { return Ea_forMdc_; }
|
|
|
00152 { return Ea_forMdc_; }
|
|
|
00148 { return Ea_last_; }
|
|
|
00148 { return Ea_last_; }
|
|
||||||||||||||||
|
Calculate total energy lost in material.
|
|
||||||||||||||||
|
Calculate total energy lost in material.
00277 {
00278 #ifdef YDEBUG
00279 cout<<"eloss():ea before: "<<Ea()<<endl;
00280 #endif
00281 HepVector v_a = a();
00282 double v_a_2_2 = v_a[2] * v_a[2];
00283 double v_a_4_2 = v_a[4] * v_a[4];
00284 double pmag = 1 / fabs(v_a[2]) * sqrt(1 + v_a_4_2);
00285 double psq = pmag * pmag;
00286 double E = sqrt(mass_ * mass_ + psq);
00287 double dE = material.dE(mass_, path, pmag);
00288 //std::cout<<" dE1: "<<dE<<std::endl;
00289
00290 if (index > 0)
00291 psq += dE * (dE + 2 * sqrt(mass_ * mass_ + psq));
00292 else
00293 psq += dE * (dE - 2 * sqrt(mass_ * mass_ + psq));
00294
00295 if (tofall_ && index < 0){
00296 // Kaon case :
00297 if (p_kaon_ > 0){
00298 double psq_kaon = p_kaon_ * p_kaon_;
00299 double dE_kaon = material.dE(MASS[3], path, p_kaon_);
00300 psq_kaon += dE_kaon * (dE_kaon -
00301 2 * sqrt(MASS[3] * MASS[3] + psq_kaon));
00302 if (psq_kaon < 0) psq_kaon = 0;
00303 p_kaon_ = sqrt(psq_kaon);
00304 }
00305 // Proton case :
00306 if (p_proton_ > 0){
00307 double psq_proton = p_proton_ * p_proton_;
00308 double dE_proton = material.dE(MASS[4], path, p_proton_);
00309 psq_proton += dE_proton * (dE_proton -
00310 2 * sqrt(MASS[4] * MASS[4] + psq_proton));
00311 if (psq_proton < 0) psq_proton = 0;
00312 p_proton_ = sqrt(psq_proton);
00313 }
00314 }
00315
00316 double dpt;
00317 if(psq < 0) dpt = 0;
00318 else dpt = v_a[2] * pmag / sqrt(psq);
00319 #ifdef YDEBUG
00320 cout<<"eloss():k: "<<v_a[2]<<" k' "<<dpt<<endl;
00321 #endif
00322 // attempt to take account of energy loss for error matrix
00323
00324 HepSymMatrix ea = Ea();
00325 double r_E_prim = (E + dE)/E;
00326
00327 // 1/ Straggling in the energy loss :
00328 if (strag_){
00329 double del_E(0);
00330 if(l_mass_==0) {
00331 del_E = dE*factor_strag_;
00332 } else {
00333 del_E = material.del_E(mass_, path, pmag);
00334 }
00335
00336 ea[2][2] += (v_a[2] * v_a[2]) * E * E * del_E * del_E / (psq*psq);
00337 }
00338
00339 // Effect of the change of curvature (just variables change):
00340 double dpt2 = dpt*dpt;
00341 double A = dpt*dpt2/(v_a_2_2*v_a[2])*r_E_prim;
00342 double B = v_a[4]/(1+v_a_4_2)*
00343 dpt*(1-dpt2/v_a_2_2*r_E_prim);
00344
00345 double ea_2_0 = A*ea[2][0] + B*ea[4][0];
00346 double ea_2_1 = A*ea[2][1] + B*ea[4][1];
00347 double ea_2_2 = A*A*ea[2][2] + 2*A*B*ea[2][4] + B*B*ea[4][4];
00348 double ea_2_3 = A*ea[2][3] + B*ea[4][3];
00349 double ea_2_4 = A*ea[2][4] + B*ea[4][4];
00350
00351 v_a[2] = dpt;
00352 a(v_a);
00353
00354 ea[2][0] = ea_2_0;
00355 //std::cout<<"ea[2][0] is "<<ea[2][0]<<" ea(3,1) is "<<ea(3,1)<<std::endl;
00356 ea[2][1] = ea_2_1;
00357 //std::cout<<"ea[2][2] is "<<ea[2][2]<<" ea(3,3) is "<<ea(3,3)<<std::endl;
00358 ea[2][2] = ea_2_2;
00359 ea[2][3] = ea_2_3;
00360 ea[2][4] = ea_2_4;
00361
00362 Ea(ea);
00363 //cout<<"eloss():dE: "<<dE<<endl;
00364 //cout<<"eloss():A: "<<A<<" B: "<<B<<endl;
00365 //cout<<"eloss():ea after: "<<Ea()<<endl;
00366 r0_ = fabs(center().perp() - fabs(radius()));
00367 }
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||||||
|
01244 {
01245 HepMatrix m_HT = m_H.T();
01246 HepPoint3D x0 = x(0);
01247 HepVector v_x(3);
01248 v_x[0] = x0.x();
01249 v_x[1] = x0.y();
01250 v_x[2] = x0.z();
01251 HepMatrix m_X = delXDelA(0);
01252 HepMatrix m_XT = m_X.T();
01253 HepMatrix m_C = m_X * Ea() * m_XT;
01254 //int err;
01255 HepVector m_K = 1 / (m_V + (m_HT * m_C * m_H)[0]) * m_H;
01256 HepVector v_a = a();
01257 HepSymMatrix ea = Ea();
01258 HepMatrix mXTmK = m_XT * m_K;
01259 double v_r = v_m - v_d - (m_HT * v_x)[0];
01260 v_a += ea * mXTmK * v_r;
01261 a(v_a);
01262 ea.assign(ea - ea * mXTmK * m_HT * m_X * ea);
01263 Ea(ea);
01264 // Record last hit included parameters :
01265 update_last();
01266 HepMatrix mCmK = m_C * m_K;
01267 v_x += mCmK * v_r;
01268 m_C -= mCmK * m_HT * m_C;
01269 v_r = v_m - v_d - (m_HT * v_x)[0];
01270 double m_R = m_V - (m_HT * m_C * m_H)[0];
01271 double chiSq = v_r * v_r / m_R;
01272 chiSq_ += chiSq;
01273 return chiSq;
01274 }
|
|
|
|
|
|
01294 {
01295 fiTerm_ = fi;
01296 }
|
|
|
|
|
|
|
|
||||||||||||||||
|
|
|
||||||||||||||||
|
00196 {
00197 int layerid = hitmdc.wire().layer().layerId();
00198 int cellid = MdcID::wire(hitmdc.rechitptr()->getMdcId());
00199 if(debug_ == 4 ){
00200 std::cout<<"the cellid is .."<<cellid<<std::endl;
00201 }
00202 double entrangle = hitmdc.rechitptr()->getEntra();
00203
00204 //std::cout<<" entrangle: "<<entrangle<<std::endl;
00205
00206 return CalibFunSvc_->driftTimeToDist(drifttime, layerid, cellid, lr, entrangle);
00207 }
|
|
||||||||||||
|
|
|
||||||||||||
|
00076 {
00077 const double vinner = 220.0e8; // cm/s
00078 const double vouter = 240.0e8; // cm/s
00079
00080 int layerid = hitmdc.wire().layer().layerId();
00081 double zhit = (hitmdc.rechitptr())->getZhit();
00082 const KalFitWire* wire = &(hitmdc.wire());
00083 HepPoint3D fPoint = wire->fwd();
00084 HepPoint3D bPoint = wire->bck();
00085
00086 // unit is centimeter
00087 double wireLen = (fPoint-bPoint).x()*(fPoint-bPoint).x()+(fPoint-bPoint).y()*(fPoint-bPoint).y()+(fPoint-bPoint).z()*(fPoint-bPoint).z();
00088 wireLen = sqrt(wireLen);
00089 double wireZLen = fabs(fPoint.z() - bPoint.z());
00090 double tp = 0.;
00091 double vp = 0.;
00092 // west readout
00093 if(0==layerid%2){
00094 // inner chamber
00095 if(layerid<8){
00096 vp = vinner;
00097 }
00098 else {
00099 vp = vouter;
00100 }
00101 tp = wireLen*fabs(zhit - bPoint.z())/wireZLen/vp;
00102 }
00103
00104 // east readout
00105 if(1==layerid%2){
00106 // inner chamber
00107 if(layerid<8){
00108 vp = vinner;
00109 }
00110 else {
00111 vp = vouter;
00112 }
00113 tp = wireLen*fabs(zhit - fPoint.z())/wireZLen/vp;
00114 }
00115
00116 // s to ns
00117 tp = 1.0e9*tp;
00118
00119 if(0==tprop_) tp = 0.;
00120
00121 //std::cout<<"propogation time: "<<tp<<std::endl;
00122
00123 int wireid = hitmdc.wire().geoID();
00124 double drifttime1(0.);
00125 double drifttime2(0.);
00126 double drifttime3(0.);
00127
00128 MdcDigiCol::iterator iter = mdcDigiCol_->begin();
00129 for(; iter != mdcDigiCol_->end(); iter++ ) {
00130 if((*iter)->identify() == (hitmdc.rechitptr())->getMdcId()) break;
00131 }
00132 //double t0 = get_T0();
00133 // see the code of wulh in /Mdc/MdcCalibFunSvc/MdcCalibFunSvc/MdcCalibFunSvc.h,
00134 // double getT0(int wireid) const { return m_t0[wireid]; }
00135
00136 //double getTimeWalk(int layid, double Q) const ;
00137 double Q = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
00138 double timewalk = CalibFunSvc_->getTimeWalk(layerid, Q);
00139
00140 if(debug_ == 4) {
00141 std::cout<<"CalibFunSvc_->getTimeWalk, timewalk ="<<timewalk<<std::endl;
00142 }
00143
00144 double timeoffset = CalibFunSvc_->getT0(wireid);
00145 if(debug_ == 4) {
00146 std::cout<<"CalibFunSvc_->getT0, timeoffset ="<<timeoffset<<std::endl;
00147 }
00148
00149 double eventt0 = getT0();
00150
00151 if(debug_ == 4) {
00152 std::cout<<"the Event T0 we get in the function getDriftTime(...) is "<<eventt0<<std::endl;
00153 }
00154
00155 // this tdc value come from MdcRecHit assigned by zhangyao
00156 double tdctime1 = hitmdc.tdc();
00157 double tdctime2(0.);
00158 double tdctime3(0.);
00159
00160 if(debug_ == 4) {
00161 std::cout<<"tdctime1 be here is .."<<tdctime1<<std::endl;
00162 }
00163
00164 // this tdc value come from MdcDigiCol time channel
00165 // 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 !
00166 if(debug_ == 4) {
00167 // std::cout<<"the size of the mdcDigiCol_ is "<<mdcDigiCol_.size()<<std::endl;
00168 }
00169 // MdcDigiCol::iterator iter = mdcDigiCol_->begin();
00170 // for(; iter != mdcDigiCol_->end(); iter++ ) {
00171 // if((*iter)->identify() == (hitmdc.rechitptr())->getMdcId()) break;
00172 // }
00173 if(debug_ == 4) {
00174 std::cout<<"the time channel be here is .."<<(*iter)->getTimeChannel()<<std::endl;
00175 }
00176 tdctime2 = RawDataUtil::MdcTime((*iter)->getTimeChannel());
00177 tdctime3 = hitmdc.rechitptr()->getTdc();
00178 drifttime1 = tdctime1 - eventt0 - toftime - timewalk -timeoffset - tp;
00179 drifttime2 = tdctime2 - eventt0 - toftime - timewalk -timeoffset - tp;
00180 drifttime3 = tdctime3 - eventt0 - toftime - timewalk -timeoffset - tp;
00181 if(debug_ == 4 ) {
00182 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;
00183 }
00184 //return drifttime3;
00185 //return drifttime1;
00186 if(drifttime_choice_ == 0)
00187 return drifttime2;
00188 if(drifttime_choice_ == 1)
00189 // use the driftT caluculated by track-finding
00190 return hitmdc.rechitptr()->getDriftT();
00191 }
|
|
|
00164 { return fiTerm_;}
|
|
|
00164 { return fiTerm_;}
|
|
|
|
|
|
00044 {
00045 return initMatrix_ ;
00046 }
|
|
|
00158 { return pathSM_;}
|
|
|
00158 { return pathSM_;}
|
|
||||||||||||||||||||
|
|
|
||||||||||||
|
|
|
||||||||||||||||||||
|
00211 {
00212 int layerid = hitmdc.wire().layer().layerId();
00213 double entrangle = hitmdc.rechitptr()->getEntra();
00214 // double tanlam = hitmdc.rechitptr()->getTanl();
00215 double z = hitmdc.rechitptr()->getZhit();
00216 double Q = hitmdc.rechitptr()->getAdc();
00217 //std::cout<<" the driftdist before getsigma is "<<dist<<" the layer is"<<layerid<<std::endl;
00218 double temp = CalibFunSvc_->getSigma(layerid, lr, dist, entrangle, tanlam, z , Q );
00219 //std::cout<<" the sigma is "<<temp<<std::endl;
00220 return temp;
00221 }
|
|
||||||||||||
|
03084 {
03085 double sigma1,sigma2,f;
03086 driftDist *= 10;//mm
03087 if(layerId<8){
03088 if(driftDist<0.5){
03089 sigma1=0.112784; sigma2=0.229274; f=0.666;
03090 }else if(driftDist<1.){
03091 sigma1=0.103123; sigma2=0.269797; f=0.934;
03092 }else if(driftDist<1.5){
03093 sigma1=0.08276; sigma2=0.17493; f=0.89;
03094 }else if(driftDist<2.){
03095 sigma1=0.070109; sigma2=0.149859; f=0.89;
03096 }else if(driftDist<2.5){
03097 sigma1=0.064453; sigma2=0.130149; f=0.886;
03098 }else if(driftDist<3.){
03099 sigma1=0.062383; sigma2=0.138806; f=0.942;
03100 }else if(driftDist<3.5){
03101 sigma1=0.061873; sigma2=0.145696; f=0.946;
03102 }else if(driftDist<4.){
03103 sigma1=0.061236; sigma2=0.119584; f=0.891;
03104 }else if(driftDist<4.5){
03105 sigma1=0.066292; sigma2=0.148426; f=0.917;
03106 }else if(driftDist<5.){
03107 sigma1=0.078074; sigma2=0.188148; f=0.911;
03108 }else if(driftDist<5.5){
03109 sigma1=0.088657; sigma2=0.27548; f=0.838;
03110 }else{
03111 sigma1=0.093089; sigma2=0.115556; f=0.367;
03112 }
03113 }else{
03114 if(driftDist<0.5){
03115 sigma1=0.112433; sigma2=0.327548; f=0.645;
03116 }else if(driftDist<1.){
03117 sigma1=0.096703; sigma2=0.305206; f=0.897;
03118 }else if(driftDist<1.5){
03119 sigma1=0.082518; sigma2=0.248913; f= 0.934;
03120 }else if(driftDist<2.){
03121 sigma1=0.072501; sigma2=0.153868; f= 0.899;
03122 }else if(driftDist<2.5){
03123 sigma1= 0.065535; sigma2=0.14246; f=0.914;
03124 }else if(driftDist<3.){
03125 sigma1=0.060497; sigma2=0.126489; f=0.918;
03126 }else if(driftDist<3.5){
03127 sigma1=0.057643; sigma2= 0.112927; f=0.892;
03128 }else if(driftDist<4.){
03129 sigma1=0.055266; sigma2=0.094833; f=0.887;
03130 }else if(driftDist<4.5){
03131 sigma1=0.056263; sigma2=0.124419; f= 0.932;
03132 }else if(driftDist<5.){
03133 sigma1=0.056599; sigma2=0.124248; f=0.923;
03134 }else if(driftDist<5.5){
03135 sigma1= 0.061377; sigma2=0.146147; f=0.964;
03136 }else if(driftDist<6.){
03137 sigma1=0.063978; sigma2=0.150591; f=0.942;
03138 }else if(driftDist<6.5){
03139 sigma1=0.072951; sigma2=0.15685; f=0.913;
03140 }else if(driftDist<7.){
03141 sigma1=0.085438; sigma2=0.255109; f=0.931;
03142 }else if(driftDist<7.5){
03143 sigma1=0.101635; sigma2=0.315529; f=0.878;
03144 }else{
03145 sigma1=0.149529; sigma2=0.374697; f=0.89;
03146 }
03147 }
03148 double sigmax = sqrt(f*sigma1*sigma1+(1 - f)*sigma2*sigma2)*0.1;
03149 return sigmax;//cm
03150 }
|
|
|
|
|
|
00059 {
00060 //------------------get event start time-----------
00061
00062 if(debug_ == 4) {
00063 std::cout<<"in function KalFitTrack::getT0 ( ), EventT0_ = "<<EventT0_<<std::endl;
00064 }
00065 return EventT0_ ;
00066 }
|
|
|
00161 { return tof2_;}
|
|
|
00161 { return tof2_;}
|
|
|
00235 { return HelixSegs_[i];}
|
|
|
00235 { return HelixSegs_[i];}
|
|
|
00234 { return HelixSegs_;}
|
|
|
|
|
|
00234 { return HelixSegs_;}
|
|
|
|
|
|
00230 { return HitsMdc_[i];}
|
|
|
00230 { return HitsMdc_[i];}
|
|
|
00229 { return HitsMdc_;}
|
|
|
|
|
|
00229 { return HitsMdc_;}
|
|
|
01824 { HitsMdc_ = lh;}
|
|
|
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.
|
|
|
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.
00722 {
00723 m_matrixValid = false;
00724 m_Ea *= 0.;
00725 }
|
|
|
00209 { insist_ = t;}
|
|
|
Extractor :.
00178 { return insist_; }
|
|
|
00209 { insist_ = t;}
|
|
|
Extractor :.
00178 { return insist_; }
|
|
|
Intersection with different geometry.
|
|
|
Intersection with different geometry.
00123 {
00124 double m_rad = radius();
00125 double l = center().perp();
00126
00127 double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
00128
00129 if(cosPhi < -1 || cosPhi > 1) return 0;
00130
00131 double dPhi = center().phi() - acos(cosPhi) - phi0();
00132
00133 if(dPhi < -M_PI) dPhi += 2 * M_PI;
00134
00135 return dPhi;
00136 }
|
|
|
|
|
|
00175 {
00176 if (tanl() != 0 && radius() != 0)
00177 return (pivot().z() + dz() - z) / (radius() * tanl());
00178 else return 0;
00179 }
|
|
||||||||||||
|
|
|
||||||||||||
|
00158 {
00159 HepPoint3D xc = plane * center();
00160 double r = radius();
00161 double d = r * r - (x - xc.x()) * (x - xc.x());
00162 if(d < 0) return 0;
00163
00164 double yy = xc.y();
00165 if(yy > 0) yy -= sqrt(d);
00166 else yy += sqrt(d);
00167
00168 double l = (plane.inverse() *
00169 HepPoint3D(x, yy, 0)).perp();
00170
00171 return intersect_cylinder(l);
00172 }
|
|
||||||||||||
|
|
|
||||||||||||
|
00140 {
00141 HepPoint3D xc = plane * center();
00142 double r = radius();
00143 double d = r * r - (y - xc.y()) * (y - xc.y());
00144 if(d < 0) return 0;
00145
00146 double xx = xc.x();
00147 if(xx > 0) xx -= sqrt(d);
00148 else xx += sqrt(d);
00149
00150 double l = (plane.inverse() *
00151 HepPoint3D(xx, y, 0)).perp();
00152
00153 return intersect_cylinder(l);
00154 }
|
|
|
|
|
|
00230 {
00231 return m_ac[2];
00232 }
|
|
|
|
|
|
Magnetic field map.
|
|
|
01815 { lead_ = i;}
|
|
|
Magnetic field map.
01816 { return lead_;}
|
|
|
|
|
|
01823 { LR_ = x;}
|
|
|
|
|
|
00182 { return mass_; }
|
|
|
01809 { return MASS[i];}
|
|
|
00182 { return mass_; }
|
|
|
00190 { return mom_; }
|
|
|
00190 { return mom_; }
|
|
||||||||||||||||||||
|
returns 4momentum vector after rotating angle dPhi in phi direction.
|
|
||||||||||||||||
|
returns 4momentum vector after rotating angle dPhi in phi direction.
|
|
||||||||||||
|
returns 4momentum vector after rotating angle dPhi in phi direction.
|
|
||||||||||||
|
returns momentum vector after rotating angle dPhi in phi direction.
|
|
|
returns momentum vector after rotating angle dPhi in phi direction.
|
|
||||||||||||||||||||
|
returns 4momentum vector after rotating angle dPhi in phi direction.
00280 {
00281 //
00282 // Calculate momentum.
00283 //
00284 // Pt = | 1/kappa | (GeV/c)
00285 //
00286 // Px = -Pt * sin(phi0 + phi)
00287 // Py = Pt * cos(phi0 + phi)
00288 // Pz = Pt * tan(lambda)
00289 //
00290 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
00291
00292 double pt = fabs(m_pt);
00293 double px = - pt * sin(m_ac[1] + phi);
00294 double py = pt * cos(m_ac[1] + phi);
00295 double pz = pt * m_ac[4];
00296 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
00297
00298 x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
00299 x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
00300 x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);
00301
00302 if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi,mass));
00303 else Emx = m_Ea;
00304
00305 return HepLorentzVector(px, py, pz, E);
00306 }
|
|
||||||||||||||||
|
returns 4momentum vector after rotating angle dPhi in phi direction.
00252 {
00253 //
00254 // Calculate momentum.
00255 //
00256 // Pt = | 1/kappa | (GeV/c)
00257 //
00258 // Px = -Pt * sin(phi0 + phi)
00259 // Py = Pt * cos(phi0 + phi)
00260 // Pz = Pt * tan(lambda)
00261 //
00262 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
00263
00264 double pt = fabs(m_pt);
00265 double px = - pt * sin(m_ac[1] + phi);
00266 double py = pt * cos(m_ac[1] + phi);
00267 double pz = pt * m_ac[4];
00268 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
00269
00270 if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi,mass));
00271 else Em = m_Ea;
00272
00273 return HepLorentzVector(px, py, pz, E);
00274 }
|
|
||||||||||||
|
returns 4momentum vector after rotating angle dPhi in phi direction.
00229 {
00230 //
00231 // Calculate momentum.
00232 //
00233 // Pt = | 1/kappa | (GeV/c)
00234 //
00235 // Px = -Pt * sin(phi0 + phi)
00236 // Py = Pt * cos(phi0 + phi)
00237 // Pz = Pt * tan(lambda)
00238 //
00239 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
00240
00241 double pt = fabs(m_pt);
00242 double px = - pt * sin(m_ac[1] + phi);
00243 double py = pt * cos(m_ac[1] + phi);
00244 double pz = pt * m_ac[4];
00245 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
00246
00247 return HepLorentzVector(px, py, pz, E);
00248 }
|
|
||||||||||||
|
returns momentum vector after rotating angle dPhi in phi direction.
00206 {
00207 //
00208 // Calculate momentum.
00209 //
00210 // Pt = | 1/kappa | (GeV/c)
00211 //
00212 // Px = -Pt * sin(phi0 + phi)
00213 // Py = Pt * cos(phi0 + phi)
00214 // Pz = Pt * tan(lambda)
00215 //
00216
00217 double pt = fabs(m_pt);
00218 double px = - pt * sin(m_ac[1] + phi);
00219 double py = pt * cos(m_ac[1] + phi);
00220 double pz = pt * m_ac[4];
00221
00222 if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
00223 else Em = m_Ea;
00224
00225 return Hep3Vector(px, py, pz);
00226 }
|
|
|
returns momentum vector after rotating angle dPhi in phi direction.
00186 {
00187 //
00188 // Calculate momentum.
00189 //
00190 // Pt = | 1/kappa | (GeV/c)
00191 //
00192 // Px = -Pt * sin(phi0 + phi)
00193 // Py = Pt * cos(phi0 + phi)
00194 // Pz = Pt * tan(lambda)
00195 //
00196
00197 double pt = fabs(m_pt);
00198 double px = - pt * sin(m_ac[1] + phi);
00199 double py = pt * cos(m_ac[1] + phi);
00200 double pz = pt * m_ac[4];
00201
00202 return Hep3Vector(px, py, pz);
00203 }
|
|
||||||||||||||||
|
|
|
||||||||||||||||
|
00183 {
00184 HepSymMatrix ea = Ea();
00185 //cout<<"ms():path "<<path<<endl;
00186 //cout<<"ms():ea before: "<<ea<<endl;
00187 double k = kappa();
00188 double t = tanl();
00189 double t2 = t * t;
00190 double pt2 = 1 + t2;
00191 double k2 = k * k;
00192
00193 double pmag = 1 / fabs(k) * sqrt(pt2);
00194 double dth = material.mcs_angle(mass_, path, pmag);
00195 double dth2 = dth * dth;
00196 double pt2dth2 = pt2 * dth2;
00197
00198 ea[1][1] += pt2dth2;
00199 ea[2][2] += k2 * t2 * dth2;
00200 ea[2][4] += k * t * pt2dth2;
00201 ea[4][4] += pt2 * pt2dth2;
00202
00203 ea[3][3] += dth2 * path * path /3 / (1 + t2);
00204 ea[3][4] += dth2 * path/2 * sqrt(1 + t2);
00205 ea[3][2] += dth2 * t / sqrt(1 + t2) * k * path/2;
00206 ea[0][0] += dth2 * path * path/3;
00207 ea[0][1] += dth2 * sqrt(1 + t2) * path/2;
00208
00209 Ea(ea);
00210 //cout<<"ms():ms angle in this: "<<dth<<endl;
00211 //cout<<"ms():ea after: "<<Ea()<<endl;
00212 if (index < 0) {
00213 double x0 = material.X0();
00214 if (x0) path_rd_ += path/x0;
00215 }
00216 }
|
|
||||||||||||
|
Calculate multiple scattering angle.
|
|
||||||||||||
|
Calculate multiple scattering angle.
00219 {
00220 double k = kappa();
00221 double t = tanl();
00222 double t2 = t * t;
00223 double k2 = k * k;
00224
00225 double pmag = ( 1 / fabs(k) ) * sqrt(1 + t2);
00226 double psq = pmag*pmag;
00227 /*
00228 double Zprims = 3/2*0.076 + 0.580/9.79*4.99*(4.99+1) +
00229 0.041/183.85*74*(74+1) + 0.302/26.98 * 13 * (13+1);
00230 double chicc = 0.00039612 * sqrt(Zprims * 0.001168);
00231 double dth = 2.557 * chicc * sqrt(path * (mass_*mass_ + psq)) / psq;
00232 */
00233
00234 //std::cout<<" mdcGasRadlen: "<<mdcGasRadlen_<<std::endl;
00235 double pathx = path/mdcGasRadlen_;
00236 double dth = 0.0136* sqrt(pathx * (mass_*mass_ + psq))/psq
00237 *(1 + 0.038 * log(pathx));;
00238 HepSymMatrix ea = Ea();
00239 #ifdef YDEBUG
00240 cout<<"msgasmdc():path "<<path<<" pathx "<<pathx<<endl;
00241 cout<<"msgasmdc():dth "<<dth<<endl;
00242 cout<<"msgasmdc():ea before: "<<ea<<endl;
00243 #endif
00244 double dth2 = dth * dth;
00245
00246 ea[1][1] += (1 + t2) * dth2;
00247 ea[2][2] += k2 * t2 * dth2;
00248 ea[2][4] += k * t * (1 + t2) * dth2;
00249 ea[4][4] += (1 + t2) * (1 + t2) * dth2;
00250
00251 // additionnal terms (terms proportional to l and l^2)
00252
00253 ea[3][3] += dth2 * path * path /3 / (1 + t2);
00254 ea[3][4] += dth2 * path/2 * sqrt(1 + t2);
00255 ea[3][2] += dth2 * t / sqrt(1 + t2) * k * path/2;
00256 ea[0][0] += dth2 * path * path/3;
00257 ea[0][1] += dth2 * sqrt(1 + t2) * path/2;
00258
00259 Ea(ea);
00260 #ifdef YDEBUG
00261 cout<<"msgasmdc():ea after: "<<Ea()<<endl;
00262 #endif
00263 if (index < 0) {
00264 pathip_ += path;
00265 // RMK : put by hand, take care !!
00266 double x0 = mdcGasRadlen_; // for the Mdc gas
00267 path_rd_ += path/x0;
00268 tof(path);
00269 #ifdef YDEBUG
00270 cout<<"ms...pathip..."<<pathip_<<endl;
00271 #endif
00272 }
00273 }
|
|
|
00200 { return ncath_; }
|
|
|
00200 { return ncath_; }
|
|
|
00217 { nchits_ = n; }
|
|
|
00198 { return nchits_; }
|
|
|
00217 { nchits_ = n; }
|
|
|
00198 { return nchits_; }
|
|
|
00216 { ndf_back_ = n; }
|
|
|
00185 { return ndf_back_; }
|
|
|
00216 { ndf_back_ = n; }
|
|
|
00185 { return ndf_back_; }
|
|
|
00203 { return nhit_r_; }
|
|
|
00203 { return nhit_r_; }
|
|
|
00204 { return nhit_z_; }
|
|
|
00204 { return nhit_z_; }
|
|
|
|
|
|
01808 { return NMASS;}
|
|
|
00218 { nster_ = n; }
|
|
|
00199 { return nster_; }
|
|
|
00218 { nster_ = n; }
|
|
|
00199 { return nster_; }
|
|
|
|
|
|
00447 {
00448 unsigned int nhit = HitsMdc_.size();
00449 int Num[50] = {0};
00450 for( unsigned i=0 ; i < nhit; i++ )
00451 Num[HitsMdc_[i].wire().layer().layerId()]++;
00452 for( unsigned i=0 ; i < nhit; i++ )
00453 if (Num[HitsMdc_[i].wire().layer().layerId()]>2)
00454 HitsMdc_[i].chi2(-2);
00455 }
|
|
|
|
|
|
|
|
|
01821 { numf_ = i;}
|
|
|
01822 { return numf_;}
|
|
|
|
|
|
00431 {
00432 for(int it=0; it<HitsMdc().size()-1; it++){
00433 if(HitsMdc_[it].wire().layer().layerId() == HitsMdc_[it+1].wire().layer().layerId())
00434 {
00435 if((kappa()<0)&&(HitsMdc_[it].wire().localId() > HitsMdc_[it+1].wire().localId())){
00436 std::swap(HitsMdc_[it], HitsMdc_[it+1]);
00437 }
00438 if((kappa()>0)&&(HitsMdc_[it].wire().localId() < HitsMdc_[it+1].wire().localId())){
00439 std::swap(HitsMdc_[it], HitsMdc_[it+1]);
00440 }
00441 }
00442 }
00443 }
|
|
|
Modifier Order the wire hits for mdc track |
|
|
Modifier Order the wire hits for mdc track 00370 {
00371 unsigned int nhit = HitsMdc_.size();
00372 Helix tracktest = *(Helix*)this;
00373 int ind = 0;
00374 double* Rad = new double[nhit];
00375 double* Ypos = new double[nhit];
00376 for( unsigned i=0 ; i < nhit; i++ ){
00377
00378 HepPoint3D fwd(HitsMdc_[i].wire().fwd());
00379 HepPoint3D bck(HitsMdc_[i].wire().bck());
00380 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
00381
00382 // Modification for stereo wires :
00383 Helix work = tracktest;
00384 work.ignoreErrorMatrix();
00385 work.pivot((fwd + bck) * .5);
00386 HepPoint3D x0 = (work.x(0).z() - bck.z())
00387 / wire.z() * wire + bck;
00388
00389 tracktest.pivot(x0);
00390 Rad[ind] = tracktest.x(0).perp();
00391 Ypos[ind] = x0.y();
00392 ind++;
00393 //cout<<"Ypos: "<<Ypos[ind-1]<<endl;
00394 }
00395
00396 // Reorder...
00397 if (index < 0)
00398 for(int j, k = nhit - 1; k >= 0; k = j){
00399 j = -1;
00400 for(int i = 1; i <= k; i++)
00401 if(Rad[i - 1] > Rad[i]){
00402 j = i - 1;
00403 std::swap(Rad[i], Rad[j]);
00404 std::swap(HitsMdc_[i], HitsMdc_[j]);
00405 }
00406 }
00407 if (index > 0)
00408 for(int j, k = nhit - 1; k >= 0; k = j){
00409 j = -1;
00410 for(int i = 1; i <= k; i++)
00411 if(Rad[i - 1] < Rad[i]){
00412 j = i - 1;
00413 std::swap(Rad[i], Rad[j]);
00414 std::swap(HitsMdc_[i], HitsMdc_[j]);
00415 }
00416 }
00417 if (index == 0)
00418 for(int j, k = nhit - 1; k >= 0; k = j){
00419 j = -1;
00420 for(int i = 1; i <= k; i++)
00421 if(Ypos[i - 1] > Ypos[i]){
00422 j = i - 1;
00423 std::swap(Ypos[i], Ypos[j]);
00424 std::swap(HitsMdc_[i], HitsMdc_[j]);
00425 }
00426 }
00427 delete [] Rad;
00428 delete [] Ypos;
00429 }
|
|
|
00212 { p_kaon_ = pl;}
|
|
|
00194 { return p_kaon_; }
|
|
|
00212 { p_kaon_ = pl;}
|
|
|
00194 { return p_kaon_; }
|
|
|
00213 { p_proton_ = pl;}
|
|
|
00195 { return p_proton_; }
|
|
|
00213 { p_proton_ = pl;}
|
|
|
00195 { return p_proton_; }
|
|
|
00201 { return pat1_; }
|
|
|
00201 { return pat1_; }
|
|
|
00202 { return pat2_; }
|
|
|
00202 { return pat2_; }
|
|
|
00188 { return path_ab_; }
|
|
|
00188 { return path_ab_; }
|
|
|
Update the path length estimation.
|
|
|
Update the path length estimation.
|
|
|
00187 { return path_rd_; }
|
|
|
00187 { return path_rd_; }
|
|
|
00211 { pathip_ = pl;}
|
|
|
00186 { return pathip_; }
|
|
|
00211 { pathip_ = pl;}
|
|
|
00186 { return pathip_; }
|
|
|
Function to calculate the path length in the layer.
|
|
|
00189 { return PathL_; }
|
|
|
Function to calculate the path length in the layer.
|
|
|
00189 { return PathL_; }
|
|
|
|
|
|
00224 {
00225 return m_ac[1];
00226 }
|
|
|
sets pivot position.
|
|
|
returns pivot position.
|
|
|
sets pivot position.
00310 {
00311
00312 //std::cout<<" in Helix::pivot:"<<std::endl;
00313 //std::cout<<" m_alpha: "<<m_alpha<<std::endl;
00314
00315 const double & dr = m_ac[0];
00316 const double & phi0 = m_ac[1];
00317 const double & kappa = m_ac[2];
00318 const double & dz = m_ac[3];
00319 const double & tanl = m_ac[4];
00320
00321 double rdr = dr + m_r;
00322 double phi = fmod(phi0 + M_PI4, M_PI2);
00323 double csf0 = cos(phi);
00324 double snf0 = (1. - csf0) * (1. + csf0);
00325 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
00326 if(phi > M_PI) snf0 = - snf0;
00327
00328 double xc = m_pivot.x() + rdr * csf0;
00329 double yc = m_pivot.y() + rdr * snf0;
00330 double csf, snf;
00331 if(m_r != 0.0) {
00332 csf = (xc - newPivot.x()) / m_r;
00333 snf = (yc - newPivot.y()) / m_r;
00334 double anrm = sqrt(csf * csf + snf * snf);
00335 if(anrm != 0.0) {
00336 csf /= anrm;
00337 snf /= anrm;
00338 phi = atan2(snf, csf);
00339 } else {
00340 csf = 1.0;
00341 snf = 0.0;
00342 phi = 0.0;
00343 }
00344 } else {
00345 csf = 1.0;
00346 snf = 0.0;
00347 phi = 0.0;
00348 }
00349 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
00350 if(phid > M_PI) phid = phid - M_PI2;
00351 double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
00352 * csf
00353 + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
00354 double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();
00355
00356 HepVector ap(5);
00357 ap[0] = drp;
00358 ap[1] = fmod(phi + M_PI4, M_PI2);
00359 ap[2] = kappa;
00360 ap[3] = dzp;
00361 ap[4] = tanl;
00362
00363 // if (m_matrixValid) m_Ea.assign(delApDelA(ap) * m_Ea * delApDelA(ap).T());
00364 if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));
00365
00366 m_a = ap;
00367 m_pivot = newPivot;
00368
00369 //...Are these needed?...iw...
00370 updateCache();
00371 return m_pivot;
00372 }
|
|
|
returns pivot position.
00200 {
00201 return m_pivot;
00202 }
|
|
|
00150 { return pivot_forMdc_;}
|
|
|
00150 { return pivot_forMdc_;}
|
|
|
returns helix parameters
00146 { return pivot_last_; }
|
|
|
returns helix parameters
00146 { return pivot_last_; }
|
|
||||||||||||
|
|
|
|
Sets pivot position in a given mag field.
|
|
||||||||||||
|
01466 {
01467
01468 Helix tracktest = *(Helix*)this;
01469 tracktest.ignoreErrorMatrix();
01470 double tl = a()[4];
01471 double th = 90.0 - 180.0*M_1_PI*atan( tl );
01472 /*
01473 int nstep(1);
01474 if (steplev_ == 1)
01475 nstep = 3;
01476 else if (steplev_ == 2 && (th > 140 || th <45))
01477 if ((pivot()-newPivot).mag()<3.)
01478 nstep = 3;
01479 else
01480 nstep = 6;
01481 */
01482 Hep3Vector delta_x((newPivot-pivot()).x()/double(outer_steps_),
01483 (newPivot-pivot()).y()/double(outer_steps_),
01484 (newPivot-pivot()).z()/double(outer_steps_));
01485 int i = 1;
01486
01487 while (i <= outer_steps_) {
01488 HepPoint3D nextPivot(pivot()+delta_x);
01489 double xnp(nextPivot.x()), ynp(nextPivot.y()), znp(nextPivot.z());
01490
01491 HepSymMatrix Ea_now = Ea();
01492 HepPoint3D piv(pivot());
01493 double xp(piv.x()), yp(piv.y()), zp(piv.z());
01494
01495 double dr = a()[0];
01496 double phi0 = a()[1];
01497 double kappa = a()[2];
01498 double dz = a()[3];
01499 double tanl = a()[4];
01500
01501 double m_rad(0);
01502 m_rad = radius();
01503
01504 double rdr = dr + m_rad;
01505 double phi = fmod(phi0 + M_PI4, M_PI2);
01506 double csf0 = cos(phi);
01507 double snf0 = (1. - csf0) * (1. + csf0);
01508 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
01509 if(phi > M_PI) snf0 = - snf0;
01510
01511 double xc = xp + rdr * csf0;
01512 double yc = yp + rdr * snf0;
01513 double csf = (xc - xnp) / m_rad;
01514 double snf = (yc - ynp) / m_rad;
01515 double anrm = sqrt(csf * csf + snf * snf);
01516 csf /= anrm;
01517 snf /= anrm;
01518 phi = atan2(snf, csf);
01519 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
01520 if(phid > M_PI) phid = phid - M_PI2;
01521 double drp = (xp + dr * csf0 + m_rad * (csf0 - csf) - xnp)
01522 * csf
01523 + (yp + dr * snf0 + m_rad * (snf0 - snf) - ynp) * snf;
01524 double dzp = zp + dz - m_rad * tanl * phid - znp;
01525
01526 HepVector ap(5);
01527 ap[0] = drp;
01528 ap[1] = fmod(phi + M_PI4, M_PI2);
01529 ap[2] = kappa;
01530 ap[3] = dzp;
01531 ap[4] = tanl;
01532
01533 //std::cout<<" numf_: "<<numf_<<std::endl;
01534
01535 // Modification due to non uniform magnetic field :
01536 if (numf_ > 10) {
01537
01538 Hep3Vector x1(xp + dr*csf0, yp + dr*snf0, zp + dz);
01539
01540 double csf0p = cos(ap[1]);
01541 double snf0p = (1. - csf0p) * (1. + csf0p);
01542 snf0p = sqrt((snf0p > 0.) ? snf0p : 0.);
01543 if(ap[1] > M_PI) snf0p = - snf0p;
01544
01545 Hep3Vector x2(xnp + drp*csf0p,
01546 ynp + drp*snf0p,
01547 znp + dzp);
01548
01549 Hep3Vector dist((x1 - x2).x()/100.0,
01550 (x1 - x2).y()/100.0,
01551 (x1 - x2).z()/100.0);
01552
01553 HepPoint3D middlebis((x1.x() + x2.x())/2,
01554 (x1.y() + x2.y())/2,
01555 (x1.z() + x2.z())/2);
01556
01557 HepVector3D field;
01558 MFSvc_->fieldVector(10.*middlebis,field);
01559 field = 1000.*field;
01560
01561 //std::cout<<"B field: "<<field<<std::endl;
01562 Hep3Vector dB(field.x(),
01563 field.y(),
01564 (field.z()-0.1*Bznom_));
01565
01566
01567 //std::cout<<" dB: "<<dB<<std::endl;
01568
01569
01570 if (field.z()) {
01571 double akappa(fabs(kappa));
01572 double sign = kappa/akappa;
01573 HepVector dp(3);
01574 dp = 0.299792458 * sign * dB.cross(dist);
01575
01576 //std::cout<<"dp: "<<dp<<std::endl;
01577
01578 HepVector dhp(3);
01579 dhp[0] = -akappa*(dp[0]*csf0p+dp[1]*snf0p);
01580 dhp[1] = kappa*akappa*(dp[0]*snf0p-dp[1]*csf0p);
01581 dhp[2] = dp[2]*akappa+dhp[1]*tanl/kappa;
01582
01583
01584 //std::cout<<"dhp: "<<dhp<<std::endl;
01585
01586
01587 ap[1] += dhp[0];
01588 ap[2] += dhp[1];
01589 ap[4] += dhp[2];
01590 }
01591 }
01592 HepMatrix m_del = delApDelA(ap);
01593 Ea_now.assign(m_del * Ea_now * m_del.T());
01594 pivot(nextPivot);
01595 a(ap);
01596 Ea(Ea_now);
01597 i++;
01598
01599 //std::cout<<" a: "<<a()<<std::endl;
01600 }
01601
01602 // Estimation of the path length :
01603 double tanl_0(tracktest.a()[4]);
01604 double phi0_0(tracktest.a()[1]);
01605 double radius_0(tracktest.radius());
01606 tracktest.pivot(newPivot);
01607
01608 double phi0_1 = tracktest.a()[1];
01609 if (fabs(phi0_1 - phi0_0) > M_PI) {
01610 if (phi0_1 > phi0_0) phi0_1 -= 2 * M_PI;
01611 else phi0_0 -= 2 * M_PI;
01612 }
01613 if(phi0_1 == phi0_0) phi0_1 = phi0_0 + 1.e-10;
01614 pathl = fabs(radius_0 * (phi0_1 - phi0_0)
01615 * sqrt(1 + tanl_0 * tanl_0));
01616 return newPivot;
01617 }
|
|
|
Sets pivot position in a given mag field.
01362 {
01363
01364 int nstep(1);
01365 HepPoint3D delta_x((newPivot-pivot()).x()/double(inner_steps_),
01366 (newPivot-pivot()).y()/double(inner_steps_),
01367 (newPivot-pivot()).z()/double(inner_steps_));
01368 int i = 1;
01369
01370 while (i <= inner_steps_) {
01371 HepPoint3D nextPivot(pivot()+delta_x);
01372 double xnp(nextPivot.x()), ynp(nextPivot.y()), znp(nextPivot.z());
01373 HepSymMatrix Ea_now = Ea();
01374 HepPoint3D piv(pivot());
01375 double xp(piv.x()), yp(piv.y()), zp(piv.z());
01376 double dr = a()[0];
01377 double phi0 = a()[1];
01378 double kappa = a()[2];
01379 double dz = a()[3];
01380 double tanl = a()[4];
01381 double m_rad(0);
01382 if (numfcor_ == 1)
01383 m_rad = radius_numf();
01384 else
01385 m_rad = radius();
01386 double rdr = dr + m_rad;
01387 double phi = fmod(phi0 + M_PI4, M_PI2);
01388 double csf0 = cos(phi);
01389 double snf0 = (1. - csf0) * (1. + csf0);
01390 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
01391 if(phi > M_PI) snf0 = - snf0;
01392
01393 double xc = xp + rdr * csf0;
01394 double yc = yp + rdr * snf0;
01395 double csf = (xc - xnp) / m_rad;
01396 double snf = (yc - ynp) / m_rad;
01397 double anrm = sqrt(csf * csf + snf * snf);
01398 csf /= anrm;
01399 snf /= anrm;
01400 phi = atan2(snf, csf);
01401 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
01402 if(phid > M_PI) phid = phid - M_PI2;
01403 double drp = (xp + dr * csf0 + m_rad * (csf0 - csf) - xnp)
01404 * csf
01405 + (yp + dr * snf0 + m_rad * (snf0 - snf) - ynp) * snf;
01406 double dzp = zp + dz - m_rad * tanl * phid - znp;
01407
01408 HepVector ap(5);
01409 ap[0] = drp;
01410 ap[1] = fmod(phi + M_PI4, M_PI2);
01411 ap[2] = kappa;
01412 ap[3] = dzp;
01413 ap[4] = tanl;
01414
01415 // Modification due to non uniform magnetic field :
01416 if (numf_ > 10) {
01417
01418 Hep3Vector x1(xp + dr*csf0, yp + dr*snf0, zp + dz);
01419 double csf0p = cos(ap[1]);
01420 double snf0p = (1. - csf0p) * (1. + csf0p);
01421 snf0p = sqrt((snf0p > 0.) ? snf0p : 0.);
01422 if(ap[1] > M_PI) snf0p = - snf0p;
01423
01424 Hep3Vector x2(xnp + drp*csf0p,
01425 ynp + drp*snf0p,
01426 znp + dzp);
01427 Hep3Vector dist((x1 - x2).x()/100.0,
01428 (x1 - x2).y()/100.0,
01429 (x1 - x2).z()/100.0);
01430 HepPoint3D middlebis((x1.x() + x2.x())/2,
01431 (x1.y() + x2.y())/2,
01432 (x1.z() + x2.z())/2);
01433 HepVector3D field;
01434 MFSvc_->fieldVector(10.*middlebis,field);
01435 field = 1000.*field;
01436 Hep3Vector dB(field.x(),
01437 field.y(),
01438 (field.z()-0.1*Bznom_));
01439 if (field.z()) {
01440 double akappa(fabs(kappa));
01441 double sign = kappa/akappa;
01442 HepVector dp(3);
01443 dp = 0.299792458 * sign * dB.cross(dist);
01444 HepVector dhp(3);
01445 dhp[0] = -akappa*(dp[0]*csf0p+dp[1]*snf0p);
01446 dhp[1] = kappa*akappa*(dp[0]*snf0p-dp[1]*csf0p);
01447 dhp[2] = dp[2]*akappa+dhp[1]*tanl/kappa;
01448 if (numfcor_ == 0){
01449 ap[1] += dhp[0];
01450 }
01451 ap[2] += dhp[1];
01452 ap[4] += dhp[2];
01453 }
01454 }
01455 HepMatrix m_del = delApDelA(ap);
01456 Ea_now.assign(m_del * Ea_now * m_del.T());
01457 pivot(nextPivot);
01458 a(ap);
01459 Ea(Ea_now);
01460 i++;
01461 }
01462 return newPivot;
01463 }
|
|
|
00142 { return point_last_;}
|
|
|
set and give out the last point of the track
00141 { point_last_ = point;}
|
|
|
00142 { return point_last_;}
|
|
|
set and give out the last point of the track
00141 { point_last_ = point;}
|
|
|
00181 { return r0_; }
|
|
|
00181 { return r0_; }
|
|
|
00197 { return r_max_; }
|
|
|
00197 { return r_max_; }
|
|
|
returns radious of helix.
|
|
|
returns radious of helix.
00206 {
00207 return m_r;
00208 }
|
|
|
Estimation of the radius in a given mag field.
|
|
|
Estimation of the radius in a given mag field.
01324 {
01325
01326 double Bz(Bznom_);
01327 //std::cout<<"Bz: "<<Bz<<std::endl;
01328 if (numf_ > 10){
01329 double dr = a()[0];
01330 double phi0 = a()[1];
01331 double dz = a()[3];
01332 double phi = fmod(phi0 + M_PI4, M_PI2);
01333 double csf0 = cos(phi);
01334 double snf0 = (1. - csf0) * (1. + csf0);
01335 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
01336 if(phi > M_PI) snf0 = - snf0;
01337 //XYZPoint
01338 HepPoint3D x0((pivot().x() + dr*csf0),
01339 (pivot().y() + dr*snf0),
01340 (pivot().z() + dz));
01341
01342 //XYZVector b;
01343 HepVector3D b;
01344
01345 //std::cout<<"b: "<<b<<std::endl;
01346
01347 MFSvc_->fieldVector(10.*x0, b);
01348
01349 //std::cout<<"b: "<<b<<std::endl;
01350
01351
01352 b = 10000.*b;
01353 Bz = b.z();
01354 }
01355 if (Bz == 0)
01356 Bz = Bznom_;
01357 double ALPHA_loc = 10000./2.99792458/Bz;
01358 return ALPHA_loc / a()[2];
01359 }
|
|
|
|
|
|
|
|
|
01819 { resolflag_ = i;}
|
|
|
01820 { return resolflag_;}
|
|
||||||||||||||||
|
sets helix pivot position, parameters, and error matrix.
|
|
||||||||||||||||
|
sets helix pivot position, parameters, and error matrix.
00377 {
00378 m_pivot = pivot;
00379 m_a = a;
00380 m_Ea = Ea;
00381 m_matrixValid = true;
00382 updateCache();
00383 }
|
|
|
|
|
|
00238 {
00239 iGeomSvc_ = igeomsvc;
00240 }
|
|
|
|
|
|
00039 {
00040 initMatrix_ = m;
00041 }
|
|
|
|
|
|
00229 {
00230 ISvcLocator* svcLocator = Gaudi::svcLocator();
00231 StatusCode sc = svcLocator->service("MagneticFieldSvc",MFSvc_);
00232 if (sc.isFailure()){
00233 std::cout << "Could not load MagneticFieldSvc!" << std::endl;
00234 }
00235 }
|
|
|
|
|
|
00224 {
00225 CalibFunSvc_ = calibsvc;
00226 }
|
|
|
|
|
|
00243 {
00244 mdcDigiCol_ = digicol;
00245 }
|
|
|
|
|
|
00049 {
00050 //------------------set event start time-----------
00051
00052 EventT0_ = eventstart;
00053 if(debug_ == 4) {
00054 std::cout<<"in function KalFitTrack::setT0(...), EventT0_ = "<<EventT0_<<std::endl;
00055 }
00056 }
|
|
|
|
|
|
00306 {
00307 return m_sp;
00308 }
|
|
||||||||||||||||||||||||
|
|
|
||||||||||||||||||||
|
Kalman smoother for Mdc.
|
|
||||||||||||||||||||||||
|
move the pivot of the helixseg to IP(0,0,0) 00459 {
00460 //
00461 double lr = HitMdc.LR();
00462 // For taking the propagation delay into account :
00463 int wire_ID = HitMdc.wire().geoID(); // from 0 to 6796
00464 int layerid = HitMdc.wire().layer().layerId();
00465 double entrangle = HitMdc.rechitptr()->getEntra();
00466 if (wire_ID<0 || wire_ID>6796){ //bes
00467 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID << std::endl;
00468 return 0;
00469 }
00470
00471 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
00472 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
00473
00474 double tofest(0);
00475 double dd(0.);
00476 double phi = fmod(phi0() + M_PI4, M_PI2);
00477 double csf0 = cos(phi);
00478 double snf0 = (1. - csf0) * (1. + csf0);
00479 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
00480 if(phi > M_PI) snf0 = - snf0;
00481
00482 if (Tof_correc_) {
00483 Hep3Vector ip(0, 0, 0);
00484 Helix work = *(Helix*)this;
00485 work.ignoreErrorMatrix();
00486 work.pivot(ip);
00487 double phi_ip = work.phi0();
00488 if (fabs(phi - phi_ip) > M_PI) {
00489 if (phi > phi_ip) phi -= 2 * M_PI;
00490 else phi_ip -= 2 * M_PI;
00491 }
00492 double t = tanl();
00493 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
00494 double pmag( sqrt( 1.0 + t*t ) / kappa());
00495 double mass_over_p( mass_ / pmag );
00496 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
00497 tofest = l / ( 29.9792458 * beta );
00498 if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
00499 }
00500
00501 const HepSymMatrix& ea = Ea();
00502 const HepVector& v_a = a();
00503
00504 HepPoint3D pivot_work = pivot();
00505
00506 double dchi2R(99999999), dchi2L(99999999);
00507
00508 HepVector v_H(5, 0);
00509 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
00510 v_H[3] = -meas.z();
00511 HepMatrix v_HT = v_H.T();
00512
00513 double estim = (v_HT * v_a)[0];
00514 HepVector ea_v_H = ea * v_H;
00515 HepMatrix ea_v_HT = (ea_v_H).T();
00516 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
00517
00518 HepSymMatrix eaNew(5);
00519 HepVector aNew(5);
00520
00521 //double time = HitMdc.tdc();
00522 // if (Tof_correc_)
00523 // time -= tofest;
00524
00525 double drifttime = getDriftTime(HitMdc , tofest);
00526 //cout<<"drifttime = "<<drifttime<<endl;
00527 seg.dt(drifttime);
00528 double ddl = getDriftDist(HitMdc, drifttime, 0);
00529 double ddr = getDriftDist(HitMdc, drifttime, 1);
00530 double erddl = getSigma( HitMdc, a()[4], 0, ddl);
00531 double erddr = getSigma( HitMdc, a()[4], 1, ddr);
00532
00533 // double dd = 1.0e-4 * 40.0 * time;
00534 double dmeas2[2] = { 0.1*ddl, 0.1*ddr }; //millimeter to centimeter
00535 double er_dmeas2[2] = {0. , 0.};
00536 if(resolflag_== 1) {
00537 er_dmeas2[0] = 0.1*erddl;
00538 er_dmeas2[1] = 0.1*erddr;
00539 }
00540 else if(resolflag_ == 0) {
00541 //int layid = HitMdc.wire().layer().layerId();
00542 //double sigma = getSigma(layid, dd);
00543 //er_dmeas2[0] = er_dmeas2[1] = sigma;
00544 }
00545
00546 // Left hypothesis :
00547 if (lr == -1) {
00548 double er_dmeasL, dmeasL;
00549 if(Tof_correc_) {
00550 dmeasL = (-1.0)*dmeas2[0]; // the dmeas&erdmeas calculated by myself
00551 er_dmeasL = er_dmeas2[0];
00552 } else {
00553 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]); // the dmeas&erdmeas calculated by track-finding algorithm
00554 er_dmeasL = HitMdc.erdist()[0];
00555 }
00556
00557 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
00558 eaNew.assign(ea - ea_v_H * AL * ea_v_HT);
00559 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
00560 if(0. == RkL) RkL = 1.e-4;
00561
00562 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
00563 aNew = v_a + diffL;
00564 double sigL = (dmeasL - (v_HT * aNew)[0]);
00565 dchi2L = (sigL*sigL) / (RkL * er_dmeasL*er_dmeasL);
00566 } else if (lr == 1) {
00567 // Right hypothesis :
00568
00569 double er_dmeasR, dmeasR;
00570 if(Tof_correc_) {
00571 dmeasR = dmeas2[1];
00572 er_dmeasR = er_dmeas2[1];
00573 } else {
00574 dmeasR = fabs(HitMdc.dist()[1]);
00575 er_dmeasR = HitMdc.erdist()[1];
00576 }
00577
00578
00579 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
00580 eaNew.assign(ea - ea_v_H * AR * ea_v_HT);
00581 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
00582 if(0. == RkR) RkR = 1.e-4;
00583
00584 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
00585 aNew = v_a + diffR;
00586 double sigR = (dmeasR- (v_HT * aNew)[0]);
00587 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
00588 }
00589
00590 // Update Kalman result :
00591 #ifdef YDEBUG
00592 cout<<"in smoother_Mdc: lr= "<<lr<<" a: "<<v_a<<" a_NEW: "<<aNew<<endl;
00593 #endif
00594 double dchi2_loc(0);
00595 if ((dchi2R < dchi2cuts_anal[layerid] && lr == 1) ||
00596 (dchi2L < dchi2cuts_anal[layerid] && lr == -1)) {
00597 ndf_back_++;
00598 int chge(1);
00599 if ( !( aNew[2] && fabs(aNew[2]-a()[2]) < 1.0 ) ) chge=0;
00600 if (lr == 1) {
00601 chiSq_back_ += dchi2R;
00602 dchi2_loc = dchi2R;
00603 dd = 0.1*ddr;
00604 } else if (lr == -1) {
00605 chiSq_back_ += dchi2L;
00606 dchi2_loc = dchi2L;
00607 dd = -0.1*ddl;
00608 }
00609 if (chge){
00610 a(aNew);
00611 Ea(eaNew);
00612 }
00613 }
00614 dchi2=dchi2_loc;
00615
00616 seg.doca_include(fabs((v_HT*a())[0]));
00617 seg.doca_exclude(fabs(estim));
00619 Hep3Vector ip(0, 0, 0);
00620 KalFitTrack helixNew1(pivot_work, a(), Ea(), 0, 0, 0);
00621 helixNew1.pivot(ip);
00622 HepVector a_temp1 = helixNew1.a();
00623 HepSymMatrix ea_temp1 = helixNew1.Ea();
00624 seg.Ea(ea_temp1);
00625 seg.a(a_temp1);
00626 seg.Ea_include(ea_temp1);
00627 seg.a_include(a_temp1);
00628
00629 KalFitTrack helixNew2(pivot_work, v_a, ea, 0, 0, 0);
00630 helixNew2.pivot(ip);
00631 HepVector a_temp2 = helixNew2.a();
00632 HepSymMatrix ea_temp2 = helixNew2.Ea();
00633 seg.Ea_exclude(ea_temp2);
00634 seg.a_exclude(a_temp2);
00635
00636 seg.tof(tofest);
00637 seg.dd(dd);
00638
00639 return chiSq_back_;
00640 }
|
|
||||||||||||||||||||
|
Kalman smoother for Mdc. move the pivot of the helixseg to IP(0,0,0) 00942 {
00943
00944
00945 HepPoint3D ip(0., 0., 0.);
00946 // attention! we should not to decide the left&right in the smoother process,
00947 // because we choose the left&right of hits from the filter process.
00948
00949 KalFitHitMdc* HitMdc = seg.HitMdc();
00950 double lr = HitMdc->LR();
00951
00952 // For taking the propagation delay into account :
00953 int layerid = (*HitMdc).wire().layer().layerId();
00954 int wire_ID = HitMdc->wire().geoID();
00955 double entrangle = HitMdc->rechitptr()->getEntra();
00956
00957 if (wire_ID<0 || wire_ID>6796){ //bes
00958 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID << std::endl;
00959 return 0;
00960 }
00961
00962 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
00963 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
00964 double dd(0.);
00965 double tofest(0);
00966 double phi = fmod(phi0() + M_PI4, M_PI2);
00967 double csf0 = cos(phi);
00968 double snf0 = (1. - csf0) * (1. + csf0);
00969 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
00970 if(phi > M_PI) snf0 = - snf0;
00971
00972 if (Tof_correc_) {
00973 Hep3Vector ip(0, 0, 0);
00974 Helix work = *(Helix*)this;
00975 work.ignoreErrorMatrix();
00976 work.pivot(ip);
00977 double phi_ip = work.phi0();
00978 if (fabs(phi - phi_ip) > M_PI) {
00979 if (phi > phi_ip) phi -= 2 * M_PI;
00980 else phi_ip -= 2 * M_PI;
00981 }
00982 double t = tanl();
00983 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
00984 double pmag( sqrt( 1.0 + t*t ) / kappa());
00985 double mass_over_p( mass_ / pmag );
00986 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
00987 tofest = l / ( 29.9792458 * beta );
00988 if(csmflag==1 && (*HitMdc).wire().y()>0.) tofest= -1. * tofest;
00989 }
00990
00991 const HepSymMatrix& ea = Ea();
00992 const HepVector& v_a = a();
00993
00994
00995
00996 HepPoint3D pivot_work = pivot();
00997
00998 /*
00999 HepVector a_work = a();
01000 HepSymMatrix ea_work = Ea();
01001
01002 KalFitTrack helix_work(pivot_work, a_work, ea_work, 0, 0, 0);
01003 helix_work.pivot(ip);
01004
01005 HepVector a_temp = helix_work.a();
01006 HepSymMatrix ea_temp = helix_work.Ea();
01007
01008 seg.Ea_pre_bwd(ea_temp);
01009 seg.a_pre_bwd(a_temp);
01010
01011 */
01012
01013 seg.a_pre_bwd(a());
01014 seg.Ea_pre_bwd(Ea());
01015
01016 double dchi2R(99999999), dchi2L(99999999);
01017 HepVector v_H(5, 0);
01018 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
01019 v_H[3] = -meas.z();
01020 HepMatrix v_HT = v_H.T();
01021
01022 double estim = (v_HT * v_a)[0];
01023 HepVector ea_v_H = ea * v_H;
01024 HepMatrix ea_v_HT = (ea_v_H).T();
01025 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
01026 HepSymMatrix eaNew(5);
01027 HepVector aNew(5);
01028 double time = (*HitMdc).tdc();
01029
01030 //seg.dt(time);
01031 // if (Tof_correc_)
01032 // {
01033 // time -= tofest;
01034 // seg.dt(time);
01035 // }
01036 // double dd = 1.0e-4 * 40.0 * time;
01037 // seg.dd(dd);
01038
01039 double drifttime = getDriftTime(*HitMdc , tofest);
01040 seg.dt(drifttime);
01041 double ddl = getDriftDist(*HitMdc, drifttime, 0 );
01042 double ddr = getDriftDist(*HitMdc, drifttime, 1 );
01043 double erddl = getSigma( *HitMdc, a()[4], 0, ddl);
01044 double erddr = getSigma( *HitMdc, a()[4], 1, ddr);
01045
01046 double dmeas2[2] = { 0.1*ddl, 0.1*ddr }; // millimeter to centimeter
01047 double er_dmeas2[2] = {0., 0.};
01048 if(resolflag_ == 1) {
01049 er_dmeas2[0] = 0.1*erddl;
01050 er_dmeas2[1] = 0.1*erddr;
01051 }else if(resolflag_ == 0)
01052 {
01053 }
01054
01055 // Left hypothesis :
01056 if (lr == -1) {
01057 double er_dmeasL, dmeasL;
01058 if(Tof_correc_) {
01059 dmeasL = (-1.0)*dmeas2[0];
01060 er_dmeasL = er_dmeas2[0];
01061 } else {
01062 dmeasL = (-1.0)*fabs((*HitMdc).dist()[0]);
01063 er_dmeasL = (*HitMdc).erdist()[0];
01064 }
01065
01066
01067 //if(layerid < 4) er_dmeasL*=2.0;
01068
01069 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
01070 eaNew.assign(ea - ea_v_H * AL * ea_v_HT);
01071 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
01072 if(0. == RkL) RkL = 1.e-4;
01073
01074 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
01075 aNew = v_a + diffL;
01076 double sigL = (dmeasL - (v_HT * aNew)[0]);
01077 dchi2L = (sigL*sigL) / (RkL * er_dmeasL*er_dmeasL);
01078 } else if (lr == 1) {
01079 // Right hypothesis :
01080
01081 double er_dmeasR, dmeasR;
01082 if(Tof_correc_) {
01083 dmeasR = dmeas2[1];
01084 er_dmeasR = er_dmeas2[1];
01085 } else {
01086 dmeasR = fabs((*HitMdc).dist()[1]);
01087 er_dmeasR = (*HitMdc).erdist()[1];
01088 }
01089
01090
01091 //if(layerid < 4) er_dmeasR*=2.0;
01092
01093
01094 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
01095 eaNew.assign(ea - ea_v_H * AR * ea_v_HT);
01096 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
01097 if(0. == RkR) RkR = 1.e-4;
01098
01099 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
01100 aNew = v_a + diffR;
01101 double sigR = (dmeasR- (v_HT * aNew)[0]);
01102 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
01103 }
01104
01105 // Update Kalman result :
01106 #ifdef YDEBUG
01107 cout<<"in smoother_Mdc: lr= "<<lr<<" a: "<<v_a<<" a_NEW: "<<aNew<<endl;
01108 #endif
01109 double dchi2_loc(0);
01110 if ((dchi2R < dchi2cuts_calib[layerid] && lr == 1) ||
01111 (dchi2L < dchi2cuts_calib[layerid] && lr == -1)) {
01112
01113 ndf_back_++;
01114 flg = 1;
01115 int chge(1);
01116 if (!(aNew[2] && fabs(aNew[2]-a()[2]) < 1.0)) chge = 0;
01117 if (lr == 1) {
01118 chiSq_back_ += dchi2R;
01119 dchi2_loc = dchi2R;
01120 dd = 0.1*ddr;
01121 // if(debug_ ==4) std::cout<<"in smoother "<<std::endl;
01122
01123 } else if (lr == -1) {
01124 chiSq_back_ += dchi2L;
01125 dchi2_loc = dchi2L;
01126 dd = -0.1*ddl;
01127
01128 }
01129 if (chge){
01130 a(aNew);
01131 Ea(eaNew);
01132 }
01133
01134 seg.a_filt_bwd(aNew);
01135 seg.Ea_filt_bwd(eaNew);
01136
01137 HepVector a_pre_fwd_temp=seg.a_pre_fwd();
01138 if( (a_pre_fwd_temp[1]-seg.a_pre_bwd()[1]) > 3. * M_PI /2.) a_pre_fwd_temp[1]-= M_PI2;
01139 if( (a_pre_fwd_temp[1]-seg.a_pre_bwd()[1]) < -3. * M_PI /2. ) a_pre_fwd_temp[1]+= M_PI2;
01140
01141 seg.a_pre_fwd(a_pre_fwd_temp);
01142
01143 HepVector a_pre_filt_temp=seg.a_filt_fwd();
01144 if( (a_pre_filt_temp[1]-seg.a_pre_bwd()[1]) > 3. * M_PI /2. ) a_pre_filt_temp[1] -= M_PI2;
01145 if( (a_pre_filt_temp[1]-seg.a_pre_bwd()[1]) < -3. * M_PI /2.) a_pre_filt_temp[1] += M_PI2;
01146 seg.a_filt_fwd(a_pre_filt_temp);
01147
01148 /*
01149 KalFitTrack helixNew(pivot_work, aNew, eaNew, 0, 0, 0);
01150 helixNew.pivot(ip);
01151 a_temp = helixNew.a();
01152 ea_temp = helixNew.Ea();
01153 seg.Ea_filt_bwd(ea_temp);
01154 seg.a_filt_bwd(a_temp);
01155 */
01156
01157 int inv;
01158
01159 if(debug_ == 4){
01160 std::cout<<"seg.Ea_pre_bwd().inverse(inv) ..."<<seg.Ea_pre_bwd().inverse(inv)<<std::endl;
01161 std::cout<<"seg.Ea_pre_fwd().inverse(inv) ..."<<seg.Ea_pre_fwd().inverse(inv)<<std::endl;
01162 }
01163
01164 HepSymMatrix Ea_pre_comb = (seg.Ea_pre_bwd().inverse(inv)+seg.Ea_pre_fwd().inverse(inv)).inverse(inv);
01165 seg.Ea_exclude(Ea_pre_comb);
01166
01167
01168 if(debug_ == 4) {
01169 std::cout<<"Ea_pre_comb_ ... "<<Ea_pre_comb<<std::endl;
01170 std::cout<<"seg.a_pre_bwd()..."<<seg.a_pre_bwd()<<std::endl;
01171 std::cout<<"seg.a_pre_fwd()..."<<seg.a_pre_fwd()<<std::endl;
01172 }
01173 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());
01174 seg.a_exclude(a_pre_comb);
01175
01176
01177 if(debug_ == 4) {
01178 std::cout<<"seg.Ea_filt_fwd()..."<<seg.Ea_filt_fwd()<<std::endl;
01179 std::cout<<"seg.Ea_filt_fwd().inverse(inv)..."<<seg.Ea_filt_fwd().inverse(inv)<<std::endl;
01180 }
01181 seg.Ea((seg.Ea_filt_fwd().inverse(inv)+seg.Ea_pre_bwd().inverse(inv)).inverse(inv));
01182 seg.Ea_include((seg.Ea_filt_fwd().inverse(inv)+seg.Ea_pre_bwd().inverse(inv)).inverse(inv));
01183
01184 if(debug_ == 4) {
01185 std::cout<<"seg.Ea() is ..."<<seg.Ea()<<std::endl;
01186 std::cout<<"seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd() ..."<<seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()<<std::endl;
01187 std::cout<<"seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd() ... "<<seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()<<std::endl;
01188 }
01189
01190 seg.a(seg.Ea()*(seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()+seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()));
01191 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()));
01192
01193 if(inv != 0) {
01194 //std::cout<<"ERROR OCCUR WHEN INVERSE MATRIX !"<<std::endl;
01195 }
01196
01197 seg.residual_exclude(dd - (v_HT*a_pre_comb)[0]);
01198 seg.residual_include(dd - (v_HT*seg.a())[0]);
01199 seg.doca_exclude(fabs((v_HT*a_pre_comb)[0]));
01200 seg.doca_include(fabs((v_HT*seg.a())[0]));
01201
01202 if(debug_ == 4){
01203 std::cout<<"the dd in smoother_Mdc is "<<dd
01204 <<" the (v_HT*a_pre_comb)[0] is "<<(v_HT*a_pre_comb)[0]<<std::endl;
01205 }
01206
01207 //compare the two method to calculate the include doca value :
01208 if(debug_ == 4){
01209 std::cout<<"method 1 ..."<<sqrt(seg.a()[0]*seg.a()[0]+seg.a()[3]*seg.a()[3])<<std::endl;
01210 std::cout<<"method 2 ..."<<fabs((v_HT*seg.a())[0])<<std::endl;
01211 }
01212
01213
01215 KalFitTrack helixNew1(pivot_work, seg.a(), seg.Ea(), 0, 0, 0);
01216 helixNew1.pivot(ip);
01217 HepVector a_temp1 = helixNew1.a();
01218 HepSymMatrix ea_temp1 = helixNew1.Ea();
01219 seg.Ea(ea_temp1);
01220 seg.a(a_temp1);
01221 seg.Ea_include(ea_temp1);
01222 seg.a_include(a_temp1);
01223
01224 KalFitTrack helixNew2(pivot_work, seg.a_exclude(), seg.Ea_exclude(), 0, 0, 0);
01225 helixNew2.pivot(ip);
01226 HepVector a_temp2 = helixNew2.a();
01227 HepSymMatrix ea_temp2 = helixNew2.Ea();
01228 seg.Ea_exclude(ea_temp2);
01229 seg.a_exclude(a_temp2);
01230
01231 seg.tof(tofest);
01232 seg.dd(dd);
01233
01234 }
01235 return chiSq_back_;
01236 }
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||||||
|
move the pivot of the helixseg to IP(0,0,0) 00645 {
00646
00647
00648 HepPoint3D ip(0., 0., 0.);
00649 // attention! we should not to decide the left&right in the smoother process,
00650 // because we choose the left&right of hits from the filter process.
00651
00652 KalFitHitMdc* HitMdc = seg.HitMdc();
00653 double lr = HitMdc->LR();
00654
00655 // For taking the propagation delay into account :
00656 int layerid = (*HitMdc).wire().layer().layerId();
00657 int wire_ID = HitMdc->wire().geoID();
00658 double entrangle = HitMdc->rechitptr()->getEntra();
00659
00660 if (wire_ID<0 || wire_ID>6796){ //bes
00661 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID << std::endl;
00662 return 0;
00663 }
00664
00665 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
00666 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
00667 double dd(0.);
00668 double tofest(0);
00669 double phi = fmod(phi0() + M_PI4, M_PI2);
00670 double csf0 = cos(phi);
00671 double snf0 = (1. - csf0) * (1. + csf0);
00672 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
00673 if(phi > M_PI) snf0 = - snf0;
00674
00675 if (Tof_correc_) {
00676 Hep3Vector ip(0, 0, 0);
00677 Helix work = *(Helix*)this;
00678 work.ignoreErrorMatrix();
00679 work.pivot(ip);
00680 double phi_ip = work.phi0();
00681 if (fabs(phi - phi_ip) > M_PI) {
00682 if (phi > phi_ip) phi -= 2 * M_PI;
00683 else phi_ip -= 2 * M_PI;
00684 }
00685 double t = tanl();
00686 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
00687 double pmag( sqrt( 1.0 + t*t ) / kappa());
00688 double mass_over_p( mass_ / pmag );
00689 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
00690 tofest = l / ( 29.9792458 * beta );
00691 if(csmflag==1 && (*HitMdc).wire().y()>0.) tofest= -1. * tofest;
00692 }
00693
00694 const HepSymMatrix& ea = Ea();
00695 const HepVector& v_a = a();
00696
00697
00698 HepPoint3D pivot_work = pivot();
00699
00700 /*
00701 HepVector a_work = a();
00702 HepSymMatrix ea_work = Ea();
00703
00704 KalFitTrack helix_work(pivot_work, a_work, ea_work, 0, 0, 0);
00705 helix_work.pivot(ip);
00706
00707 HepVector a_temp = helix_work.a();
00708 HepSymMatrix ea_temp = helix_work.Ea();
00709
00710 seg.Ea_pre_bwd(ea_temp);
00711 seg.a_pre_bwd(a_temp);
00712
00713 */
00714
00715 seg.a_pre_bwd(a());
00716 seg.Ea_pre_bwd(Ea());
00717
00718 double dchi2R(99999999), dchi2L(99999999);
00719 HepVector v_H(5, 0);
00720 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
00721 v_H[3] = -meas.z();
00722 HepMatrix v_HT = v_H.T();
00723
00724 double estim = (v_HT * v_a)[0];
00725 HepVector ea_v_H = ea * v_H;
00726 HepMatrix ea_v_HT = (ea_v_H).T();
00727 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
00728 HepSymMatrix eaNew(5);
00729 HepVector aNew(5);
00730 double time = (*HitMdc).tdc();
00731
00732 //seg.dt(time);
00733 // if (Tof_correc_)
00734 // {
00735 // time -= tofest;
00736 // seg.dt(time);
00737 // }
00738 // double dd = 1.0e-4 * 40.0 * time;
00739 // seg.dd(dd);
00740
00741 double drifttime = getDriftTime(*HitMdc , tofest);
00742 seg.dt(drifttime);
00743 double ddl = getDriftDist(*HitMdc, drifttime, 0 );
00744 double ddr = getDriftDist(*HitMdc, drifttime, 1 );
00745 double erddl = getSigma( *HitMdc, a()[4], 0, ddl);
00746 double erddr = getSigma( *HitMdc, a()[4], 1, ddr);
00747
00748 double dmeas2[2] = { 0.1*ddl, 0.1*ddr }; // millimeter to centimeter
00749 double er_dmeas2[2] = {0., 0.};
00750 if(resolflag_ == 1) {
00751 er_dmeas2[0] = 0.1*erddl;
00752 er_dmeas2[1] = 0.1*erddr;
00753 }else if(resolflag_ == 0)
00754 {
00755 }
00756
00757 // Left hypothesis :
00758 if (lr == -1) {
00759 double er_dmeasL, dmeasL;
00760 if(Tof_correc_) {
00761 dmeasL = (-1.0)*dmeas2[0];
00762 er_dmeasL = er_dmeas2[0];
00763 } else {
00764 dmeasL = (-1.0)*fabs((*HitMdc).dist()[0]);
00765 er_dmeasL = (*HitMdc).erdist()[0];
00766 }
00767
00768
00769 //if(layerid < 4) er_dmeasL*=2.0;
00770
00771 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
00772 eaNew.assign(ea - ea_v_H * AL * ea_v_HT);
00773 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
00774 if(0. == RkL) RkL = 1.e-4;
00775
00776 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
00777 aNew = v_a + diffL;
00778 double sigL = (dmeasL - (v_HT * aNew)[0]);
00779 dchi2L = (sigL*sigL) / (RkL * er_dmeasL*er_dmeasL);
00780 } else if (lr == 1) {
00781 // Right hypothesis :
00782
00783 double er_dmeasR, dmeasR;
00784 if(Tof_correc_) {
00785 dmeasR = dmeas2[1];
00786 er_dmeasR = er_dmeas2[1];
00787 } else {
00788 dmeasR = fabs((*HitMdc).dist()[1]);
00789 er_dmeasR = (*HitMdc).erdist()[1];
00790 }
00791
00792
00793 //if(layerid < 4) er_dmeasR*=2.0;
00794
00795
00796 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
00797 eaNew.assign(ea - ea_v_H * AR * ea_v_HT);
00798 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
00799 if(0. == RkR) RkR = 1.e-4;
00800
00801 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
00802 aNew = v_a + diffR;
00803 double sigR = (dmeasR- (v_HT * aNew)[0]);
00804 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
00805 }
00806
00807 // Update Kalman result :
00808 #ifdef YDEBUG
00809 cout<<"in smoother_Mdc: lr= "<<lr<<" a: "<<v_a<<" a_NEW: "<<aNew<<endl;
00810 #endif
00811 double dchi2_loc(0);
00812 if ((dchi2R < dchi2cuts_calib[layerid] && lr == 1) ||
00813 (dchi2L < dchi2cuts_calib[layerid] && lr == -1)) {
00814
00815 ndf_back_++;
00816 flg = 1;
00817 int chge(1);
00818 if (!(aNew[2] && fabs(aNew[2]-a()[2]) < 1.0)) chge = 0;
00819 if (lr == 1) {
00820 chiSq_back_ += dchi2R;
00821 dchi2_loc = dchi2R;
00822 dd = 0.1*ddr;
00823 // if(debug_ ==4) std::cout<<"in smoother "<<std::endl;
00824
00825 } else if (lr == -1) {
00826 chiSq_back_ += dchi2L;
00827 dchi2_loc = dchi2L;
00828 dd = -0.1*ddl;
00829
00830 }
00831 if (chge){
00832 a(aNew);
00833 Ea(eaNew);
00834 }
00835
00836 seg.a_filt_bwd(aNew);
00837 seg.Ea_filt_bwd(eaNew);
00838
00839 HepVector a_pre_fwd_temp=seg.a_pre_fwd();
00840 if( (a_pre_fwd_temp[1]-seg.a_pre_bwd()[1]) > 3. * M_PI /2.) a_pre_fwd_temp[1]-= M_PI2;
00841 if( (a_pre_fwd_temp[1]-seg.a_pre_bwd()[1]) < -3. * M_PI /2. ) a_pre_fwd_temp[1]+= M_PI2;
00842
00843 seg.a_pre_fwd(a_pre_fwd_temp);
00844
00845 HepVector a_pre_filt_temp=seg.a_filt_fwd();
00846 if( (a_pre_filt_temp[1]-seg.a_pre_bwd()[1]) > 3. * M_PI /2. ) a_pre_filt_temp[1] -= M_PI2;
00847 if( (a_pre_filt_temp[1]-seg.a_pre_bwd()[1]) < -3. * M_PI /2.) a_pre_filt_temp[1] += M_PI2;
00848 seg.a_filt_fwd(a_pre_filt_temp);
00849
00850 /*
00851 KalFitTrack helixNew(pivot_work, aNew, eaNew, 0, 0, 0);
00852 helixNew.pivot(ip);
00853 a_temp = helixNew.a();
00854 ea_temp = helixNew.Ea();
00855 seg.Ea_filt_bwd(ea_temp);
00856 seg.a_filt_bwd(a_temp);
00857 */
00858
00859 int inv;
00860
00861 if(debug_ == 4){
00862 std::cout<<"seg.Ea_pre_bwd().inverse(inv) ..."<<seg.Ea_pre_bwd().inverse(inv)<<std::endl;
00863 std::cout<<"seg.Ea_pre_fwd().inverse(inv) ..."<<seg.Ea_pre_fwd().inverse(inv)<<std::endl;
00864 }
00865
00866 HepSymMatrix Ea_pre_comb = (seg.Ea_pre_bwd().inverse(inv)+seg.Ea_pre_fwd().inverse(inv)).inverse(inv);
00867 seg.Ea_exclude(Ea_pre_comb);
00868
00869
00870 if(debug_ == 4) {
00871 std::cout<<"Ea_pre_comb_ ... "<<Ea_pre_comb<<std::endl;
00872 std::cout<<"seg.a_pre_bwd()..."<<seg.a_pre_bwd()<<std::endl;
00873 std::cout<<"seg.a_pre_fwd()..."<<seg.a_pre_fwd()<<std::endl;
00874 }
00875 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());
00876 seg.a_exclude(a_pre_comb);
00877
00878 if(debug_ == 4) {
00879 std::cout<<"seg.Ea_filt_fwd()..."<<seg.Ea_filt_fwd()<<std::endl;
00880 std::cout<<"seg.Ea_filt_fwd().inverse(inv)..."<<seg.Ea_filt_fwd().inverse(inv)<<std::endl;
00881 }
00882 seg.Ea((seg.Ea_filt_fwd().inverse(inv)+seg.Ea_pre_bwd().inverse(inv)).inverse(inv));
00883 seg.Ea_include((seg.Ea_filt_fwd().inverse(inv)+seg.Ea_pre_bwd().inverse(inv)).inverse(inv));
00884
00885 if(debug_ == 4) {
00886 std::cout<<"seg.Ea() is ..."<<seg.Ea()<<std::endl;
00887 std::cout<<"seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd() ..."<<seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()<<std::endl;
00888 std::cout<<"seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd() ... "<<seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()<<std::endl;
00889 }
00890
00891 seg.a(seg.Ea()*(seg.Ea_filt_fwd().inverse(inv)*seg.a_filt_fwd()+seg.Ea_pre_bwd().inverse(inv)*seg.a_pre_bwd()));
00892 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()));
00893
00894 if(inv != 0) {
00895 //std::cout<<"ERROR OCCUR WHEN INVERSE MATRIX !"<<std::endl;
00896 }
00897
00898 seg.residual_exclude(dd - (v_HT*a_pre_comb)[0]);
00899 seg.residual_include(dd - (v_HT*seg.a())[0]);
00900 seg.doca_exclude(fabs((v_HT*a_pre_comb)[0]));
00901 seg.doca_include(fabs((v_HT*seg.a())[0]));
00902
00903 if(debug_ == 4){
00904 std::cout<<"the dd in smoother_Mdc is "<<dd
00905 <<" the (v_HT*a_pre_comb)[0] is "<<(v_HT*a_pre_comb)[0]<<std::endl;
00906 }
00907
00908 //compare the two method to calculate the include doca value :
00909 if(debug_ == 4){
00910 std::cout<<"method 1 ..."<<sqrt(seg.a()[0]*seg.a()[0]+seg.a()[3]*seg.a()[3])<<std::endl;
00911 std::cout<<"method 2 ..."<<fabs((v_HT*seg.a())[0])<<std::endl;
00912 }
00913
00914
00916 KalFitTrack helixNew1(pivot_work, seg.a(), seg.Ea(), 0, 0, 0);
00917 helixNew1.pivot(ip);
00918 HepVector a_temp1 = helixNew1.a();
00919 HepSymMatrix ea_temp1 = helixNew1.Ea();
00920 seg.Ea(ea_temp1);
00921 seg.a(a_temp1);
00922 seg.Ea_include(ea_temp1);
00923 seg.a_include(a_temp1);
00924
00925 KalFitTrack helixNew2(pivot_work, seg.a_exclude(), seg.Ea_exclude(), 0, 0, 0);
00926 helixNew2.pivot(ip);
00927 HepVector a_temp2 = helixNew2.a();
00928 HepSymMatrix ea_temp2 = helixNew2.Ea();
00929 seg.Ea_exclude(ea_temp2);
00930 seg.a_exclude(a_temp2);
00931
00932 seg.tof(tofest);
00933 seg.dd(dd);
00934
00935 }
00936 return chiSq_back_;
00937 }
|
|
|
|
|
|
00242 {
00243 return m_ac[4];
00244 }
|
|
|
00191 { return tof_; }
|
|
|
Update the tof estimation.
|
|
|
00191 { return tof_; }
|
|
|
Update the tof estimation.
01300 {
01301 double light_speed( 29.9792458 ); // light speed in cm/nsec
01302 double t = tanl();
01303 double pmag( sqrt( 1.0 + t*t ) / kappa());
01304 if (pmag !=0) {
01305 double mass_over_p( mass_ / pmag );
01306 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
01307 tof_ += path / ( light_speed * beta );
01308 }
01309
01310 if (tofall_) {
01311 if (p_kaon_ > 0){
01312 double massk_over_p( MASS[3] / p_kaon_ );
01313 double beta_kaon( 1.0 / sqrt( 1.0 + massk_over_p * massk_over_p ) );
01314 tof_kaon_ += path / (light_speed * beta_kaon);
01315 }
01316 if (p_proton_ > 0){
01317 double massp_over_p( MASS[4] / p_proton_ );
01318 double beta_proton( 1.0 / sqrt( 1.0 + massp_over_p * massp_over_p ) );
01319 tof_proton_ += path / (light_speed * beta_proton);
01320 }
01321 }
01322 }
|
|
|
00192 { return tof_kaon_; }
|
|
|
00192 { return tof_kaon_; }
|
|
|
00193 { return tof_proton_; }
|
|
|
00193 { return tof_proton_; }
|
|
|
00208 { trasan_id_ = t;}
|
|
|
00180 { return trasan_id_; }
|
|
|
00208 { trasan_id_ = t;}
|
|
|
00180 { return trasan_id_; }
|
|
|
Reinitialize (modificator).
00207 { type_ = t;}
|
|
|
00179 { return type_; }
|
|
|
Reinitialize (modificator).
00207 { type_ = t;}
|
|
|
00179 { return type_; }
|
|
|
|
|
|
01795 {
01796 int j(0);
01797 if (i < 31){
01798 j = (int) pow(2.,i);
01799 if (!(pat1_ & j))
01800 pat1_ = pat1_ | j;
01801 } else if (i < 50) {
01802 j = (int) pow(2.,(i-31));
01803 if (!(pat2_ & j))
01804 pat2_ = pat2_ | j;
01805 }
01806 }
|
|
|
|
|
|
00116 {
00117 pivot_forMdc_ = pivot();
00118 a_forMdc_ = a();
00119 Ea_forMdc_ = Ea();
00120 }
|
|
||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||||||||||||||
|
Include the Mdc wire hits.
|
|
||||||||||||||||||||||||||||
|
02446 {
02447
02448
02449 HepPoint3D ip(0.,0.,0.);
02450
02451 KalFitHitMdc* HitMdc = HelixSeg.HitMdc();
02452 double lr = HitMdc->LR();
02453 int layerid = (*HitMdc).wire().layer().layerId();
02454 double entrangle = HitMdc->rechitptr()->getEntra();
02455
02456 if(debug_ == 4) {
02457 std::cout<<"in update_hits(kalFitHelixSeg,...) : the layerid ="<<layerid<<std::endl;
02458 std::cout<<" update_hits: lr= " << lr <<std::endl;
02459 }
02460 // For taking the propagation delay into account :
02461 const KalFitWire& Wire = HitMdc->wire();
02462 int wire_ID = Wire.geoID();
02463
02464 if (wire_ID<0 || wire_ID>6796){ //bes
02465 std::cout << " KalFitTrack: wire_ID problem : " << wire_ID
02466 << std::endl;
02467 return 0;
02468 }
02469 int flag_ster = Wire.stereo();
02470 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
02471 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
02472 double tofest(0);
02473 double phi = fmod(phi0() + M_PI4, M_PI2);
02474 double csf0 = cos(phi);
02475 double snf0 = (1. - csf0) * (1. + csf0);
02476 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
02477 if(phi > M_PI) snf0 = - snf0;
02478 if (Tof_correc_){
02479 double t = tanl();
02480 double dphi(0);
02481
02482 Helix work = *(Helix*)this;
02483 work.ignoreErrorMatrix();
02484 double phi_ip(0);
02485 Hep3Vector ip(0, 0, 0);
02486 work.pivot(ip);
02487 phi_ip = work.phi0();
02488 if (fabs(phi - phi_ip) > M_PI) {
02489 if (phi > phi_ip) phi -= 2 * M_PI;
02490 else phi_ip -= 2 * M_PI;
02491 }
02492 dphi = phi - phi_ip;
02493
02494 double l = fabs(radius() * dphi * sqrt(1 + t * t));
02495 double pmag( sqrt( 1.0 + t*t ) / kappa());
02496 double mass_over_p( mass_ / pmag );
02497 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
02498 tofest = l / ( 29.9792458 * beta );
02499 if(csmflag==1 && (*HitMdc).wire().y()>0.) tofest= -1. * tofest;
02500 }
02501
02502 const HepSymMatrix& ea = Ea();
02503 HelixSeg.Ea_pre_fwd(ea);
02504 HelixSeg.Ea_filt_fwd(ea);
02505
02506
02507 const HepVector& v_a = a();
02508 HelixSeg.a_pre_fwd(v_a);
02509 HelixSeg.a_filt_fwd(v_a);
02510
02511 /*
02512
02513 HepPoint3D pivot_work = pivot();
02514 HepVector a_work = a();
02515 HepSymMatrix ea_work = Ea();
02516 KalFitTrack helix_work(pivot_work, a_work, ea_work, 0, 0, 0);
02517 helix_work.pivot(ip);
02518
02519 HepVector a_temp = helix_work.a();
02520 HepSymMatrix ea_temp = helix_work.Ea();
02521
02522 HelixSeg.Ea_pre_fwd(ea_temp);
02523 HelixSeg.a_pre_fwd(a_temp);
02524
02525 // My God, for the protection purpose
02526 HelixSeg.Ea_filt_fwd(ea_temp);
02527 HelixSeg.a_filt_fwd(a_temp);
02528
02529 */
02530
02531
02532 double dchi2R(99999999), dchi2L(99999999);
02533
02534 HepVector v_H(5, 0);
02535 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
02536 v_H[3] = -meas.z();
02537 HepMatrix v_HT = v_H.T();
02538
02539 double estim = (v_HT * v_a)[0];
02540 HepVector ea_v_H = ea * v_H;
02541 HepMatrix ea_v_HT = (ea_v_H).T();
02542 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
02543
02544 HepSymMatrix eaNewL(5), eaNewR(5);
02545 HepVector aNewL(5), aNewR(5);
02546 #ifdef YDEBUG
02547 cout<<"phi "<<phi<<" snf0 "<<snf0<<" csf0 "<<csf0<<endl
02548 <<" meas: "<<meas <<endl;
02549 cout<<"the related matrix:"<<endl;
02550 cout<<"v_H "<<v_H<<endl;
02551 cout<<"v_a "<<v_a<<endl;
02552 cout<<"ea "<<ea<<endl;
02553 cout<<"v_H_T_ea_v_H "<<v_H_T_ea_v_H<<endl;
02554 cout<<"LR_= "<< LR_ << "lr= " << lr <<endl;
02555 #endif
02556
02557 double time = (*HitMdc).tdc();
02558 // if (Tof_correc_)
02559 // time = time - way*tofest;
02560
02561 // double dd = 1.0e-4 * 40.0 * time;
02562 //the following 3 line : tempory
02563
02564 double drifttime = getDriftTime(*HitMdc , tofest);
02565 double ddl = getDriftDist(*HitMdc, drifttime, 0 );
02566 double ddr = getDriftDist(*HitMdc, drifttime, 1 );
02567 double erddl = getSigma( *HitMdc, a()[4], 0, ddl);
02568 double erddr = getSigma( *HitMdc, a()[4], 1, ddr);
02569
02570 if(debug_ == 4){
02571 std::cout<<"the pure drift time is "<<drifttime<<std::endl;
02572 std::cout<<"the drift dist left hypothesis is "<<ddl<<std::endl;
02573 std::cout<<"the drift dist right hypothesis is "<<ddr<<std::endl;
02574 std::cout<<"the error sigma left hypothesis is "<<erddl<<std::endl;
02575 std::cout<<"the error sigma right hypothesis is "<<erddr<<std::endl;
02576 }
02577 double dmeas2[2] = { 0.1*ddl , 0.1*ddr }; //millimeter to centimeter
02578 double er_dmeas2[2] = {0. , 0.} ;
02579
02580 if(resolflag_ == 1) {
02581 er_dmeas2[0] = 0.1*erddl;
02582 er_dmeas2[1] = 0.1*erddr;
02583 }
02584 else if(resolflag_ == 0) {
02585 // int layid = (*HitMdc).wire().layer().layerId();
02586 // double sigma = getSigma(layid, dd);
02587 // er_dmeas2[0] = er_dmeas2[1] = sigma;
02588 }
02589
02590 // Left hypothesis :
02591 if ((LR_==0 && lr != 1.0) ||
02592 (LR_==1 && lr == -1.0)) {
02593
02594 double er_dmeasL, dmeasL;
02595 if(Tof_correc_) {
02596 dmeasL = (-1.0)*fabs(dmeas2[0]);
02597 er_dmeasL = er_dmeas2[0];
02598 } else {
02599 dmeasL = (-1.0)*fabs((*HitMdc).dist()[0]);
02600 er_dmeasL = (*HitMdc).erdist()[0];
02601 }
02602
02603 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
02604 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
02605 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
02606 if(0. == RkL) RkL = 1.e-4;
02607
02608 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
02609
02610 aNewL = v_a + diffL;
02611 double sigL = dmeasL -(v_HT * aNewL)[0];
02612 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
02613 }
02614
02615 if ((LR_==0 && lr != -1.0) ||
02616 (LR_==1 && lr == 1.0)) {
02617
02618 double er_dmeasR, dmeasR;
02619 if(Tof_correc_) {
02620 dmeasR = dmeas2[1];
02621 er_dmeasR = er_dmeas2[1];
02622 } else {
02623 dmeasR = fabs((*HitMdc).dist()[1]);
02624 er_dmeasR = (*HitMdc).erdist()[1];
02625 }
02626
02627 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
02628 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
02629 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
02630 if(0. == RkR) RkR = 1.e-4;
02631
02632 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
02633
02634 aNewR = v_a + diffR;
02635 double sigR = dmeasR -(v_HT * aNewR)[0];
02636 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
02637 }
02638
02639 // Update Kalman result :
02640 double dchi2_loc(0);
02641
02642 #ifdef YDEBUG
02643 cout<<" track::update_hits......"<<endl;
02644 std::cout << " dchi2R = " << dchi2R << ", dchi2L = " << dchi2L
02645 << " estimate = "<<estim<< std::endl;
02646 std::cout << " dmeasR = " << dmeas2[1]
02647 << ", dmeasL = " << (-1.0)*fabs(dmeas2[0]) << ", check HitMdc.ddl = "
02648 << (*HitMdc).dist()[0] << std::endl;
02649 #endif
02650
02651 if (fabs(dchi2R-dchi2L)<10. && inext>0) {
02652
02653 KalFitHitMdc & HitMdc_next = HitsMdc_[inext];
02654
02655 Helix H_fromR(pivot(), aNewR, eaNewR);
02656 double chi2_fromR(chi2_next(H_fromR, HitMdc_next,csmflag));
02657
02658 Helix H_fromL(pivot(), aNewL, eaNewL);
02659 double chi2_fromL(chi2_next(H_fromL, HitMdc_next,csmflag));
02660 #ifdef YDEBUG
02661 std::cout << " chi2_fromL = " << chi2_fromL
02662 << ", chi2_fromR = " << chi2_fromR << std::endl;
02663 #endif
02664 if (fabs(chi2_fromR-chi2_fromL)<10.){
02665 int inext2(-1);
02666 if (inext+1<HitsMdc_.size())
02667 for( int k=inext+1 ; k < HitsMdc_.size(); k++ )
02668 if (!(HitsMdc_[k].chi2()<0)) {
02669 inext2 = k;
02670 break;
02671 }
02672
02673 if (inext2 != -1){
02674 KalFitHitMdc & HitMdc_next2 = HitsMdc_[inext2];
02675 double chi2_fromR2(chi2_next(H_fromR, HitMdc_next2, csmflag));
02676 double chi2_fromL2(chi2_next(H_fromL, HitMdc_next2, csmflag));
02677 #ifdef YDEBUG
02678 std::cout << " chi2_fromL2 = " << chi2_fromL2
02679 << ", chi2_fromR2 = " << chi2_fromR2 << std::endl;
02680 #endif
02681 if (fabs(dchi2R+chi2_fromR+chi2_fromR2-(dchi2L+chi2_fromL+chi2_fromL2))<2) {
02682 if (chi2_fromR2<chi2_fromL2)
02683 dchi2L = DBL_MAX;
02684 else
02685 dchi2R = DBL_MAX;
02686 }
02687 }
02688 }
02689
02690 if (!(dchi2L == DBL_MAX && dchi2R == DBL_MAX)) {
02691 if (dchi2R+chi2_fromR < dchi2L+chi2_fromL){
02692 dchi2L = DBL_MAX;
02693 #ifdef YDEBUG
02694 std::cout << " choose right..." << std::endl;
02695 #endif
02696 } else {
02697 dchi2R = DBL_MAX;
02698 #ifdef YDEBUG
02699 std::cout << " choose left..." << std::endl;
02700 #endif
02701 }
02702 }
02703 }
02704
02705 if ((0 < dchi2R && dchi2R < dchi2cutf_calib[layerid]) ||
02706 (0 < dchi2L && dchi2L < dchi2cutf_calib[layerid])) {
02707
02708 if (((LR_==0 && dchi2R < dchi2L && lr != -1) ||
02709 (LR_==1 && lr == 1)) &&
02710 fabs(aNewR[2]-a()[2]) < 1000. && aNewR[2]) {
02711 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2R < dchi2cutf_calib[layerid]){
02712 nchits_++;
02713 if (flag_ster) nster_++;
02714
02715 Ea(eaNewR);
02716 HelixSeg.Ea_filt_fwd(eaNewR);
02717 a(aNewR);
02718 HelixSeg.a_filt_fwd(aNewR);
02719
02720 /*
02721 Ea(eaNewR);
02722 a(aNewR);
02723
02724 KalFitTrack helixR(pivot_work, aNewR, eaNewR, 0, 0, 0);
02725 helixR.pivot(ip);
02726
02727 a_temp = helixR.a();
02728 ea_temp = helixR.Ea();
02729
02730 HelixSeg.Ea_filt_fwd(ea_temp);
02731 HelixSeg.a_filt_fwd(a_temp);
02732 //HelixSeg.filt_include(1);
02733
02734 */
02735
02736 chiSq_ += dchi2R;
02737 dchi2_loc = dchi2R;
02738 if (l_mass_ == lead_){
02739 (*HitMdc).LR(1);
02740 HelixSeg.LR(1);
02741 }
02742 update_bit((*HitMdc).wire().layer().layerId());
02743 }
02744 } else if (((LR_==0 && dchi2L <= dchi2R && lr != 1) ||
02745 (LR_==1 && lr == -1)) &&
02746 fabs(aNewL[2]-a()[2]) < 1000. && aNewL[2]){
02747 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2L < dchi2cutf_calib[layerid]){
02748 nchits_++;
02749 if (flag_ster) nster_++;
02750 Ea(eaNewL);
02751 HelixSeg.Ea_filt_fwd(eaNewL);
02752 a(aNewL);
02753 HelixSeg.a_filt_fwd(aNewL);
02754
02755
02756 /*
02757 Ea(eaNewL);
02758 a(aNewL);
02759
02760 KalFitTrack helixL(pivot_work, aNewL, eaNewL, 0, 0, 0);
02761 helixL.pivot(ip);
02762 a_temp = helixL.a();
02763 ea_temp = helixL.Ea();
02764 HelixSeg.Ea_filt_fwd(ea_temp);
02765 HelixSeg.a_filt_fwd(a_temp);
02766 //HelixSeg.filt_include(1);
02767
02768 */
02769
02770 chiSq_ += dchi2L;
02771 dchi2_loc = dchi2L;
02772 if (l_mass_ == lead_){
02773 (*HitMdc).LR(-1);
02774 HelixSeg.LR(-1);
02775 }
02776 update_bit((*HitMdc).wire().layer().layerId());
02777 }
02778 }
02779 }
02780
02781 if (dchi2_loc > dchi2_max_) {
02782 dchi2_max_ = dchi2_loc ;
02783 r_max_ = pivot().perp();
02784 }
02785 dchi2out = dchi2_loc;
02786 // if(dchi2out == 0) {
02787 // dchi2out = ( (dchi2L < dchi2R ) ? dchi2L : dchi2R ) ;
02788 // }
02789 return chiSq_;
02790 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Include the Mdc wire hits.
01834 {
01835
01836 double lr = HitMdc.LR();
01837 //std::cout<<" in update_hits: lr= " << lr <<std::endl;
01838 // For taking the propagation delay into account :
01839 const KalFitWire& Wire = HitMdc.wire();
01840 int wire_ID = Wire.geoID();
01841 int layerid = HitMdc.wire().layer().layerId();
01842 //std::cout<<" when in layer: "<<layerid<<std::endl;
01843 double entrangle = HitMdc.rechitptr()->getEntra();
01844
01845 if (wire_ID<0 || wire_ID>6796){ //bes
01846 std::cout << " KalFitTrack: wire_ID problem : " << wire_ID
01847 << std::endl;
01848 return 0;
01849 }
01850 int flag_ster = Wire.stereo();
01851 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
01852 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
01853 double tofest(0);
01854 double phi = fmod(phi0() + M_PI4, M_PI2);
01855 double csf0 = cos(phi);
01856 double snf0 = (1. - csf0) * (1. + csf0);
01857 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
01858 if(phi > M_PI) snf0 = - snf0;
01859 if (Tof_correc_){
01860 double t = tanl();
01861 double dphi(0);
01862
01863 Helix work = *(Helix*)this;
01864 work.ignoreErrorMatrix();
01865 double phi_ip(0);
01866 Hep3Vector ip(0, 0, 0);
01867 work.pivot(ip);
01868 phi_ip = work.phi0();
01869 if (fabs(phi - phi_ip) > M_PI) {
01870 if (phi > phi_ip) phi -= 2 * M_PI;
01871 else phi_ip -= 2 * M_PI;
01872 }
01873 dphi = phi - phi_ip;
01874 double l = fabs(radius() * dphi * sqrt(1 + t * t));
01875 double pmag( sqrt( 1.0 + t*t ) / kappa());
01876 double mass_over_p( mass_ / pmag );
01877 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
01878 tofest = l / ( 29.9792458 * beta );
01879 if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
01880 }
01881
01882 const HepSymMatrix& ea = Ea();
01883 const HepVector& v_a = a();
01884 double dchi2R(99999999), dchi2L(99999999);
01885
01886 HepVector v_H(5, 0);
01887 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
01888 v_H[3] = -meas.z();
01889 HepMatrix v_HT = v_H.T();
01890
01891 double estim = (v_HT * v_a)[0];
01892 dtrack = estim;
01893 // std::cout<<" in update_hits estim is "<<estim<<std::endl;
01894 HepVector ea_v_H = ea * v_H;
01895 HepMatrix ea_v_HT = (ea_v_H).T();
01896 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
01897
01898 HepSymMatrix eaNewL(5), eaNewR(5);
01899 HepVector aNewL(5), aNewR(5);
01900
01901 double drifttime = getDriftTime(HitMdc , tofest);
01902 double ddl = getDriftDist(HitMdc, drifttime, 0 );
01903 double ddr = getDriftDist(HitMdc, drifttime, 1 );
01904 double erddl = getSigma( HitMdc, a()[4], 0, ddl);
01905 double erddr = getSigma( HitMdc, a()[4], 1, ddr);
01906
01907 if(debug_ == 4) {
01908 std::cout<<"drifttime in update_hits() for ananlysis is "<<drifttime<<std::endl;
01909 std::cout<<"ddl in update_hits() for ananlysis is "<<ddl<<std::endl;
01910 std::cout<<"ddr in update_hits() for ananlysis is "<<ddr<<std::endl;
01911 std::cout<<"erddl in update_hits() for ananlysis is "<<erddl<<std::endl;
01912 std::cout<<"erddr in update_hits() for ananlysis is "<<erddr<<std::endl;
01913 }
01914
01915 //the following 3 line : tempory
01916 double dmeas2[2] = { 0.1*ddl, 0.1*ddr }; // millimeter to centimeter
01917 double er_dmeas2[2] = {0.,0.};
01918 if(resolflag_ == 1) {
01919 er_dmeas2[0] = 0.1*erddl;
01920 er_dmeas2[1] = 0.1*erddr;
01921 }
01922 else if(resolflag_ == 0) {
01923 // int layid = HitMdc.wire().layer().layerId();
01924 // double sigma = getSigma(layid, dd);
01925 // er_dmeas2[0] = er_dmeas2[1] = sigma;
01926 }
01927
01928 // Left hypothesis :
01929 if ((LR_==0 && lr != 1.0) ||
01930 (LR_==1 && lr == -1.0)) {
01931 double er_dmeasL, dmeasL;
01932 if(Tof_correc_) {
01933 dmeasL = (-1.0)*fabs(dmeas2[0]);
01934 er_dmeasL = er_dmeas2[0];
01935 } else {
01936 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
01937 er_dmeasL = HitMdc.erdist()[0];
01938 }
01939 dtdc = dmeasL;
01940 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
01941 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
01942 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
01943 if(debug_ == 4){
01944 std::cout<<" ea_v_H * AL * ea_v_HT is: "<<ea_v_H * AL * ea_v_HT<<std::endl;
01945 std::cout<<" v_H is: "<<v_H<<" Ea is: "<<ea<<" v_H_T_ea_v_H is: "<<v_H_T_ea_v_H<<std::endl;
01946 std::cout<<" AL is: "<<AL<<" (v_H_T_ea_v_H)[0] * AL is: "<<(v_H_T_ea_v_H)[0]*AL<<std::endl;
01947 }
01948
01949 if(0. == RkL) RkL = 1.e-4;
01950
01951 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
01952 aNewL = v_a + diffL;
01953 double sigL = dmeasL -(v_HT * aNewL)[0];
01954 dtracknew = (v_HT * aNewL)[0];
01955 dchi2L = (sigL * sigL)/(RkL * er_dmeasL*er_dmeasL);
01956
01957 if(debug_ == 4){
01958 std::cout<<" fwd_filter : in left hypothesis dchi2L is "<<dchi2L<<std::endl;
01959 }
01960 }
01961
01962 if ((LR_==0 && lr != -1.0) ||
01963 (LR_==1 && lr == 1.0)) {
01964 double er_dmeasR, dmeasR;
01965 if(Tof_correc_) {
01966 dmeasR = dmeas2[1];
01967 er_dmeasR = er_dmeas2[1];
01968 } else {
01969 dmeasR = fabs(HitMdc.dist()[1]);
01970 er_dmeasR = HitMdc.erdist()[1];
01971 }
01972 dtdc = dmeasR;
01973 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
01974 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
01975 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
01976
01977 if(debug_ == 4){
01978 std::cout<<" ea_v_H * AR * ea_v_HT is: "<<ea_v_H * AR * ea_v_HT<<std::endl;
01979 std::cout<<" v_H is: "<<v_H<<" Ea is: "<<ea<<" v_H_T_ea_v_H is: "<<v_H_T_ea_v_H<<std::endl;
01980 std::cout<<" AR is: "<<AR<<" (v_H_T_ea_v_H)[0] * AR is: "<<(v_H_T_ea_v_H)[0]*AR<<std::endl;
01981 }
01982
01983 if(0. == RkR) RkR = 1.e-4;
01984 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
01985 aNewR = v_a + diffR;
01986 double sigR = dmeasR -(v_HT * aNewR)[0];
01987 dtracknew = (v_HT * aNewR)[0];
01988 dchi2R = (sigR*sigR)/(RkR * er_dmeasR*er_dmeasR);
01989 if(debug_ == 4){
01990 std::cout<<" fwd_filter : in right hypothesis dchi2R is "<<dchi2R<<std::endl;
01991 }
01992 }
01993 // Update Kalman result :
01994 double dchi2_loc(0);
01995 if (fabs(dchi2R-dchi2L)<10. && inext>0) {
01996 KalFitHitMdc & HitMdc_next = HitsMdc_[inext];
01997 Helix H_fromR(pivot(), aNewR, eaNewR);
01998 double chi2_fromR(chi2_next(H_fromR, HitMdc_next, csmflag));
01999 Helix H_fromL(pivot(), aNewL, eaNewL);
02000 double chi2_fromL(chi2_next(H_fromL, HitMdc_next, csmflag));
02001 #ifdef YDEBUG
02002 std::cout << " chi2_fromL = " << chi2_fromL
02003 << ", chi2_fromR = " << chi2_fromR << std::endl;
02004 #endif
02005 if (fabs(chi2_fromR-chi2_fromL)<10.){
02006 int inext2(-1);
02007 if (inext+1<HitsMdc_.size())
02008 for( int k=inext+1 ; k < HitsMdc_.size(); k++ )
02009 if (!(HitsMdc_[k].chi2()<0)) {
02010 inext2 = k;
02011 break;
02012 }
02013
02014 if (inext2 != -1){
02015 KalFitHitMdc & HitMdc_next2 = HitsMdc_[inext2];
02016 double chi2_fromR2(chi2_next(H_fromR, HitMdc_next2, csmflag));
02017 double chi2_fromL2(chi2_next(H_fromL, HitMdc_next2, csmflag));
02018 #ifdef YDEBUG
02019 std::cout << " chi2_fromL2 = " << chi2_fromL2
02020 << ", chi2_fromR2 = " << chi2_fromR2 << std::endl;
02021 #endif
02022 if (fabs(dchi2R+chi2_fromR+chi2_fromR2-(dchi2L+chi2_fromL+chi2_fromL2))<2) {
02023 if (chi2_fromR2<chi2_fromL2)
02024 dchi2L = DBL_MAX;
02025 else
02026 dchi2R = DBL_MAX;
02027 }
02028 }
02029 }
02030
02031 if (!(dchi2L == DBL_MAX && dchi2R == DBL_MAX)) {
02032 if (dchi2R+chi2_fromR < dchi2L+chi2_fromL){
02033 dchi2L = DBL_MAX;
02034 #ifdef YDEBUG
02035 std::cout << " choose right..." << std::endl;
02036 #endif
02037 } else {
02038 dchi2R = DBL_MAX;
02039 #ifdef YDEBUG
02040 std::cout << " choose left..." << std::endl;
02041 #endif
02042 }
02043 }
02044 }
02045 if ((0 < dchi2R && dchi2R < dchi2cutf_anal[layerid]) ||
02046 (0 < dchi2L && dchi2L < dchi2cutf_anal[layerid])) {
02047
02048 if (((LR_==0 && dchi2R < dchi2L && lr != -1) ||
02049 (LR_==1 && lr == 1)) &&
02050 fabs(aNewR[2]-a()[2]) < 1000. && aNewR[2]) {
02051 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2R < dchi2cutf_anal[layerid]){
02052 nchits_++;
02053 if (flag_ster) nster_++;
02054 Ea(eaNewR);
02055 a(aNewR);
02056 chiSq_ += dchi2R;
02057 dchi2_loc = dchi2R;
02058 if (l_mass_ == lead_){
02059 HitMdc.LR(1);
02060 }
02061 update_bit(HitMdc.wire().layer().layerId());
02062 }
02063 } else if (((LR_==0 && dchi2L <= dchi2R && lr != 1) ||
02064 (LR_==1 && lr == -1)) &&
02065 fabs(aNewL[2]-a()[2]) < 1000. && aNewL[2]){
02066 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2L < dchi2cutf_anal[layerid]){
02067 nchits_++;
02068 if (flag_ster) nster_++;
02069 Ea(eaNewL);
02070 a(aNewL);
02071 chiSq_ += dchi2L;
02072 dchi2_loc = dchi2L;
02073 if (l_mass_ == lead_){
02074 HitMdc.LR(-1);
02075 }
02076 update_bit(HitMdc.wire().layer().layerId());
02077 }
02078 }
02079 }
02080 if (dchi2_loc > dchi2_max_) {
02081 dchi2_max_ = dchi2_loc ;
02082 r_max_ = pivot().perp();
02083 }
02084 dchi2out = dchi2_loc;
02085 //if(dchi2out == 0) {
02086 // dchi2out = ( (dchi2L < dchi2R ) ? dchi2L : dchi2R ) ;
02087 // }
02088 return chiSq_;
02089 }
|
|
||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||
|
02092 {
02093
02094
02095 HepPoint3D ip(0.,0.,0.);
02096
02097 KalFitHitMdc* HitMdc = HelixSeg.HitMdc();
02098 double lr = HitMdc->LR();
02099 int layerid = (*HitMdc).wire().layer().layerId();
02100 // cout<<"layerid: "<<layerid<<endl;
02101 double entrangle = HitMdc->rechitptr()->getEntra();
02102
02103 if(debug_ == 4) {
02104 std::cout<<"in update_hits(kalFitHelixSeg,...) : the layerid ="<<layerid<<std::endl;
02105 std::cout<<" update_hits: lr= " << lr <<std::endl;
02106 }
02107 // For taking the propagation delay into account :
02108 const KalFitWire& Wire = HitMdc->wire();
02109 int wire_ID = Wire.geoID();
02110
02111 if (wire_ID<0 || wire_ID>6796){ //bes
02112 std::cout << " KalFitTrack: wire_ID problem : " << wire_ID
02113 << std::endl;
02114 return 0;
02115 }
02116 int flag_ster = Wire.stereo();
02117 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
02118 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
02119 double tofest(0);
02120 double phi = fmod(phi0() + M_PI4, M_PI2);
02121 double csf0 = cos(phi);
02122 double snf0 = (1. - csf0) * (1. + csf0);
02123 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
02124 if(phi > M_PI) snf0 = - snf0;
02125 if (Tof_correc_){
02126 double t = tanl();
02127 double dphi(0);
02128
02129 Helix work = *(Helix*)this;
02130 work.ignoreErrorMatrix();
02131 double phi_ip(0);
02132 Hep3Vector ip(0, 0, 0);
02133 work.pivot(ip);
02134 phi_ip = work.phi0();
02135 if (fabs(phi - phi_ip) > M_PI) {
02136 if (phi > phi_ip) phi -= 2 * M_PI;
02137 else phi_ip -= 2 * M_PI;
02138 }
02139 dphi = phi - phi_ip;
02140
02141 double l = fabs(radius() * dphi * sqrt(1 + t * t));
02142 double pmag( sqrt( 1.0 + t*t ) / kappa());
02143 double mass_over_p( mass_ / pmag );
02144 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
02145 tofest = l / ( 29.9792458 * beta );
02146 if(csmflag==1 && (*HitMdc).wire().y()>0.) tofest= -1. * tofest;
02147 }
02148
02149 const HepSymMatrix& ea = Ea();
02150 HelixSeg.Ea_pre_fwd(ea);
02151 HelixSeg.Ea_filt_fwd(ea);
02152
02153
02154 const HepVector& v_a = a();
02155 HelixSeg.a_pre_fwd(v_a);
02156 HelixSeg.a_filt_fwd(v_a);
02157
02158 /*
02159
02160 HepPoint3D pivot_work = pivot();
02161 HepVector a_work = a();
02162 HepSymMatrix ea_work = Ea();
02163 KalFitTrack helix_work(pivot_work, a_work, ea_work, 0, 0, 0);
02164 helix_work.pivot(ip);
02165
02166 HepVector a_temp = helix_work.a();
02167 HepSymMatrix ea_temp = helix_work.Ea();
02168
02169 HelixSeg.Ea_pre_fwd(ea_temp);
02170 HelixSeg.a_pre_fwd(a_temp);
02171
02172 // My God, for the protection purpose
02173 HelixSeg.Ea_filt_fwd(ea_temp);
02174 HelixSeg.a_filt_fwd(a_temp);
02175
02176 */
02177
02178 double dchi2R(99999999), dchi2L(99999999);
02179
02180 HepVector v_H(5, 0);
02181 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
02182 v_H[3] = -meas.z();
02183 HepMatrix v_HT = v_H.T();
02184
02185 double estim = (v_HT * v_a)[0];
02186 HepVector ea_v_H = ea * v_H;
02187 HepMatrix ea_v_HT = (ea_v_H).T();
02188 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
02189
02190 HepSymMatrix eaNewL(5), eaNewR(5);
02191 HepVector aNewL(5), aNewR(5);
02192 #ifdef YDEBUG
02193 cout<<"phi "<<phi<<" snf0 "<<snf0<<" csf0 "<<csf0<<endl
02194 <<" meas: "<<meas <<endl;
02195 cout<<"the related matrix:"<<endl;
02196 cout<<"v_H "<<v_H<<endl;
02197 cout<<"v_a "<<v_a<<endl;
02198 cout<<"ea "<<ea<<endl;
02199 cout<<"v_H_T_ea_v_H "<<v_H_T_ea_v_H<<endl;
02200 cout<<"LR_= "<< LR_ << "lr= " << lr <<endl;
02201 #endif
02202
02203 double time = (*HitMdc).tdc();
02204 // if (Tof_correc_)
02205 // time = time - way*tofest;
02206
02207 // double dd = 1.0e-4 * 40.0 * time;
02208 //the following 3 line : tempory
02209
02210 double drifttime = getDriftTime(*HitMdc , tofest);
02211 double ddl = getDriftDist(*HitMdc, drifttime, 0 );
02212 double ddr = getDriftDist(*HitMdc, drifttime, 1 );
02213 double erddl = getSigma( *HitMdc, a()[4], 0, ddl);
02214 double erddr = getSigma( *HitMdc, a()[4], 1, ddr);
02215
02216 if(debug_ == 4){
02217 std::cout<<"the pure drift time is "<<drifttime<<std::endl;
02218 std::cout<<"the drift dist left hypothesis is "<<ddl<<std::endl;
02219 std::cout<<"the drift dist right hypothesis is "<<ddr<<std::endl;
02220 std::cout<<"the error sigma left hypothesis is "<<erddl<<std::endl;
02221 std::cout<<"the error sigma right hypothesis is "<<erddr<<std::endl;
02222 }
02223 double dmeas2[2] = { 0.1*ddl , 0.1*ddr }; //millimeter to centimeter
02224 double er_dmeas2[2] = {0. , 0.} ;
02225
02226 if(resolflag_ == 1) {
02227 er_dmeas2[0] = 0.1*erddl;
02228 er_dmeas2[1] = 0.1*erddr;
02229 }
02230 else if(resolflag_ == 0) {
02231 // int layid = (*HitMdc).wire().layer().layerId();
02232 // double sigma = getSigma(layid, dd);
02233 // er_dmeas2[0] = er_dmeas2[1] = sigma;
02234 }
02235
02236 // Left hypothesis :
02237 if ((LR_==0 && lr != 1.0) ||
02238 (LR_==1 && lr == -1.0)) {
02239
02240 double er_dmeasL, dmeasL;
02241 if(Tof_correc_) {
02242 dmeasL = (-1.0)*fabs(dmeas2[0]);
02243 er_dmeasL = er_dmeas2[0];
02244 } else {
02245 dmeasL = (-1.0)*fabs((*HitMdc).dist()[0]);
02246 er_dmeasL = (*HitMdc).erdist()[0];
02247 }
02248
02249 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
02250 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
02251 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
02252 if(0. == RkL) RkL = 1.e-4;
02253
02254 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
02255
02256 aNewL = v_a + diffL;
02257 double sigL = dmeasL -(v_HT * aNewL)[0];
02258 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
02259 }
02260
02261 if ((LR_==0 && lr != -1.0) ||
02262 (LR_==1 && lr == 1.0)) {
02263
02264 double er_dmeasR, dmeasR;
02265 if(Tof_correc_) {
02266 dmeasR = dmeas2[1];
02267 er_dmeasR = er_dmeas2[1];
02268 } else {
02269 dmeasR = fabs((*HitMdc).dist()[1]);
02270 er_dmeasR = (*HitMdc).erdist()[1];
02271 }
02272
02273 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
02274 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
02275 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
02276 if(0. == RkR) RkR = 1.e-4;
02277
02278 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
02279
02280 aNewR = v_a + diffR;
02281 double sigR = dmeasR -(v_HT * aNewR)[0];
02282 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
02283 }
02284
02285 // Update Kalman result :
02286 double dchi2_loc(0);
02287
02288 #ifdef YDEBUG
02289 cout<<" track::update_hits......"<<endl;
02290 std::cout << " dchi2R = " << dchi2R << ", dchi2L = " << dchi2L
02291 << " estimate = "<<estim<< std::endl;
02292 std::cout << " dmeasR = " << dmeas2[1]
02293 << ", dmeasL = " << (-1.0)*fabs(dmeas2[0]) << ", check HitMdc.ddl = "
02294 << (*HitMdc).dist()[0] << std::endl;
02295 std::cout<<"dchi2L : "<<dchi2L<<" ,dchi2R: "<<dchi2R<<std::endl;
02296 #endif
02297
02298
02299 if (fabs(dchi2R-dchi2L)<10. && inext>0) {
02300
02301 KalFitHitMdc & HitMdc_next = HitsMdc_[inext];
02302
02303 Helix H_fromR(pivot(), aNewR, eaNewR);
02304 double chi2_fromR(chi2_next(H_fromR, HitMdc_next,csmflag));
02305 //double chi2_fromR(chi2_next(H_fromR, HitMdc_next));
02306
02307 Helix H_fromL(pivot(), aNewL, eaNewL);
02308 double chi2_fromL(chi2_next(H_fromL, HitMdc_next,csmflag));
02309 //double chi2_fromL(chi2_next(H_fromL, HitMdc_next));
02310 #ifdef YDEBUG
02311 std::cout << " chi2_fromL = " << chi2_fromL
02312 << ", chi2_fromR = " << chi2_fromR << std::endl;
02313 #endif
02314 if (fabs(chi2_fromR-chi2_fromL)<10.){
02315 int inext2(-1);
02316 if (inext+1<HitsMdc_.size())
02317 for( int k=inext+1 ; k < HitsMdc_.size(); k++ )
02318 if (!(HitsMdc_[k].chi2()<0)) {
02319 inext2 = k;
02320 break;
02321 }
02322
02323 if (inext2 != -1){
02324 KalFitHitMdc & HitMdc_next2 = HitsMdc_[inext2];
02325 // double chi2_fromR2(chi2_next(H_fromR, HitMdc_next2));
02326 // double chi2_fromL2(chi2_next(H_fromL, HitMdc_next2));
02327 double chi2_fromR2(chi2_next(H_fromR, HitMdc_next2, csmflag));
02328 double chi2_fromL2(chi2_next(H_fromL, HitMdc_next2, csmflag));
02329 #ifdef YDEBUG
02330 std::cout << " chi2_fromL2 = " << chi2_fromL2
02331 << ", chi2_fromR2 = " << chi2_fromR2 << std::endl;
02332 #endif
02333 if (fabs(dchi2R+chi2_fromR+chi2_fromR2-(dchi2L+chi2_fromL+chi2_fromL2))<2) {
02334 if (chi2_fromR2<chi2_fromL2)
02335 dchi2L = DBL_MAX;
02336 else
02337 dchi2R = DBL_MAX;
02338 }
02339 }
02340 }
02341
02342 if (!(dchi2L == DBL_MAX && dchi2R == DBL_MAX)) {
02343 if (dchi2R+chi2_fromR < dchi2L+chi2_fromL){
02344 dchi2L = DBL_MAX;
02345 #ifdef YDEBUG
02346 std::cout << " choose right..." << std::endl;
02347 #endif
02348 } else {
02349 dchi2R = DBL_MAX;
02350 #ifdef YDEBUG
02351 std::cout << " choose left..." << std::endl;
02352 #endif
02353 }
02354 }
02355 }
02356
02357 if ((0 < dchi2R && dchi2R < dchi2cutf_calib[layerid]) ||
02358 (0 < dchi2L && dchi2L < dchi2cutf_calib[layerid])) {
02359
02360 if (((LR_==0 && dchi2R < dchi2L && lr != -1) ||
02361 (LR_==1 && lr == 1)) &&
02362 fabs(aNewR[2]-a()[2]) < 1000. && aNewR[2]) {
02363 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2R < dchi2cutf_calib[layerid]){
02364 nchits_++;
02365 if (flag_ster) nster_++;
02366 if(layerid>19){
02367 Ea(eaNewR);
02368 HelixSeg.Ea_filt_fwd(eaNewR);
02369 a(aNewR);
02370 HelixSeg.a_filt_fwd(aNewR);
02371 }
02372 /*
02373 Ea(eaNewR);
02374 a(aNewR);
02375
02376 KalFitTrack helixR(pivot_work, aNewR, eaNewR, 0, 0, 0);
02377 helixR.pivot(ip);
02378
02379 a_temp = helixR.a();
02380 ea_temp = helixR.Ea();
02381
02382 HelixSeg.Ea_filt_fwd(ea_temp);
02383 HelixSeg.a_filt_fwd(a_temp);
02384 //HelixSeg.filt_include(1);
02385
02386 */
02387
02388 chiSq_ += dchi2R;
02389 dchi2_loc = dchi2R;
02390 if (l_mass_ == lead_){
02391 (*HitMdc).LR(1);
02392 HelixSeg.LR(1);
02393 }
02394 update_bit((*HitMdc).wire().layer().layerId());
02395 }
02396 } else if (((LR_==0 && dchi2L <= dchi2R && lr != 1) ||
02397 (LR_==1 && lr == -1)) &&
02398 fabs(aNewL[2]-a()[2]) < 1000. && aNewL[2]){
02399 if (nchits_ < (unsigned) nmdc_hit2_ || dchi2L < dchi2cutf_calib[layerid]){
02400 nchits_++;
02401 if (flag_ster) nster_++;
02402 if(layerid>19){
02403 Ea(eaNewL);
02404 HelixSeg.Ea_filt_fwd(eaNewL);
02405 a(aNewL);
02406 HelixSeg.a_filt_fwd(aNewL);
02407 }
02408
02409 /*
02410 Ea(eaNewL);
02411 a(aNewL);
02412
02413 KalFitTrack helixL(pivot_work, aNewL, eaNewL, 0, 0, 0);
02414 helixL.pivot(ip);
02415 a_temp = helixL.a();
02416 ea_temp = helixL.Ea();
02417 HelixSeg.Ea_filt_fwd(ea_temp);
02418 HelixSeg.a_filt_fwd(a_temp);
02419 //HelixSeg.filt_include(1);
02420
02421 */
02422
02423 chiSq_ += dchi2L;
02424 dchi2_loc = dchi2L;
02425 if (l_mass_ == lead_){
02426 (*HitMdc).LR(-1);
02427 HelixSeg.LR(-1);
02428 }
02429 update_bit((*HitMdc).wire().layer().layerId());
02430 }
02431 }
02432 }
02433
02434 if (dchi2_loc > dchi2_max_) {
02435 dchi2_max_ = dchi2_loc ;
02436 r_max_ = pivot().perp();
02437 }
02438 dchi2out = dchi2_loc;
02439 // if(dchi2out == 0) {
02440 // dchi2out = ( (dchi2L < dchi2R ) ? dchi2L : dchi2R ) ;
02441 // }
02442 return chiSq_;
02443 }
|
|
|
Record the current parameters as ..._last information :.
|
|
|
Record the current parameters as ..._last information :.
|
|
||||||||||||
|
returns position and convariance matrix(Ex) after rotation.
|
|
||||||||||||
|
|
|
|
returns position after rotating angle dPhi in phi direction.
|
|
||||||||||||
|
returns position and convariance matrix(Ex) after rotation.
00165 {
00166 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
00167 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
00168 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
00169
00170 //
00171 // Calculate position error matrix.
00172 // Ex(phi) = (@x/@a)(Ea)(@x/@a)^T, phi is deflection angle to specify the
00173 // point to be calcualted.
00174 //
00175 // HepMatrix dXDA(3, 5, 0);
00176 // dXDA = delXDelA(phi);
00177 // Ex.assign(dXDA * m_Ea * dXDA.T());
00178
00179 if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
00180 else Ex = m_Ea;
00181
00182 return HepPoint3D(x, y, z);
00183 }
|
|
||||||||||||
|
00148 {
00149 //
00150 // Calculate position (x,y,z) along helix.
00151 //
00152 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
00153 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
00154 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
00155 //
00156
00157 p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
00158 p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
00159 p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
00160
00161 return p;
00162 }
|
|
|
returns position after rotating angle dPhi in phi direction.
00131 {
00132 //
00133 // Calculate position (x,y,z) along helix.
00134 //
00135 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
00136 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
00137 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
00138 //
00139
00140 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
00141 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
00142 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
00143
00144 return HepPoint3D(x, y, z);
00145 }
|
|
|
|
|
|
|
|
|
Flags.
|
|
|
|
|
|
|
|
|
|
|
|
Cut chi2 for each hit.
|
|
|
Cut chi2 for each hit.
|
|
|
|
|
|
|
|
|
|
|
|
Constant alpha for uniform field.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
for debug
|
|
|
the drifttime choice
|
|
|
|
|
|
|
|
|
|
|
|
factor of energy loss straggling for electron
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Flags.
|
|
|
Use L/R decision from MdcRecHit information :.
|
|
|
Initial value: { 0.000510999,
0.105658,
0.139568,
0.493677,
0.938272 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Cut chi2 for each hit.
|
|
|
|
|
|
Flag for treatment of non-uniform mag field.
|
|
|
NUMF treatment improved.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
wire resoltion flag
|
|
|
Level of precision (= 1 : 5steps for all tracks; = 2: 5 steps for very non uniform part) |
|
|
Flag to take account of energy loss straggling :.
|
|
|
|
|
|
|
|
|
Flag for TOF correction.
|
|
|
|
|
|
|
|
|
|
|
|
for signal propagation correction
|
|
|
|
|
|
|
1.3.9.1