CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
KalFitTrack Class Reference

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

#include <KalFitTrack.h>

+ Inheritance diagram for KalFitTrack:

Public Member Functions

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

Static Public Member Functions

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

Static Public Attributes

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

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

Detailed Description

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

Definition at line 35 of file KalFitTrack.h.

Constructor & Destructor Documentation

◆ KalFitTrack()

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

constructor

Definition at line 84 of file KalFitTrack.cxx.

88: KalmanFit::Helix(pivot, a, Ea), type_(0), insist_(0), chiSq_(0),
89 nchits_(0), nster_(0), ncath_(0),
90 ndf_back_(0), chiSq_back_(0), pathip_(0),
91 path_rd_(0), path_ab_(0), tof_(0), dchi2_max_(0), r_max_(0),
92 tof_kaon_(0), tof_proton_(0), pat1_(0), pat2_(0), layer_prec_(0),
93 trasan_id_(0), nhit_r_(0), nhit_z_(0),pathSM_(0),tof2_(0)
94{
95 memset(PathL_, 0, sizeof(PathL_));
96 l_mass_ = m;
97 if (m < NMASS) mass_ = MASS[m];
98 else mass_ = MASS[2];
99 r0_ = fabs(center().perp() - fabs(radius()));
100 //bFieldZ(Bznom_);
101 Bznom_=bFieldZ(); // 2012-09-13 wangll
102 update_last();
103
104 resetLayerUsed();// myLayerUsed
105}
static double Bznom_
void resetLayerUsed()
void update_last(void)
Record the current parameters as ..._last information :
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
const HepSymMatrix & Ea(void) const
returns error matrix.
const HepVector & a(void) const
returns helix parameters.
const HepPoint3D & pivot(void) const
returns pivot position.

◆ ~KalFitTrack()

KalFitTrack::~KalFitTrack ( void )

destructor

Definition at line 108 of file KalFitTrack.cxx.

109{
110 // delete all objects
111
112}

Member Function Documentation

◆ a_forMdc()

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

Definition at line 156 of file KalFitTrack.h.

156{ return a_forMdc_; }

◆ a_last()

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

Definition at line 152 of file KalFitTrack.h.

152{ return a_last_; }

◆ add_nhit_r()

void KalFitTrack::add_nhit_r ( void )
inline

Definition at line 225 of file KalFitTrack.h.

225{ nhit_r_++;}

◆ add_nhit_z()

void KalFitTrack::add_nhit_z ( void )
inline

Definition at line 226 of file KalFitTrack.h.

226{ nhit_z_++;}

◆ addPathSM()

void KalFitTrack::addPathSM ( double path)

Definition at line 1308 of file KalFitTrack.cxx.

1308 {
1309 pathSM_ += path;
1310}

Referenced by KalFitAlg::smoother_anal().

◆ addTofSM()

void KalFitTrack::addTofSM ( double time)

Definition at line 1313 of file KalFitTrack.cxx.

1313 {
1314 tof2_ += time;
1315}
Double_t time

Referenced by KalFitAlg::smoother_anal().

◆ appendHelixSegs()

void KalFitTrack::appendHelixSegs ( KalFitHelixSeg s)

◆ appendHitsMdc()

void KalFitTrack::appendHitsMdc ( KalFitHitMdc h)

Functions for Mdc hits list.

Definition at line 1849 of file KalFitTrack.cxx.

1849{ HitsMdc_.push_back(h);}

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

◆ back() [1/2]

void KalFitTrack::back ( int i)
static

Definition at line 1841 of file KalFitTrack.cxx.

1841{ back_ = i;}

◆ back() [2/2]

int KalFitTrack::back ( void )
static

Definition at line 1842 of file KalFitTrack.cxx.

1842{ return back_;}

Referenced by KalFitAlg::initialize().

◆ chgmass()

void KalFitTrack::chgmass ( int i)

Definition at line 1834 of file KalFitTrack.cxx.

1834 {
1835 mass_=MASS[i];
1836 l_mass_=i;
1837}

Referenced by KalFitAlg::complete_track(), and KalFitAlg::complete_track().

◆ chi2_next() [1/2]

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

Definition at line 2995 of file KalFitTrack.cxx.

2995 {
2996
2997 double lr = HitMdc.LR();
2998 const KalFitWire& Wire = HitMdc.wire();
2999 int wire_ID = Wire.geoID();
3000 int layerid = HitMdc.wire().layer().layerId();
3001 double entrangle = HitMdc.rechitptr()->getEntra();
3002
3003 HepPoint3D fwd(Wire.fwd());
3004 HepPoint3D bck(Wire.bck());
3005 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
3006 KalmanFit::Helix work = H;
3007 work.ignoreErrorMatrix();
3008 work.pivot((fwd + bck) * .5);
3009 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
3010 H.pivot(x0kal);
3011
3012 Hep3Vector meas = H.momentum(0).cross(wire).unit();
3013
3014 if (wire_ID<0 || wire_ID>6796){ //bes
3015 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
3016 << std::endl;
3017 return DBL_MAX;
3018 }
3019
3020 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
3021 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
3022 double tofest(0);
3023 double phi = fmod(phi0() + M_PI4, M_PI2);
3024 double csf0 = cos(phi);
3025 double snf0 = (1. - csf0) * (1. + csf0);
3026 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
3027 if(phi > M_PI) snf0 = - snf0;
3028
3029 if (Tof_correc_) {
3030 Hep3Vector ip(0, 0, 0);
3031 KalmanFit::Helix work = *(KalmanFit::Helix*)this;
3032 work.ignoreErrorMatrix();
3033 work.pivot(ip);
3034 double phi_ip = work.phi0();
3035 if (fabs(phi - phi_ip) > M_PI) {
3036 if (phi > phi_ip) phi -= 2 * M_PI;
3037 else phi_ip -= 2 * M_PI;
3038 }
3039 double t = tanl();
3040 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
3041 double pmag( sqrt( 1.0 + t*t ) / kappa());
3042 double mass_over_p( mass_ / pmag );
3043 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
3044 tofest = l / ( 29.9792458 * beta );
3045 // if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
3046 }
3047
3048 const HepSymMatrix& ea = H.Ea();
3049 const HepVector& v_a = H.a();
3050 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
3051
3052 HepVector v_H(5, 0);
3053 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
3054 v_H[3] = -meas.z();
3055 HepMatrix v_HT = v_H.T();
3056
3057 double estim = (v_HT * v_a)[0];
3058 HepVector ea_v_H = ea * v_H;
3059 HepMatrix ea_v_HT = (ea_v_H).T();
3060 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
3061
3062 HepSymMatrix eaNewL(5), eaNewR(5);
3063 HepVector aNewL(5), aNewR(5);
3064
3065 //double time = HitMdc.tdc();
3066 //if (Tof_correc_)
3067 // time = time - tofest;
3068 double drifttime = getDriftTime(HitMdc , tofest);
3069 double ddl = getDriftDist(HitMdc, drifttime , 0 );
3070 double ddr = getDriftDist(HitMdc, drifttime , 1 );
3071 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
3072 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
3073
3074 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
3075 double er_dmeas2[2] = {0. , 0.};
3076 if(resolflag_ == 1) {
3077 er_dmeas2[0] = 0.1*erddl;
3078 er_dmeas2[1] = 0.1*erddr;
3079 }
3080 else if(resolflag_ == 0) {
3081 // int layid = HitMdc.wire().layer().layerId();
3082 // double sigma = getSigma(layid, dd);
3083 // er_dmeas2[0] = er_dmeas2[1] = sigma;
3084 }
3085
3086 if ((LR_==0 && lr != 1.0) ||
3087 (LR_==1 && lr == -1.0)) {
3088
3089 double er_dmeasL, dmeasL;
3090 if(Tof_correc_) {
3091 dmeasL = (-1.0)*fabs(dmeas2[0]);
3092 er_dmeasL = er_dmeas2[0];
3093 } else {
3094 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
3095 er_dmeasL = HitMdc.erdist()[0];
3096 }
3097
3098 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
3099 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
3100 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
3101 if(0. == RkL) RkL = 1.e-4;
3102
3103 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
3104 aNewL = v_a + diffL;
3105 double sigL = dmeasL -(v_HT * aNewL)[0];
3106 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
3107 }
3108
3109 if ((LR_==0 && lr != -1.0) ||
3110 (LR_==1 && lr == 1.0)) {
3111
3112 double er_dmeasR, dmeasR;
3113 if(Tof_correc_) {
3114 dmeasR = dmeas2[1];
3115 er_dmeasR = er_dmeas2[1];
3116 } else {
3117 dmeasR = fabs(HitMdc.dist()[1]);
3118 er_dmeasR = HitMdc.erdist()[1];
3119 }
3120
3121 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
3122 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
3123 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
3124 if(0. == RkR) RkR = 1.e-4;
3125
3126 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
3127 aNewR = v_a + diffR;
3128 double sigR = dmeasR -(v_HT * aNewR)[0];
3129 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
3130 }
3131
3132 if (dchi2R < dchi2L){
3133 H.a(aNewR); H.Ea(eaNewR);
3134 } else {
3135 H.a(aNewL); H.Ea(eaNewL);
3136 }
3137 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
3138}
double cos(const BesAngle a)
Definition BesAngle.h:213
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
#define DBL_MAX
Definition KalFitTrack.h:4
#define M_PI
Definition TConstant.h:4
int LR(void) const
RecMdcHit * rechitptr(void)
const double * erdist(void) const
const double * dist(void) const
const KalFitWire & wire(void) const
const int layerId(void) const
returns layer ID
double getDriftTime(KalFitHitMdc &hitmdc, double toftime) const
double getSigma(int layerId, double driftDist) const
static int resolflag_
wire resoltion flag
static int Tof_correc_
Flag for TOF correction.
double getDriftDist(KalFitHitMdc &hitmdc, double drifttime, int lr) const
static int LR_
Use L/R decision from MdcRecHit information :
KalFitHitMdc & HitMdc(int i)
Description of a Wire class.
Definition KalFitWire.h:46
unsigned int geoID(void) const
Definition KalFitWire.h:74
HepPoint3D bck(void) const
Definition KalFitWire.h:93
HepPoint3D fwd(void) const
Geometry :
Definition KalFitWire.h:92
const KalFitLayer_Mdc & layer(void) const
Definition KalFitWire.h:70
void ignoreErrorMatrix(void)
unsets error matrix. Error calculations will be ignored after this function call until an error matri...
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
const double getEntra(void) const
Definition RecMdcHit.h:54
IMPLICIT REAL *A H
Definition myXsection.h:1
int t()
Definition t.c:1

