CGEM BOSS 6.6.5.g
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_
Definition: KalFitTrack.h:306
void resetLayerUsed()
Definition: KalFitTrack.h:350
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().

◆ chi2_next() [1/2]

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

Definition at line 2968 of file KalFitTrack.cxx.

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

3113 {
3114
3115 double lr = HitMdc.LR();
3116 const KalFitWire& Wire = HitMdc.wire();
3117 int wire_ID = Wire.geoID();
3118 int layerid = HitMdc.wire().layer().layerId();
3119 double entrangle = HitMdc.rechitptr()->getEntra();
3120
3121 HepPoint3D fwd(Wire.fwd());
3122 HepPoint3D bck(Wire.bck());
3123 Hep3Vector wire = (Hep3Vector)fwd -(Hep3Vector)bck;
3124 KalmanFit::Helix work = H;
3125 work.ignoreErrorMatrix();
3126 work.pivot((fwd + bck) * .5);
3127 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
3128 H.pivot(x0kal);
3129
3130 Hep3Vector meas = H.momentum(0).cross(wire).unit();
3131
3132 if (wire_ID<0 || wire_ID>6796){ //bes
3133 std::cout << "KalFitTrack : wire_ID problem : " << wire_ID
3134 << std::endl;
3135 return DBL_MAX;
3136 }
3137
3138 double x[3] ={pivot().x(), pivot().y(), pivot().z()};
3139 double pmom[3] ={momentum().x(), momentum().y(), momentum().z()};
3140 double tofest(0);
3141 double phi = fmod(phi0() + M_PI4, M_PI2);
3142 double csf0 = cos(phi);
3143 double snf0 = (1. - csf0) * (1. + csf0);
3144 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
3145 if(phi > M_PI) snf0 = - snf0;
3146
3147 if (Tof_correc_) {
3148 Hep3Vector ip(0, 0, 0);
3149 KalmanFit::Helix work = *(KalmanFit::Helix*)this;
3150 work.ignoreErrorMatrix();
3151 work.pivot(ip);
3152 double phi_ip = work.phi0();
3153 if (fabs(phi - phi_ip) > M_PI) {
3154 if (phi > phi_ip) phi -= 2 * M_PI;
3155 else phi_ip -= 2 * M_PI;
3156 }
3157 double t = tanl();
3158 double l = fabs(radius() * (phi - phi_ip) * sqrt(1 + t * t));
3159 double pmag( sqrt( 1.0 + t*t ) / kappa());
3160 double mass_over_p( mass_ / pmag );
3161 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
3162 tofest = l / ( 29.9792458 * beta );
3163 if(csmflag==1 && HitMdc.wire().y()>0.) tofest= -1. * tofest;
3164 }
3165
3166 const HepSymMatrix& ea = H.Ea();
3167 const HepVector& v_a = H.a();
3168 double dchi2R(DBL_MAX), dchi2L(DBL_MAX);
3169
3170 HepVector v_H(5, 0);
3171 v_H[0] = -csf0 * meas.x() - snf0 * meas.y();
3172 v_H[3] = -meas.z();
3173 HepMatrix v_HT = v_H.T();
3174
3175 double estim = (v_HT * v_a)[0];
3176 HepVector ea_v_H = ea * v_H;
3177 HepMatrix ea_v_HT = (ea_v_H).T();
3178 HepVector v_H_T_ea_v_H = v_HT * ea_v_H;
3179
3180 HepSymMatrix eaNewL(5), eaNewR(5);
3181 HepVector aNewL(5), aNewR(5);
3182
3183 //double time = HitMdc.tdc();
3184 //if (Tof_correc_)
3185 // time = time - tofest;
3186 double drifttime = getDriftTime(HitMdc , tofest);
3187 double ddl = getDriftDist(HitMdc, drifttime , 0 );
3188 double ddr = getDriftDist(HitMdc, drifttime , 1 );
3189 double erddl = getSigma( HitMdc, H.a()[4], 0, ddl);
3190 double erddr = getSigma( HitMdc, H.a()[4], 1, ddr);
3191
3192 double dmeas2[2] = { 0.1*ddl, 0.1*ddr };
3193 double er_dmeas2[2] = {0. , 0.};
3194 if(resolflag_ == 1) {
3195 er_dmeas2[0] = 0.1*erddl;
3196 er_dmeas2[1] = 0.1*erddr;
3197 }
3198 else if(resolflag_ == 0) {
3199 // int layid = HitMdc.wire().layer().layerId();
3200 // double sigma = getSigma(layid, dd);
3201 // er_dmeas2[0] = er_dmeas2[1] = sigma;
3202 }
3203
3204 if ((LR_==0 && lr != 1.0) ||
3205 (LR_==1 && lr == -1.0)) {
3206
3207 double er_dmeasL, dmeasL;
3208 if(Tof_correc_) {
3209 dmeasL = (-1.0)*fabs(dmeas2[0]);
3210 er_dmeasL = er_dmeas2[0];
3211 } else {
3212 dmeasL = (-1.0)*fabs(HitMdc.dist()[0]);
3213 er_dmeasL = HitMdc.erdist()[0];
3214 }
3215
3216 double AL = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasL*er_dmeasL);
3217 eaNewL.assign(ea - ea_v_H * AL * ea_v_HT);
3218 double RkL = 1 - (v_H_T_ea_v_H)[0] * AL;
3219 if(0. == RkL) RkL = 1.e-4;
3220
3221 HepVector diffL = ea_v_H * AL * (dmeasL - estim);
3222 aNewL = v_a + diffL;
3223 double sigL = dmeasL -(v_HT * aNewL)[0];
3224 dchi2L = (sigL * sigL) / (RkL * er_dmeasL*er_dmeasL);
3225 }
3226
3227 if ((LR_==0 && lr != -1.0) ||
3228 (LR_==1 && lr == 1.0)) {
3229
3230 double er_dmeasR, dmeasR;
3231 if(Tof_correc_) {
3232 dmeasR = dmeas2[1];
3233 er_dmeasR = er_dmeas2[1];
3234 } else {
3235 dmeasR = fabs(HitMdc.dist()[1]);
3236 er_dmeasR = HitMdc.erdist()[1];
3237 }
3238
3239 double AR = 1 / ((v_H_T_ea_v_H)[0] + er_dmeasR*er_dmeasR);
3240 eaNewR.assign(ea - ea_v_H * AR * ea_v_HT);
3241 double RkR = 1 - (v_H_T_ea_v_H)[0] * AR;
3242 if(0. == RkR) RkR = 1.e-4;
3243
3244 HepVector diffR = ea_v_H * AR * (dmeasR - estim);
3245 aNewR = v_a + diffR;
3246 double sigR = dmeasR -(v_HT * aNewR)[0];
3247 dchi2R = (sigR*sigR) / (RkR * er_dmeasR*er_dmeasR);
3248 }
3249
3250 if (dchi2R < dchi2L){
3251 H.a(aNewR); H.Ea(eaNewR);
3252 } else {
3253 H.a(aNewL); H.Ea(eaNewL);
3254 }
3255 return ((dchi2R < dchi2L) ? dchi2R : dchi2L);
3256}
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::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_
Definition: KalFitTrack.h:321
static double factor_strag_
factor of energy loss straggling for electron
Definition: KalFitTrack.h:317
static int strag_
Flag to take account of energy loss straggling :
Definition: KalFitTrack.h:315

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_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
Definition: KalFitTrack.h:188