◆ chi2_next() [2/2]

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

Definition at line 3140 of file KalFitTrack.cxx.

3140 {
3141
3142 double lr = HitMdc.LR();
3143 const KalFitWire& Wire = HitMdc.wire();
3144 int wire_ID = Wire.geoID();
3145 int layerid = HitMdc.wire().layer().layerId();
3146 double entrangle = HitMdc.rechitptr()->getEntra();
3147
3148 HepPoint3D fwd(Wire.fwd());
3149 HepPoint3D bck(Wire.bck());
3150 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
3151 KalmanFit::Helix work = H;
3152 work.ignoreErrorMatrix();
3153 work.pivot((fwd + bck) * .5);
3154 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
3155 H.pivot(x0kal);
3156
3157 Hep3Vector meas = H.momentum(0).cross(wire).unit();
3158
3159 if (wire_ID<0 || wire_ID>6796){ //bes
3160 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
3161 << std::endl;
3162 return DBL_MAX;
3163 }
3164
3165 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
3166 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
3167 double tofest(0);
3168 double phi = fmod(phi0() + M_PI4, M_PI2);
3169 double csf0 = cos(phi);
3170 double snf0 = (1. - csf0) * (1. + csf0);
3171 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
3172 if(phi > M_PI) snf0 = - snf0;
3173
3174 if (Tof_correc_) {
3175 Hep3Vector ip(0, 0, 0);
3176 KalmanFit::Helix work = *(KalmanFit::Helix*)this;
3177 work.ignoreErrorMatrix();
3178 work.pivot(ip);
3179 double phi_ip = work.phi0();
3180 if (fabs(phi - phi_ip) > M_PI) {
3181 if (phi > phi_ip) phi -= 2 * M_PI;
3182 else phi_ip -= 2 * M_PI;
3183 }
3184 double t = tanl();
3185 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
3186 double pmag( sqrt( 1.0 + t*t ) / kappa());
3187 double mass_over_p( mass_ / pmag );
3188 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
3189 tofest = l / ( 29.9792458 * beta );
3190 if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
3191 }
3192
3193 const HepSymMatrix& ea = H.Ea();
3194 const HepVector& v_a = H.a();
3195 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
3196
3197 HepVector v_H(5, 0);
3198 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
3199 v_H[3] = -meas.z();
3200 HepMatrix v_HT = v_H.T();
3201
3202 double estim = (v_HT * v_a)[0];
3203 HepVector ea_v_H = ea * v_H;
3204 HepMatrix ea_v_HT = (ea_v_H).T();
3205 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
3206
3207 HepSymMatrix eaNewL(5), eaNewR(5);
3208 HepVector aNewL(5), aNewR(5);
3209
3210 //double time = HitMdc.tdc();
3211 //if (Tof_correc_)
3212 // time = time - tofest;
3213 double drifttime = getDriftTime(HitMdc , tofest);
3214 double ddl = getDriftDist(HitMdc, drifttime , 0 );
3215 double ddr = getDriftDist(HitMdc, drifttime , 1 );
3216 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
3217 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
3218
3219 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
3220 double er_dmeas2[2] = {0. , 0.};
3221 if(resolflag_ == 1) {
3222 er_dmeas2[0] = 0.1*erddl;
3223 er_dmeas2[1] = 0.1*erddr;
3224 }
3225 else if(resolflag_ == 0) {
3226 // int layid = HitMdc.wire().layer().layerId();
3227 // double sigma = getSigma(layid, dd);
3228 // er_dmeas2[0] = er_dmeas2[1] = sigma;
3229 }
3230
3231 if ((LR_==0 && lr != 1.0) ||
3232 (LR_==1 && lr == -1.0)) {
3233
3234 double er_dmeasL, dmeasL;
3235 if(Tof_correc_) {
3236 dmeasL = (-1.0)*fabs(dmeas2[0]);
3237 er_dmeasL = er_dmeas2[0];
3238 } else {
3239 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
3240 er_dmeasL = HitMdc.erdist()[0];
3241 }
3242
3243 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
3244 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
3245 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
3246 if(0. == RkL) RkL = 1.e-4;
3247
3248 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
3249 aNewL = v_a + diffL;
3250 double sigL = dmeasL -(v_HT * aNewL)[0];
3251 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
3252 }
3253
3254 if ((LR_==0 && lr != -1.0) ||
3255 (LR_==1 && lr == 1.0)) {
3256
3257 double er_dmeasR, dmeasR;
3258 if(Tof_correc_) {
3259 dmeasR = dmeas2[1];
3260 er_dmeasR = er_dmeas2[1];
3261 } else {
3262 dmeasR = fabs(HitMdc.dist()[1]);
3263 er_dmeasR = HitMdc.erdist()[1];
3264 }
3265
3266 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
3267 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
3268 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
3269 if(0. == RkR) RkR = 1.e-4;
3270
3271 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
3272 aNewR = v_a + diffR;
3273 double sigR = dmeasR -(v_HT * aNewR)[0];
3274 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
3275 }
3276
3277 if (dchi2R < dchi2L){
3278 H.a(aNewR); H.Ea(eaNewR);
3279 } else {
3280 H.a(aNewL); H.Ea(eaNewL);
3281 }
3282 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
3283}
double y(void) const
Definition KalFitWire.h:101

Referenced by update_hits_csmalign().

◆ chiSq() [1/2]

void KalFitTrack::chiSq ( double c)
inline

Definition at line 219 of file KalFitTrack.h.

219{ chiSq_ = c; }

◆ chiSq() [2/2]

double KalFitTrack::chiSq ( void ) const
inline

◆ chiSq_back() [1/2]

void KalFitTrack::chiSq_back ( double c)
inline

Definition at line 220 of file KalFitTrack.h.

220{ chiSq_back_ = c; }

◆ chiSq_back() [2/2]

double KalFitTrack::chiSq_back ( void ) const
inline

Definition at line 189 of file KalFitTrack.h.

189{ return chiSq_back_; }

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

◆ cor_tanldep()

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

Correct the error according the current tanl value :

◆ dchi2_max()

double KalFitTrack::dchi2_max ( void ) const
inline

Definition at line 201 of file KalFitTrack.h.

201{ return dchi2_max_; }

◆ Ea_forMdc()

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

Definition at line 157 of file KalFitTrack.h.

157{ return Ea_forMdc_; }

◆ Ea_last()

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

Definition at line 153 of file KalFitTrack.h.

153{ return Ea_last_; }

◆ eloss()

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

Calculate total energy lost in material.

Definition at line 281 of file KalFitTrack.cxx.

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

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

◆ filter()

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

Definition at line 1266 of file KalFitTrack.cxx.

1268{
1269 HepMatrix m_HT = m_H.T();
1270 HepPoint3D x0 = x(0);
1271 HepVector v_x(3);
1272 v_x[0] = x0.x();
1273 v_x[1] = x0.y();
1274 v_x[2] = x0.z();
1275 HepMatrix m_X = delXDelA(0);
1276 HepMatrix m_XT = m_X.T();
1277 HepMatrix m_C = m_X * Ea() * m_XT;
1278 //int err;
1279 HepVector m_K = 1 / (m_V + (m_HT * m_C * m_H)[0]) * m_H;
1280 HepVector v_a = a();
1281 HepSymMatrix ea = Ea();
1282 HepMatrix mXTmK = m_XT * m_K;
1283 double v_r = v_m - v_d - (m_HT * v_x)[0];
1284 v_a += ea * mXTmK * v_r;
1285 a(v_a);
1286 ea.assign(ea - ea * mXTmK * m_HT * m_X * ea);
1287 Ea(ea);
1288 // Record last hit included parameters :
1289 update_last();
1290 HepMatrix mCmK = m_C * m_K;
1291 v_x += mCmK * v_r;
1292 m_C -= mCmK * m_HT * m_C;
1293 v_r = v_m - v_d - (m_HT * v_x)[0];
1294 double m_R = m_V - (m_HT * m_C * m_H)[0];
1295 double chiSq = v_r * v_r / m_R;
1296 chiSq_ += chiSq;
1297 return chiSq;
1298}
********INTEGER modcns REAL m_C
Definition PseuMar.h:13
double chiSq(void) const

◆ fiTerm()

void KalFitTrack::fiTerm ( double fi)

Definition at line 1318 of file KalFitTrack.cxx.

1318 {
1319 fiTerm_ = fi;
1320}

Referenced by KalFitAlg::smoother_anal().

◆ getClusterRefVec()

ClusterRefVec KalFitTrack::getClusterRefVec ( )
inline

Definition at line 243 of file KalFitTrack.h.

243{ return ClusterRefVec_;} //added hy Huang Zhen

◆ getDigi()

double KalFitTrack::getDigi ( ) const

◆ getDriftDist()

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

Definition at line 194 of file KalFitTrack2.cxx.

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

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

◆ getDriftTime()

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

Definition at line 74 of file KalFitTrack2.cxx.

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

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

◆ getFiTerm()

double KalFitTrack::getFiTerm ( void )
inline

Definition at line 169 of file KalFitTrack.h.

169{ return fiTerm_;}

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

◆ getInitMatrix()

HepSymMatrix KalFitTrack::getInitMatrix ( void ) const

Definition at line 42 of file KalFitTrack2.cxx.

43{
44 return initMatrix_ ;
45}

Referenced by KalFitAlg::smoother_calib().

◆ getPathSM()

double KalFitTrack::getPathSM ( void )
inline

Definition at line 163 of file KalFitTrack.h.

163{ return pathSM_;}

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

◆ getSigma() [1/2]

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

Definition at line 3286 of file KalFitTrack.cxx.

3286 {
3287 double sigma1,sigma2,f;
3288 driftDist *= 10;//mm
3289 if(layerId<8){
3290 if(driftDist<0.5){
3291 sigma1=0.112784; sigma2=0.229274; f=0.666;
3292 }else if(driftDist<1.){
3293 sigma1=0.103123; sigma2=0.269797; f=0.934;
3294 }else if(driftDist<1.5){
3295 sigma1=0.08276; sigma2=0.17493; f=0.89;
3296 }else if(driftDist<2.){
3297 sigma1=0.070109; sigma2=0.149859; f=0.89;
3298 }else if(driftDist<2.5){
3299 sigma1=0.064453; sigma2=0.130149; f=0.886;
3300 }else if(driftDist<3.){
3301 sigma1=0.062383; sigma2=0.138806; f=0.942;
3302 }else if(driftDist<3.5){
3303 sigma1=0.061873; sigma2=0.145696; f=0.946;
3304 }else if(driftDist<4.){
3305 sigma1=0.061236; sigma2=0.119584; f=0.891;
3306 }else if(driftDist<4.5){
3307 sigma1=0.066292; sigma2=0.148426; f=0.917;
3308 }else if(driftDist<5.){
3309 sigma1=0.078074; sigma2=0.188148; f=0.911;
3310 }else if(driftDist<5.5){
3311 sigma1=0.088657; sigma2=0.27548; f=0.838;
3312 }else{
3313 sigma1=0.093089; sigma2=0.115556; f=0.367;
3314 }
3315 }else{
3316 if(driftDist<0.5){
3317 sigma1=0.112433; sigma2=0.327548; f=0.645;
3318 }else if(driftDist<1.){
3319 sigma1=0.096703; sigma2=0.305206; f=0.897;
3320 }else if(driftDist<1.5){
3321 sigma1=0.082518; sigma2=0.248913; f= 0.934;
3322 }else if(driftDist<2.){
3323 sigma1=0.072501; sigma2=0.153868; f= 0.899;
3324 }else if(driftDist<2.5){
3325 sigma1= 0.065535; sigma2=0.14246; f=0.914;
3326 }else if(driftDist<3.){
3327 sigma1=0.060497; sigma2=0.126489; f=0.918;
3328 }else if(driftDist<3.5){
3329 sigma1=0.057643; sigma2= 0.112927; f=0.892;
3330 }else if(driftDist<4.){
3331 sigma1=0.055266; sigma2=0.094833; f=0.887;
3332 }else if(driftDist<4.5){
3333 sigma1=0.056263; sigma2=0.124419; f= 0.932;
3334 }else if(driftDist<5.){
3335 sigma1=0.056599; sigma2=0.124248; f=0.923;
3336 }else if(driftDist<5.5){
3337 sigma1= 0.061377; sigma2=0.146147; f=0.964;
3338 }else if(driftDist<6.){
3339 sigma1=0.063978; sigma2=0.150591; f=0.942;
3340 }else if(driftDist<6.5){
3341 sigma1=0.072951; sigma2=0.15685; f=0.913;
3342 }else if(driftDist<7.){
3343 sigma1=0.085438; sigma2=0.255109; f=0.931;
3344 }else if(driftDist<7.5){
3345 sigma1=0.101635; sigma2=0.315529; f=0.878;
3346 }else{
3347 sigma1=0.149529; sigma2=0.374697; f=0.89;
3348 }
3349 }
3350 double sigmax = sqrt(f*sigma1*sigma1+(1 - f)*sigma2*sigma2)*0.1;
3351 return sigmax;//cm
3352}
double sigma2(0)

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

◆ getSigma() [2/2]

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

Definition at line 209 of file KalFitTrack2.cxx.

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

◆ getT0()

double KalFitTrack::getT0 ( void ) const

Definition at line 57 of file KalFitTrack2.cxx.

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

Referenced by getDriftTime().

◆ getTofSM()

double KalFitTrack::getTofSM ( void )
inline

Definition at line 166 of file KalFitTrack.h.

166{ return tof2_;}

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

◆ HelixSeg()

KalFitHelixSeg & KalFitTrack::HelixSeg ( int i)
inline

Definition at line 240 of file KalFitTrack.h.

240{ return HelixSegs_[i];}

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

◆ HelixSegs() [1/2]

◆ HelixSegs() [2/2]

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

Definition at line 239 of file KalFitTrack.h.

239{ return HelixSegs_;}

◆ HitMdc()

◆ HitsMdc() [1/2]

◆ HitsMdc() [2/2]

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

Definition at line 234 of file KalFitTrack.h.

234{ return HitsMdc_;}

Referenced by order_hits().

◆ insist() [1/2]

void KalFitTrack::insist ( int t)
inline

Definition at line 214 of file KalFitTrack.h.

214{ insist_ = t;}

◆ insist() [2/2]

int KalFitTrack::insist ( void ) const
inline

Extractor :

Definition at line 183 of file KalFitTrack.h.

183{ return insist_; }

◆ intersect_cylinder()

double KalFitTrack::intersect_cylinder ( double r) const

Intersection with different geometry.

Definition at line 128 of file KalFitTrack.cxx.

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

Referenced by KalFitAlg::Cgem_filter_anal(), KalFitAlg::fitGemHits(), KalFitAlg::fromFHitToInnerWall(), KalFitAlg::innerwall(), KalFitCylinder::intersect(), KalFitCylinder::intersect(), KalFitCylinder::intersect(), intersect_yz_plane(), intersect_zx_plane(), and update_hits().

◆ intersect_xy_plane()

double KalFitTrack::intersect_xy_plane ( double z) const

Definition at line 180 of file KalFitTrack.cxx.

181{
182 if (tanl() != 0 && radius() != 0)
183 return (pivot().z() + dz() - z) / (radius() * tanl());
184 else return 0;
185}

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

◆ intersect_yz_plane()

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

Definition at line 162 of file KalFitTrack.cxx.

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

◆ intersect_zx_plane()

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

Definition at line 144 of file KalFitTrack.cxx.

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

◆ lead() [1/2]

void KalFitTrack::lead ( int i)
static

Definition at line 1839 of file KalFitTrack.cxx.

1839{ lead_ = i;}

◆ lead() [2/2]

int KalFitTrack::lead ( void )
static

Magnetic field map.

Definition at line 1840 of file KalFitTrack.cxx.

1840{ return lead_;}

Referenced by KalFitAlg::initialize().

◆ LR()

void KalFitTrack::LR ( int x)
static

Definition at line 1847 of file KalFitTrack.cxx.

1847{ LR_ = x;}

Referenced by KalFitAlg::initialize().

◆ mass() [1/2]

double KalFitTrack::mass ( int i)
static

Definition at line 1833 of file KalFitTrack.cxx.

1833{ return MASS[i];}

◆ mass() [2/2]

double KalFitTrack::mass ( void ) const
inline

◆ mom()

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

Definition at line 195 of file KalFitTrack.h.

195{ return mom_; }

◆ ms()

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

Definition at line 187 of file KalFitTrack.cxx.

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

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

◆ msgasmdc()

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

Calculate multiple scattering angle.

Definition at line 224 of file KalFitTrack.cxx.

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

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

◆ ncath()

unsigned int KalFitTrack::ncath ( void ) const
inline

Definition at line 205 of file KalFitTrack.h.

205{ return ncath_; }

◆ nchits() [1/2]

void KalFitTrack::nchits ( int n)
inline

Definition at line 222 of file KalFitTrack.h.

222{ nchits_ = n; }
const Int_t n

◆ nchits() [2/2]

unsigned int KalFitTrack::nchits ( void ) const
inline

◆ ndf_back() [1/2]

void KalFitTrack::ndf_back ( int n)
inline

Definition at line 221 of file KalFitTrack.h.

221{ ndf_back_ = n; }

◆ ndf_back() [2/2]

int KalFitTrack::ndf_back ( void ) const
inline

Definition at line 190 of file KalFitTrack.h.

190{ return ndf_back_; }

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

◆ nhit_r()

int KalFitTrack::nhit_r ( void ) const
inline

Definition at line 208 of file KalFitTrack.h.

208{ return nhit_r_; }

◆ nhit_z()

int KalFitTrack::nhit_z ( void ) const
inline

Definition at line 209 of file KalFitTrack.h.

209{ return nhit_z_; }

◆ nLayerUsed()

int KalFitTrack::nLayerUsed ( )
inline

Definition at line 343 of file KalFitTrack.h.

344 {
345 int n=0;
346 for(int i=0; i<43; i++) n+=myLayerUsed[i];
347 return n;
348 }

◆ nmass()

int KalFitTrack::nmass ( void )
static

Definition at line 1832 of file KalFitTrack.cxx.

1832{ return NMASS;}

Referenced by KalFitAlg::complete_track(), and KalFitAlg::complete_track().

◆ nster() [1/2]

void KalFitTrack::nster ( int n)
inline

Definition at line 223 of file KalFitTrack.h.

223{ nster_ = n; }

◆ nster() [2/2]

unsigned int KalFitTrack::nster ( void ) const
inline

Definition at line 204 of file KalFitTrack.h.

204{ return nster_; }

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

◆ number_wirhit()

void KalFitTrack::number_wirhit ( void )

Definition at line 464 of file KalFitTrack.cxx.

465{
466 unsigned int nhit = HitsMdc_.size();
467 int Num[50] = {0};
468 for( unsigned i=0 ; i < nhit; i++ )
469 Num[HitsMdc_[i].wire().layer().layerId()]++;
470 for( unsigned i=0 ; i < nhit; i++ )
471 if (Num[HitsMdc_[i].wire().layer().layerId()]>2)
472 HitsMdc_[i].chi2(-2);
473}

◆ numf() [1/2]

void KalFitTrack::numf ( int i)
static

Definition at line 1845 of file KalFitTrack.cxx.

1845{ numf_ = i;}
static int numf_
Flag for treatment of non-uniform mag field.

◆ numf() [2/2]

int KalFitTrack::numf ( void )
static

Definition at line 1846 of file KalFitTrack.cxx.

1846{ return numf_;}

Referenced by KalFitAlg::initialize().

◆ order_hits()

void KalFitTrack::order_hits ( void )

Definition at line 449 of file KalFitTrack.cxx.

449 {
450 for(int it=0; it<HitsMdc().size()-1; it++){
451 if(HitsMdc_[it].wire().layer().layerId() == HitsMdc_[it+1].wire().layer().layerId())
452 {
453 if((kappa()<0)&&(HitsMdc_[it].wire().localId() > HitsMdc_[it+1].wire().localId())){
454 std::swap(HitsMdc_[it], HitsMdc_[it+1]);
455 }
456 if((kappa()>0)&&(HitsMdc_[it].wire().localId() < HitsMdc_[it+1].wire().localId())){
457 std::swap(HitsMdc_[it], HitsMdc_[it+1]);
458 }
459 }
460 }
461}
vector< KalFitHitMdc > & HitsMdc(void)

◆ order_wirhit()

void KalFitTrack::order_wirhit ( int index)

Modifier Order the wire hits for mdc track

Definition at line 380 of file KalFitTrack.cxx.

381{
382 unsigned int nhit = HitsMdc_.size();
383 KalmanFit::Helix tracktest = *(KalmanFit::Helix*)this;
384
385 KalmanFit::Helix helixOrigin = *(KalmanFit::Helix*)this;
386 HepPoint3D origin(0,0,0);
387 helixOrigin.pivot(origin);
388
389 int ind = 0;
390 double* Rad = new double[nhit];
391 double* Ypos = new double[nhit];
392 for( unsigned i=0 ; i < nhit; i++ ){
393
394 HepPoint3D fwd(HitsMdc_[i].wire().fwd());
395 HepPoint3D bck(HitsMdc_[i].wire().bck());
396 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
397
398 // Modification for stereo wires :
399 KalmanFit::Helix work = tracktest;
400 work.ignoreErrorMatrix();
401 work.pivot((fwd + bck) * .5);
402 HepPoint3D x0 = (work.x(0).z() - bck.z())
403 / wire.z() * wire + bck;
404
405 tracktest.pivot(x0);
406 Rad[ind] = tracktest.x(0).perp();// radius !
407 HepPoint3D posOnTrk = tracktest.x(0);
408 //Rad[ind] = helixOrigin.flightArc(posOnTrk);// wangll added 2020-05-22
409 Ypos[ind] = x0.y();
410 ind++;
411 //cout<<"Ypos: "<<Ypos[ind-1]<<endl;
412 }
413
414 // Reorder...
415 if (index < 0)
416 for(int j, k = nhit - 1; k >= 0; k = j){
417 j = -1;
418 for(int i = 1; i <= k; i++)
419 if(Rad[i - 1] > Rad[i]){
420 j = i - 1;
421 std::swap(Rad[i], Rad[j]);
422 std::swap(HitsMdc_[i], HitsMdc_[j]);
423 }
424 }
425 if (index > 0)
426 for(int j, k = nhit - 1; k >= 0; k = j){
427 j = -1;
428 for(int i = 1; i <= k; i++)
429 if(Rad[i - 1] < Rad[i]){
430 j = i - 1;
431 std::swap(Rad[i], Rad[j]);
432 std::swap(HitsMdc_[i], HitsMdc_[j]);
433 }
434 }
435 if (index == 0)
436 for(int j, k = nhit - 1; k >= 0; k = j){
437 j = -1;
438 for(int i = 1; i <= k; i++)
439 if(Ypos[i - 1] > Ypos[i]){
440 j = i - 1;
441 std::swap(Ypos[i], Ypos[j]);
442 std::swap(HitsMdc_[i], HitsMdc_[j]);
443 }
444 }
445 delete [] Rad;
446 delete [] Ypos;
447}
const DifPoint origin

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

◆ p_kaon() [1/2]

void KalFitTrack::p_kaon ( double pl)
inline

Definition at line 217 of file KalFitTrack.h.

217{ p_kaon_ = pl;}

◆ p_kaon() [2/2]

double KalFitTrack::p_kaon ( void ) const
inline

Definition at line 199 of file KalFitTrack.h.

199{ return p_kaon_; }

Referenced by KalFitAlg::complete_track(), and KalFitAlg::complete_track().

◆ p_proton() [1/2]

void KalFitTrack::p_proton ( double pl)
inline

Definition at line 218 of file KalFitTrack.h.

218{ p_proton_ = pl;}

◆ p_proton() [2/2]

double KalFitTrack::p_proton ( void ) const
inline

Definition at line 200 of file KalFitTrack.h.

200{ return p_proton_; }

Referenced by KalFitAlg::complete_track(), and KalFitAlg::complete_track().

◆ pat1()

int KalFitTrack::pat1 ( void ) const
inline

Definition at line 206 of file KalFitTrack.h.

206{ return pat1_; }

◆ pat2()

int KalFitTrack::pat2 ( void ) const
inline

Definition at line 207 of file KalFitTrack.h.

207{ return pat2_; }

◆ path_ab()

double KalFitTrack::path_ab ( void ) const
inline

Definition at line 193 of file KalFitTrack.h.

193{ return path_ab_; }

◆ path_add()

void KalFitTrack::path_add ( double path)

Update the path length estimation.

Definition at line 1301 of file KalFitTrack.cxx.

1302{
1303 pathip_ += path;
1304 tof(path);
1305}

◆ path_rd()

double KalFitTrack::path_rd ( void ) const
inline

Definition at line 192 of file KalFitTrack.h.

192{ return path_rd_; }

◆ pathip() [1/2]

void KalFitTrack::pathip ( double pl)
inline

Definition at line 216 of file KalFitTrack.h.

216{ pathip_ = pl;}

◆ pathip() [2/2]

double KalFitTrack::pathip ( void ) const
inline

Definition at line 191 of file KalFitTrack.h.

191{ return pathip_; }

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

◆ PathL()

double KalFitTrack::PathL ( int layer)

Function to calculate the path length in the layer.

◆ pathl()

double * KalFitTrack::pathl ( void )
inline

Definition at line 194 of file KalFitTrack.h.

194{ return PathL_; }

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

◆ pivot_forMdc()

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

Definition at line 155 of file KalFitTrack.h.

155{ return pivot_forMdc_;}

◆ pivot_last()

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

returns helix parameters

Definition at line 151 of file KalFitTrack.h.

151{ return pivot_last_; }

◆ pivot_numf() [1/2]

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

Sets pivot position in a given mag field.

Definition at line 1386 of file KalFitTrack.cxx.

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

Referenced by KalFitAlg::Cgem_filter_anal(), KalFitAlg::filter_fwd_anal(), KalFitAlg::filter_fwd_calib(), KalFitAlg::fitGemHits(), KalFitAlg::fromFHitToInnerWall(), KalFitAlg::innerwall(), KalFitAlg::smoother_anal(), KalFitAlg::smoother_calib(), update_hits(), KalFitElement::updateTrack(), KalFitElement::updateTrack(), KalFitElement::updateTrack_alreadyfound(), and KalFitElement::updateTrack_rphi().