◆ 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
Definition: KalFitTrack.h:287
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(), 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
Definition: KalFitHitMdc.h:32
double getT0(void) const
static int drifttime_choice_
the drifttime choice
Definition: KalFitTrack.h:328
static int tprop_
for signal propagation correction
Definition: KalFitTrack.h:284
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(), 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().

◆ 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().

◆ getSigma() [1/2]

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

Definition at line 3259 of file KalFitTrack.cxx.

3259 {
3260 double sigma1,sigma2,f;
3261 driftDist *= 10;//mm
3262 if(layerId<8){
3263 if(driftDist<0.5){
3264 sigma1=0.112784; sigma2=0.229274; f=0.666;
3265 }else if(driftDist<1.){
3266 sigma1=0.103123; sigma2=0.269797; f=0.934;
3267 }else if(driftDist<1.5){
3268 sigma1=0.08276; sigma2=0.17493; f=0.89;
3269 }else if(driftDist<2.){
3270 sigma1=0.070109; sigma2=0.149859; f=0.89;
3271 }else if(driftDist<2.5){
3272 sigma1=0.064453; sigma2=0.130149; f=0.886;
3273 }else if(driftDist<3.){
3274 sigma1=0.062383; sigma2=0.138806; f=0.942;
3275 }else if(driftDist<3.5){
3276 sigma1=0.061873; sigma2=0.145696; f=0.946;
3277 }else if(driftDist<4.){
3278 sigma1=0.061236; sigma2=0.119584; f=0.891;
3279 }else if(driftDist<4.5){
3280 sigma1=0.066292; sigma2=0.148426; f=0.917;
3281 }else if(driftDist<5.){
3282 sigma1=0.078074; sigma2=0.188148; f=0.911;
3283 }else if(driftDist<5.5){
3284 sigma1=0.088657; sigma2=0.27548; f=0.838;
3285 }else{
3286 sigma1=0.093089; sigma2=0.115556; f=0.367;
3287 }
3288 }else{
3289 if(driftDist<0.5){
3290 sigma1=0.112433; sigma2=0.327548; f=0.645;
3291 }else if(driftDist<1.){
3292 sigma1=0.096703; sigma2=0.305206; f=0.897;
3293 }else if(driftDist<1.5){
3294 sigma1=0.082518; sigma2=0.248913; f= 0.934;
3295 }else if(driftDist<2.){
3296 sigma1=0.072501; sigma2=0.153868; f= 0.899;
3297 }else if(driftDist<2.5){
3298 sigma1= 0.065535; sigma2=0.14246; f=0.914;
3299 }else if(driftDist<3.){
3300 sigma1=0.060497; sigma2=0.126489; f=0.918;
3301 }else if(driftDist<3.5){
3302 sigma1=0.057643; sigma2= 0.112927; f=0.892;
3303 }else if(driftDist<4.){
3304 sigma1=0.055266; sigma2=0.094833; f=0.887;
3305 }else if(driftDist<4.5){
3306 sigma1=0.056263; sigma2=0.124419; f= 0.932;
3307 }else if(driftDist<5.){
3308 sigma1=0.056599; sigma2=0.124248; f=0.923;
3309 }else if(driftDist<5.5){
3310 sigma1= 0.061377; sigma2=0.146147; f=0.964;
3311 }else if(driftDist<6.){
3312 sigma1=0.063978; sigma2=0.150591; f=0.942;
3313 }else if(driftDist<6.5){
3314 sigma1=0.072951; sigma2=0.15685; f=0.913;
3315 }else if(driftDist<7.){
3316 sigma1=0.085438; sigma2=0.255109; f=0.931;
3317 }else if(driftDist<7.5){
3318 sigma1=0.101635; sigma2=0.315529; f=0.878;
3319 }else{
3320 sigma1=0.149529; sigma2=0.374697; f=0.89;
3321 }
3322 }
3323 double sigmax = sqrt(f*sigma1*sigma1+(1 - f)*sigma2*sigma2)*0.1;
3324 return sigmax;//cm
3325}
double sigma2(0)

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

◆ getSigma() [2/2]

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

Definition at line 209 of file KalFitTrack2.cxx.

210{
211 int layerid = hitmdc.wire().layer().layerId();
212 double entrangle = hitmdc.rechitptr()->getEntra();
213 // double tanlam = hitmdc.rechitptr()->getTanl();
214 double z = hitmdc.rechitptr()->getZhit();
215 double Q = hitmdc.rechitptr()->getAdc();
216 //std::cout<<" the driftdist before getsigma is "<<dist<<" the layer is"<<layerid<<std::endl;
217 //cout<<"layerid, lr, dist, entrangle, tanlam, z , Q = "<<layerid<<", "<<lr<<", "<<dist<<", "<<entrangle<<", "<<tanlam<<", "<<z<<", "<<Q<<endl;//wangll
218 double temp = CalibFunSvc_->getSigma(layerid, lr, dist, entrangle, tanlam, z , Q );
219 //std::cout<<" the sigma is "<<temp<<std::endl;
220 return temp;
221}
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().

◆ 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()

KalFitHitMdc & KalFitTrack::HitMdc ( int  i)
inline

◆ 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(), 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().

◆ 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

Definition at line 187 of file KalFitTrack.h.

187{ return mass_; }

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

◆ mom()

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

Definition at line 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_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
Definition: KalFitTrack.h:196
static double mdcGasRadlen_
Definition: KalFitTrack.h:281

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(), 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().