◆ pivot_numf() [2/2]

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

Definition at line 1490 of file KalFitTrack.cxx.

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

◆ point_last() [1/2]

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

set and give out the last point of the track

Definition at line 146 of file KalFitTrack.h.

146{ point_last_ = point;}

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

◆ point_last() [2/2]

const HepPoint3D & KalFitTrack::point_last ( void )
inline

Definition at line 147 of file KalFitTrack.h.

147{ return point_last_;}

◆ r0()

double KalFitTrack::r0 ( void ) const
inline

Definition at line 186 of file KalFitTrack.h.

186{ return r0_; }

◆ r_max()

double KalFitTrack::r_max ( void ) const
inline

Definition at line 202 of file KalFitTrack.h.

202{ return r_max_; }

◆ radius_numf()

double KalFitTrack::radius_numf ( void ) const

Estimation of the radius in a given mag field.

Definition at line 1348 of file KalFitTrack.cxx.

1348 {
1349
1350 double Bz(Bznom_);
1351 //std::cout<<"Bz: "<<Bz<<std::endl;
1352 if (numf_ > 10){
1353 double dr = a()[0];
1354 double phi0 = a()[1];
1355 double dz = a()[3];
1356 double phi = fmod(phi0 + M_PI4, M_PI2);
1357 double csf0 = cos(phi);
1358 double snf0 = (1. - csf0) * (1. + csf0);
1359 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1360 if(phi > M_PI) snf0 = - snf0;
1361 //XYZPoint
1362 HepPoint3D x0((pivot().x() + dr*csf0),
1363 (pivot().y() + dr*snf0),
1364 (pivot().z() + dz));
1365
1366 //XYZVector b;
1367 HepVector3D b;
1368
1369 //std::cout<<"b: "<<b<<std::endl;
1370
1371 MFSvc_->fieldVector(10.*x0, b);
1372
1373 //std::cout<<"b: "<<b<<std::endl;
1374
1375
1376 b = 10000.*b;
1377 Bz = b.z();
1378 }
1379 if (Bz == 0)
1380 Bz = Bznom_;
1381 double ALPHA_loc = 10000./2.99792458/Bz;
1382 return ALPHA_loc / a()[2];
1383}

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

◆ resetLayerUsed()

void KalFitTrack::resetLayerUsed ( )
inline

Definition at line 350 of file KalFitTrack.h.

350 {
351 for(int i=0; i<43; i++) myLayerUsed[i]=0;
352 }

Referenced by KalFitTrack().

◆ resol() [1/2]

void KalFitTrack::resol ( int i)
static

Definition at line 1843 of file KalFitTrack.cxx.

1843{ resolflag_ = i;}

◆ resol() [2/2]

int KalFitTrack::resol ( void )
static

Definition at line 1844 of file KalFitTrack.cxx.

1844{ return resolflag_;}

Referenced by KalFitAlg::initialize().

◆ setClusterRefVec()

void KalFitTrack::setClusterRefVec ( ClusterRefVec clusterRefVec)
inline

Definition at line 242 of file KalFitTrack.h.

242{ ClusterRefVec_ = clusterRefVec; } //added hy Huang Zhen

◆ setIMdcGeomSvc()

void KalFitTrack::setIMdcGeomSvc ( IMdcGeomSvc * igeomsvc)
static

Definition at line 239 of file KalFitTrack2.cxx.

240{
241 iGeomSvc_ = igeomsvc;
242}

Referenced by KalFitAlg::setGeomSvc_init().

◆ setInitMatrix()

void KalFitTrack::setInitMatrix ( HepSymMatrix m)
static

◆ setMagneticFieldSvc()

void KalFitTrack::setMagneticFieldSvc ( IMagneticFieldSvc * mf)
static

Definition at line 228 of file KalFitTrack2.cxx.

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

Referenced by KalFitAlg::initialize().

◆ setMdcCalibFunSvc()

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

Definition at line 223 of file KalFitTrack2.cxx.

224{
225 CalibFunSvc_ = calibsvc;
226}

Referenced by KalFitAlg::setCalibSvc_init().

◆ setMdcDigiCol()

void KalFitTrack::setMdcDigiCol ( MdcDigiCol * digicol)
static

Definition at line 244 of file KalFitTrack2.cxx.

245{
246 mdcDigiCol_ = digicol;
247}

Referenced by KalFitAlg::execute().

◆ setT0()

void KalFitTrack::setT0 ( double t0)
static

Definition at line 47 of file KalFitTrack2.cxx.

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

Referenced by KalFitAlg::execute().

◆ smoother_Mdc() [1/2]

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

Kalman smoother for Mdc.

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

◆ smoother_Mdc() [2/2]

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

◆ smoother_Mdc_csmalign()

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

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

Definition at line 668 of file KalFitTrack.cxx.

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

Referenced by KalFitAlg::smoother_calib().

◆ tof() [1/2]

void KalFitTrack::tof ( double path)

Update the tof estimation.

Definition at line 1323 of file KalFitTrack.cxx.

1324{
1325 double light_speed( 29.9792458 ); // light speed in cm/nsec
1326 double t = tanl();
1327 double pmag( sqrt( 1.0 + t*t ) / kappa());
1328 if (pmag !=0) {
1329 double mass_over_p( mass_ / pmag );
1330 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
1331 tof_ += path / ( light_speed * beta );
1332 }
1333
1334 if (tofall_) {
1335 if (p_kaon_ > 0){
1336 double massk_over_p( MASS[3] / p_kaon_ );
1337 double beta_kaon( 1.0 / sqrt( 1.0 + massk_over_p * massk_over_p ) );
1338 tof_kaon_ += path / (light_speed * beta_kaon);
1339 }
1340 if (p_proton_ > 0){
1341 double massp_over_p( MASS[4] / p_proton_ );
1342 double beta_proton( 1.0 / sqrt( 1.0 + massp_over_p * massp_over_p ) );
1343 tof_proton_ += path / (light_speed * beta_proton);
1344 }
1345 }
1346}

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

◆ tof() [2/2]

double KalFitTrack::tof ( void ) const
inline

Definition at line 196 of file KalFitTrack.h.

196{ return tof_; }

Referenced by msgasmdc(), and path_add().

◆ tof_kaon()

double KalFitTrack::tof_kaon ( void ) const
inline

Definition at line 197 of file KalFitTrack.h.

197{ return tof_kaon_; }

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

◆ tof_proton()

double KalFitTrack::tof_proton ( void ) const
inline

Definition at line 198 of file KalFitTrack.h.

198{ return tof_proton_; }

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

◆ trasan_id() [1/2]

void KalFitTrack::trasan_id ( int t)
inline

Definition at line 213 of file KalFitTrack.h.

213{ trasan_id_ = t;}

◆ trasan_id() [2/2]

int KalFitTrack::trasan_id ( void ) const
inline

◆ type() [1/2]

void KalFitTrack::type ( int t)
inline

Reinitialize (modificator)

Definition at line 212 of file KalFitTrack.h.

212{ type_ = t;}

◆ type() [2/2]

int KalFitTrack::type ( void ) const
inline

◆ update_bit()

void KalFitTrack::update_bit ( int i)

Definition at line 1819 of file KalFitTrack.cxx.

1819 {
1820 int j(0);
1821 if (i < 31){
1822 j = (int) pow(2.,i);
1823 if (!(pat1_ & j))
1824 pat1_ = pat1_ | j;
1825 } else if (i < 50) {
1826 j = (int) pow(2.,(i-31));
1827 if (!(pat2_ & j))
1828 pat2_ = pat2_ | j;
1829 }
1830}

Referenced by update_hits_csmalign().

◆ update_forMdc()

void KalFitTrack::update_forMdc ( void )

Definition at line 121 of file KalFitTrack.cxx.

122{
123 pivot_forMdc_ = pivot();
124 a_forMdc_ = a();
125 Ea_forMdc_ = Ea();
126}

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

◆ update_hits() [1/3]

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

◆ update_hits() [2/3]

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

◆ update_hits() [3/3]

double KalFitTrack::update_hits ( RecCgemCluster * Cluster,
double recR,
int way,
int csmflag )

Definition at line 2841 of file KalFitTrack.cxx.

2842{
2843 //const double a_stero[3] = {(45.94*3.1415926/180),(31.10*3.1415926/180),(32.99*3.1415926/180)};
2844 //const double r_layer[3] = {87.5026,132.7686,175.2686};
2845 //const double x_reso[3]={0.1372,0.1476,0.1412};
2846 //const double v_reso[3]={0.1273,0.1326,0.1378};
2847
2848 // --- measurement
2849 double recZ = Cluster->getRecZ()/10;// cm
2850 double recPhi = Cluster->getrecphi();
2851 while(recPhi > M_PI) recPhi-=2*M_PI;
2852 while(recPhi <-M_PI) recPhi+=2*M_PI;
2853 HepVector v_measu(2,0);
2854 v_measu(1) = recPhi;
2855 v_measu(2) = recZ;
2856
2857 // --- estimation after prediction, but before filter
2858 //double Phi = intersect_cylinder(recR);
2859 HepPoint3D x0kal = x(intersect_cylinder(recR));
2860 pivot_numf(x0kal);
2861 HepVector v_estim(2,0);
2862 v_estim(1) = x0kal.phi();
2863 v_estim(2) = x0kal.z();
2864
2865 // --- change pivot to (0,0,0)
2866 // HepPoint3D origin(0,0,0); pivot(origin);
2867 // --- get a and Ea
2868 const HepSymMatrix& ea = Ea();
2869 HepVector v_a = a();
2870
2871 // --- derivative matrix
2872 double dPhi=intersect_cylinder(recR);
2873 double dr = v_a(1);
2874 double phi0 = v_a(2);
2875 double kappa = v_a(3);
2876 double tanl = v_a(5);
2877 double x0 = x0kal.x();
2878 double y0 = x0kal.y();
2879 double Alpha = alpha();
2880 HepMatrix H(2, 5, 0);
2881 H(1,1) = -y0*cos(phi0)/(y0*y0+x0*x0)+x0*sin(phi0)/(x0*x0+y0*y0);
2882 H(1,2) = -y0/(y0*y0+x0*x0)*((-1)*dr*sin(phi0)+Alpha/kappa*(sin(phi0+dPhi)-sin(phi0)))+x0/(x0*x0+y0*y0)*(dr*cos(phi0)+Alpha/kappa*(cos(phi0)-cos(phi0+dPhi)));
2883 H(1,3) = -y0/(y0*y0+x0*x0)*(-1)*Alpha/(kappa*kappa)*(cos(phi0)-cos(phi0+dPhi))+x0/(x0*x0+y0*y0)*(-1)*Alpha/(kappa*kappa)*(sin(phi0)-sin(phi0+dPhi));
2884 H(2,3) = Alpha/(kappa*kappa)*tanl*dPhi;
2885 H(2,4) = 1.0;
2886 H(2,5) = -1*(Alpha/kappa)*dPhi;
2887
2888 // --- error matrix of Cgem
2889 int layer = Cluster->getlayerid();
2890 ISvcLocator* svcLocator = Gaudi::svcLocator();
2891 ICgemGeomSvc* ISvc;
2892 StatusCode Cgem_sc=svcLocator->service("CgemGeomSvc", ISvc);
2893 CgemGeomSvc* CgemGeomSvc_=dynamic_cast<CgemGeomSvc *>(ISvc);
2894 CgemGeoLayer* CgemLayer = CgemGeomSvc_->getCgemLayer(layer);
2895 double R_x = (CgemLayer->getInnerROfAnodeCu2() + CgemLayer->getOuterROfAnodeCu2())/2;// in mm
2896 //R_v = (CgemLayer->getInnerROfAnodeCu1() + CgemLayer->getOuterROfAnodeCu1())/2;
2897 double a_stero = (CgemLayer->getAngleOfStereo())*M_PI/180;
2898
2899 ICgemCalibFunSvc* CgemCalibSvc;
2900 StatusCode sc = svcLocator->service ("CgemCalibFunSvc", CgemCalibSvc);
2901 int iView(0), mode(2);
2902 double Q(100), T(100);
2903 double Phi_momentum = momentum(dPhi).phi();
2904 double Phi_position = x0kal.phi();
2905 double delta_phi=Phi_momentum - Phi_position;
2906 while(delta_phi>M_PI) delta_phi-=CLHEP::twopi;
2907 while(delta_phi<-M_PI) delta_phi+=CLHEP::twopi;
2908 double sigma_X = CgemCalibSvc->getSigma(layer,iView,mode,delta_phi,Q,T);// in mm
2909 double sigma_V = CgemCalibSvc->getSigma(layer,iView,mode,delta_phi,Q,T);// in mm
2910
2911 HepSymMatrix V(2,0);
2912 //V(1,1) = pow(x_reso[layer]/r_layer[layer],2);
2913 //V(2,2) = pow(sqrt(0.01*v_reso[layer]*v_reso[layer]/(sin(a_stero[layer])*sin(a_stero[layer]))+0.01*x_reso[layer]*x_reso[layer]/(tan(a_stero[layer])*tan(a_stero[layer]))),2);
2914 V(1,1) = pow(sigma_X/R_x,2);
2915 //V(2,2) = pow(sqrt(0.01*sigma_V*sigma_V/(sin(a_stero)*sin(a_stero))+0.01*sigma_X*sigma_X/(tan(a_stero)*tan(a_stero))),2);
2916 V(2,2) = 0.01*sigma_V*sigma_V/(sin(a_stero)*sin(a_stero))+0.01*sigma_X*sigma_X/(tan(a_stero)*tan(a_stero));
2917 V(1,2) = 0.1*V(1,1)*R_x/tan(a_stero);
2918 V(2,1) = V(1,2);
2919
2920 // --- gain matrix
2921 int ierr=-1;
2922 HepMatrix K = ea*H.T()*(V+H*ea*H.T()).inverse(ierr);
2923 if(ierr != 0) cout<<"KalFitTrack::update_hits(cluster) error in inverse operation for gain matrix!"<<endl;
2924
2925 // --- new error matrix
2926 HepSymMatrix EaNew(5,0);
2927 EaNew.assign(ea-K*H*ea);
2928 //Ea(EaNew);
2929
2930 // --- diff before filter
2931 HepVector v_diff = v_measu - v_estim;
2932 while(v_diff(1) > M_PI) v_diff(1)-=2*M_PI;
2933 while(v_diff(1) <- M_PI) v_diff(1)+=2*M_PI;
2934
2935 // --- new parameters
2936 HepVector aNew = v_a + K*v_diff;
2937
2938 // --- new v_estim
2939 //a(v_a + K*v_diff);// temporary for new v_estim
2940 a(aNew);// temporary for new v_estim
2941 HepPoint3D x0kal_new = x(intersect_cylinder(recR));
2942 HepVector v_estim_new(2,0);
2943 v_estim_new(1) = x0kal_new.phi();
2944 v_estim_new(2) = x0kal_new.z();
2945
2946 // --- difference/residual between measurement and updated estimation
2947 v_diff = v_measu - v_estim_new;
2948 while(v_diff(1) > M_PI) v_diff(1)-=2*M_PI;
2949 while(v_diff(1) <- M_PI) v_diff(1)+=2*M_PI;
2950
2951 // --- new derivative matrix (needed?)
2952 //track_new.pivot_numf(x0kal_new,pathl);
2953 //HepVector a_new = track_new.a();
2954 //HepVector Ea_new = track_new.ea();
2955
2956 // --- R matrix and Chisuqre
2957 HepSymMatrix R(2,0);
2958 R.assign(V-H*EaNew*H.T());
2959 HepVector dChi2 = v_diff.T()*R.inverse(ierr)*v_diff;
2960 if(ierr != 0) cout<<"KalFitTrack::update_hits(cluster) error in inverse operation of R matrix!"<<endl;
2961
2962 // --- check if update the track
2963 int layerid = Cluster->getlayerid();
2964 if(dChi2(1)>0 && fabs(aNew[2]-a()[2]) < 1000. && aNew[2] && dChi2(1) < 2.0*dchi2cutf_anal[layerid])
2965 {
2966 nchits_++;
2967 nchits_++;
2968 nster_++;
2969 // update track parameters and errors
2970 Ea(EaNew);
2971 a(aNew);
2972 // change pivot to x0kal_new
2973 //pivot(x0kal_new);
2974 // --- update Chisqure of track
2975 if(way>0) chiSq(chiSq()+dChi2(1));
2976 if(way<0) chiSq_back(chiSq_back()+dChi2(1));
2977 }
2978 else {
2979 a(v_a);
2980 // change pivot back to x0kal
2981 //pivot(x0kal);
2982 //cout<<"KalTrack drop a cgem cluster, layer="<<layerid
2983 // <<", v_measure="<<v_measu(1)<<", "<<v_measu(2)
2984 // <<", v_estim="<<v_estim(1)<<", "<<v_estim(2)
2985 // <<", v_estim_new="<<v_estim_new(1)<<", "<<v_estim_new(2)
2986 // <<", v_diff="<<v_diff(1)<<", "<<v_diff(2)
2987 // <<", chi2="<<dChi2(1)<<", kappa="<<a()[2]<<", new_kappa="<<aNew[2]<<", chi2_cut="<<2.0*dchi2cutf_anal[layerid]<<endl;
2988 }
2989
2990
2991 return dChi2(1);
2992}
double tan(const BesAngle a)
Definition BesAngle.h:216
double sin(const BesAngle a)
Definition BesAngle.h:210
const double a_stero[3]
double getAngleOfStereo() const
double getOuterROfAnodeCu2() const
double getInnerROfAnodeCu2() const
CgemGeoLayer * getCgemLayer(int i) const
Definition CgemGeomSvc.h:48
virtual double getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const =0
const HepPoint3D & pivot_numf(const HepPoint3D &newPivot)
Sets pivot position in a given mag field.
double chiSq_back(void) const
static double dchi2cutf_anal[43]
Definition KalFitTrack.h:57
double getRecZ(void) const
int getlayerid(void) const
double getrecphi(void) const
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:27

◆ update_hits_csmalign()

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

Definition at line 2140 of file KalFitTrack.cxx.

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

Referenced by KalFitAlg::filter_fwd_calib().

◆ update_last()

void KalFitTrack::update_last ( void )

Record the current parameters as ..._last information :

Definition at line 114 of file KalFitTrack.cxx.

115{
116 pivot_last_ = pivot();
117 a_last_ = a();
118 Ea_last_ = Ea();
119}

Referenced by filter(), and KalFitTrack().

◆ useLayer()

void KalFitTrack::useLayer ( int iLay)
inline

Definition at line 354 of file KalFitTrack.h.

354 {
355 if(iLay>=0 && iLay<=43) myLayerUsed[iLay]=1;
356 }

Member Data Documentation

◆ Bznom_

◆ chi2_hitf_

double KalFitTrack::chi2_hitf_ = 1000
static

Cut chi2 for each hit.

Definition at line 289 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize().

◆ chi2_hits_

double KalFitTrack::chi2_hits_ = 1000
static

Definition at line 289 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize().

◆ chi2mdc_hit2_

double KalFitTrack::chi2mdc_hit2_
static

Definition at line 320 of file KalFitTrack.h.

◆ dchi2cutf_anal

double KalFitTrack::dchi2cutf_anal = {0.}
static

Definition at line 57 of file KalFitTrack.h.

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

◆ dchi2cutf_calib

double KalFitTrack::dchi2cutf_calib = {0.}
static

Definition at line 59 of file KalFitTrack.h.

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

◆ dchi2cuts_anal

double KalFitTrack::dchi2cuts_anal = {0.}
static

Definition at line 58 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut().

◆ dchi2cuts_calib

double KalFitTrack::dchi2cuts_calib = {0.}
static

Definition at line 60 of file KalFitTrack.h.

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

◆ debug_

int KalFitTrack::debug_ = 0
static

◆ drifttime_choice_

int KalFitTrack::drifttime_choice_ = 0
static

the drifttime choice

Definition at line 328 of file KalFitTrack.h.

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

◆ factor_strag_

double KalFitTrack::factor_strag_ = 0.4
static

factor of energy loss straggling for electron

Definition at line 317 of file KalFitTrack.h.

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

◆ inner_steps_

int KalFitTrack::inner_steps_ = 3
static

Definition at line 292 of file KalFitTrack.h.

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

◆ LR_

◆ mdcGasRadlen_

double KalFitTrack::mdcGasRadlen_ = 0.
static

Definition at line 281 of file KalFitTrack.h.

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

◆ nmdc_hit2_

int KalFitTrack::nmdc_hit2_ = 500
static

Cut chi2 for each hit.

Definition at line 319 of file KalFitTrack.h.

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

◆ numf_

int KalFitTrack::numf_ = 0
static

Flag for treatment of non-uniform mag field.

Definition at line 291 of file KalFitTrack.h.

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

◆ numfcor_

int KalFitTrack::numfcor_ = 1
static

NUMF treatment improved.

Definition at line 304 of file KalFitTrack.h.

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

◆ outer_steps_

int KalFitTrack::outer_steps_ = 3
static

Definition at line 293 of file KalFitTrack.h.

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

◆ resolflag_

int KalFitTrack::resolflag_ = 0
static

wire resoltion flag

Definition at line 323 of file KalFitTrack.h.

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

◆ steplev_

int KalFitTrack::steplev_ = 0
static

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

Definition at line 309 of file KalFitTrack.h.

Referenced by KalFitAlg::initialize().

◆ strag_

int KalFitTrack::strag_ = 1
static

Flag to take account of energy loss straggling :

Definition at line 315 of file KalFitTrack.h.

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

◆ Tof_correc_

int KalFitTrack::Tof_correc_ = 1
static

Flag for TOF correction.

Definition at line 312 of file KalFitTrack.h.

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

◆ tofall_

◆ tprop_

int KalFitTrack::tprop_ = 1
static

for signal propagation correction

Definition at line 284 of file KalFitTrack.h.

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


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