◆ 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.
Definition: KalFitTrack.h:291

◆ 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)
Definition: KalFitTrack.h:234

◆ 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().

◆ 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().

◆ 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().

◆ 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(), 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_
Definition: KalFitTrack.h:292
static int numfcor_
NUMF treatment improved.
Definition: KalFitTrack.h:304
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_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)
Definition: KalFitTrack.h:194
static int outer_steps_
Definition: KalFitTrack.h:293

◆ 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(), 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(), 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.
Definition: KalFitHitMdc.h:17
Description of a track class (<- Helix.cc)
Definition: KalFitTrack.h:36
static double dchi2cuts_calib[43]
Definition: KalFitTrack.h:301

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().

◆ 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().

◆ 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().

◆ 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;
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
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 const HepSymMatrix& ea = Ea();
2866 HepVector v_a = a();
2867
2868 // --- derivative matrix
2869 double dPhi=intersect_cylinder(recR);
2870 double dr = v_a(1);
2871 double phi0 = v_a(2);
2872 double kappa = v_a(3);
2873 double tanl = v_a(5);
2874 double x0 = x0kal.x();
2875 double y0 = x0kal.y();
2876 double Alpha = alpha();
2877 HepMatrix H(2, 5, 0);
2878 H(1,1) = -y0*cos(phi0)/(y0*y0+x0*x0)+x0*sin(phi0)/(x0*x0+y0*y0);
2879 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)));
2880 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));
2881 H(2,3) = Alpha/(kappa*kappa)*tanl*dPhi;
2882 H(2,4) = 1.0;
2883 H(2,5) = -1*(Alpha/kappa)*dPhi;
2884
2885 // --- error matrix of Cgem
2886 int layer = Cluster->getlayerid();
2887 ISvcLocator* svcLocator = Gaudi::svcLocator();
2888 ICgemGeomSvc* ISvc;
2889 StatusCode Cgem_sc=svcLocator->service("CgemGeomSvc", ISvc);
2890 CgemGeomSvc* CgemGeomSvc_=dynamic_cast<CgemGeomSvc *>(ISvc);
2891 CgemGeoLayer* CgemLayer = CgemGeomSvc_->getCgemLayer(layer);
2892 double R_x = (CgemLayer->getInnerROfAnodeCu2() + CgemLayer->getOuterROfAnodeCu2())/2;
2893 //R_v = (CgemLayer->getInnerROfAnodeCu1() + CgemLayer->getOuterROfAnodeCu1())/2;
2894 double a_stero = (CgemLayer->getAngleOfStereo())*M_PI/180;
2895
2896 ICgemCalibFunSvc* CgemCalibSvc;
2897 StatusCode sc = svcLocator->service ("CgemCalibFunSvc", CgemCalibSvc);
2898 int iView(0), mode(2);
2899 double Q(100), T(100);
2900 double Phi_momentum = momentum(dPhi).phi();
2901 double Phi_position = x0kal.phi();
2902 double delta_phi=Phi_momentum - Phi_position;
2903 while(delta_phi>M_PI) delta_phi-=CLHEP::twopi;
2904 while(delta_phi<-M_PI) delta_phi+=CLHEP::twopi;
2905 double sigma_X = CgemCalibSvc->getSigma(layer,iView,mode,delta_phi,Q,T);// in mm
2906 double sigma_V = CgemCalibSvc->getSigma(layer,iView,mode,delta_phi,Q,T);// in mm
2907
2908 HepSymMatrix V(2,0);
2909 //V(1,1) = pow(x_reso[layer]/r_layer[layer],2);
2910 //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);
2911 V(1,1) = pow(sigma_X/R_x,2);
2912 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);
2913 //V(1,2) = V(1,1)*R_x/tan(a_stero);
2914 //V(2,1) = V(1,2);
2915
2916 // --- gain matrix
2917 int ierr=-1;
2918 HepMatrix K = ea*H.T()*(V+H*ea*H.T()).inverse(ierr);
2919 if(ierr != 0) cout<<"KalFitTrack::update_hits(cluster) errer in inverse operation of gain matrix!"<<endl;
2920
2921 // --- new error matrix
2922 HepSymMatrix EaNew(5,0);
2923 EaNew.assign(ea-K*H*ea);
2924 Ea(EaNew);
2925
2926 // --- diff
2927 HepVector v_diff = v_measu - v_estim;
2928 while(v_diff(1) > M_PI) v_diff(1)-=2*M_PI;
2929 while(v_diff(1) <- M_PI) v_diff(1)+=2*M_PI;
2930
2931 // --- new parameters
2932 //HepVector aNew = v_a + K*v_diff;
2933 a(v_a + K*v_diff);
2934
2935 nchits_++;
2936 nchits_++;
2937 nster_++;
2938
2939 // --- new v_estim
2940 HepPoint3D x0kal_new = x(intersect_cylinder(recR));
2941 HepVector v_estim_new(2,0);
2942 v_estim_new(1) = x0kal_new.phi();
2943 v_estim_new(2) = x0kal_new.z();
2944
2945 // --- difference between measurement and updated estimation
2946 v_diff = v_measu - v_estim_new;
2947 while(v_diff(1) > M_PI) v_diff(1)-=2*M_PI;
2948 while(v_diff(1) <- M_PI) v_diff(1)+=2*M_PI;
2949
2950 // --- new derivative matrix
2951 //track_new.pivot_numf(x0kal_new,pathl);
2952 //HepVector a_new = track_new.a();
2953 //HepVector Ea_new = track_new.ea();
2954
2955 // --- R matrix and Chisuqre
2956 HepSymMatrix R(2,0);
2957 R.assign(V-H*EaNew*H.T());
2958 HepVector dChi2 = v_diff.T()*R.inverse(ierr)*v_diff;
2959 if(ierr != 0)cout<<"KalFitTrack::update_hits(cluster) errer in inverse operation of R matrix!"<<endl;
2960
2961 // --- update Chisqure of track
2962 if(way>0) chiSq(chiSq()+dChi2(1));
2963 if(way<0) chiSq_back(chiSq_back()+dChi2(1));
2964 return dChi2(1);
2965}
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
Definition: CgemGeoLayer.h:218
double getOuterROfAnodeCu2() const
Definition: CgemGeoLayer.h:168
double getInnerROfAnodeCu2() const
Definition: CgemGeoLayer.h:167
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
Definition: KalFitTrack.h:189
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.
Definition: KalFitTrack.h:319
KalFitHelixSeg & HelixSeg(int i)
Definition: KalFitTrack.h:240
double chi2_next(KalmanFit::Helix &H, KalFitHitMdc &HitMdc, int csmflag)
static double dchi2cutf_calib[43]
Definition: KalFitTrack.h:299
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 295 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut().

◆ dchi2cutf_calib

double KalFitTrack::dchi2cutf_calib = {0.}
static

Definition at line 299 of file KalFitTrack.h.

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

◆ dchi2cuts_anal

double KalFitTrack::dchi2cuts_anal = {0.}
static

Definition at line 297 of file KalFitTrack.h.

Referenced by KalFitAlg::setDchisqCut().

◆ dchi2cuts_calib

double KalFitTrack::dchi2cuts_calib = {0.}
static

Definition at line 301 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(), 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(), 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(), KalFitAlg::initialize(), smoother_Mdc_csmalign(), and update_hits_csmalign().

◆ tofall_

int KalFitTrack::tofall_ = 1
static

◆ 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: