12 const TMatrixT<double>& sigma)
16 for(
int i=0;i<sigma.GetNrows();++i){
17 fCov[i][i]=sigma[i][0];
37 const TMatrixT<double>& _state,
38 const TMatrixT<double>& sigma)
42 for(
int i=0;i<sigma.GetNrows();++i){
43 fCov[i][i]=sigma[i][0];
49 const TMatrixT<double>& sigma,
55 for(
int i=0;i<sigma.GetNrows();++i){
56 fCov[i][i]=sigma[i][0];
67 TMatrixT<double>& statePred){
72 TVector3
w=u.Cross(v);
76 TVector3 wfrom=ufrom.Cross(vfrom);
82 std::cerr<<
"track paralell to detector plane"<<std::endl
83 <<
"extrapolation impossible, aborting"<<std::endl;
87 double t= ((( w * o) - (w * p1)) /
90 double dist=t*dir.Mag();
94 double state0=(p2-o)*u;
95 double state1=(p2-o)*v;
96 double state2=(dir*u)/(dir*w);
97 double state3=(dir*v)/(dir*w);
99 statePred[0].Assign(state0);
100 statePred[1].Assign(state1);
101 statePred[2].Assign(state2);
102 statePred[3].Assign(state3);
106 TMatrixT<double>& statePred,
107 TMatrixT<double>& covPred)
112 TVector3 o=pl.
getO();
113 TVector3 u=pl.
getU();
114 TVector3 v=pl.
getV();
115 TVector3
w=u.Cross(v);
119 TVector3 wfrom=ufrom.Cross(vfrom);
128 TVector3 p1=ofrom+
fState[0][0]*ufrom+
fState[1][0]*vfrom;
139 std::cerr<<
"track paralell to detector plane"<<std::endl
140 <<
"extrapolation impossible, aborting"<<std::endl;
144 double t= ((( w * o) - (w * p1)) /
147 double dist=t*dir.Mag();
149 TVector3 p2=p1+t*
dir;
154 double state0=(p2-o)*u;
155 double state1=(p2-o)*v;
156 double state2=(dir*u)/(dir*w);
157 double state3=(dir*v)/(dir*w);
170 TMatrixT<double> jacobian(4,4);
248 double jacobian0 = (ufrom
249 - 1 / (w *
dir) * (w * ufrom) *
fState[2][0] * ufrom
250 - 1 / (w *
dir) * (w * ufrom) *
fState[3][0] * vfrom
251 - 1 / (w *
dir) * (w * ufrom) * wfrom) * u;
253 double jacobian1 = (vfrom
254 - 1 / (w *
dir) * (w * vfrom) *
fState[2][0] * ufrom
255 - 1 / (w *
dir) * (w * vfrom) *
fState[3][0] * vfrom
256 - 1 / (w *
dir) * (w * vfrom) * wfrom)* u;
258 double jacobian2 = (- pow( (w * dir), -2) * (w * o) *
fState[2][0] * ufrom * (w * ufrom)
259 + 1 / (w *
dir) * (w * o) * ufrom
260 - pow((w * dir), -2) * (w * o) *
fState[3][0] * vfrom * (w * ufrom)
261 - pow((w * dir), -2) * (w * o) * wfrom * (w * ufrom)
262 + pow((w * dir), -2) * (w * p1) *
fState[2][0] * ufrom * (w * ufrom)
263 - 1 / (w *
dir) * (w * p1) * ufrom
264 + pow((w * dir), -2) * (w * p1) *
fState[3][0] * vfrom * (w * ufrom)
265 + pow((w * dir), -2) * (w * p1) * wfrom * (w * ufrom))* u;
267 double jacobian3 = (-pow( ( w * dir), -2) * (w * o) *
fState[2][0] * ufrom * (w * vfrom)
268 - pow((w * dir), -2) * (w * o) *
fState[3][0] * vfrom * (w * vfrom)
269 + 1 / (w *
dir) * (w * o) * vfrom
270 - pow((w * dir), -2) * (w * o) * wfrom * (w * vfrom)
271 + pow((w * dir), -2) * (w * p1) *
fState[2][0] * ufrom * (w * vfrom)
272 + pow((w * dir), -2) * (w * p1) *
fState[3][0] * vfrom * (w * vfrom)
273 - 1 / (w *
dir) * (w * p1) * vfrom
274 + pow((w * dir), -2) * (w * p1) * wfrom * (w * vfrom))* u;
276 double jacobian4 = (ufrom
277 - 1 / (w *
dir) * (w * ufrom) *
fState[2][0] * ufrom
278 - 1 / (w *
dir) * (w * ufrom) *
fState[3][0] * vfrom
279 - 1 / (w *
dir) * (w * ufrom) * wfrom)* v;
281 double jacobian5 = (vfrom
282 - 1 / (w *
dir) * (w * vfrom) *
fState[2][0] * ufrom
283 - 1 / (w *
dir) * (w * vfrom) *
fState[3][0] * vfrom
284 - 1 / (w *
dir) * (w * vfrom) * wfrom)* v;
286 double jacobian6 = (- pow( ( w * dir), -2) * (w * o) *
fState[2][0] * ufrom * (w * ufrom)
287 + 1 / (w *
dir) * (w * o) * ufrom
288 - pow((w * dir), -2) * (w * o) *
fState[3][0] * vfrom * (w * ufrom)
289 - pow((w * dir), -2) * (w * o) * wfrom * (w * ufrom)
290 + pow((w * dir), -2) * (w * p1) *
fState[2][0] * ufrom * (w * ufrom)
291 - 1 / (w *
dir) * (w * p1) * ufrom
292 + pow((w * dir), -2) * (w * p1) *
fState[3][0] * vfrom * (w * ufrom)
293 + pow((w * dir), -2) * (w * p1) * wfrom * (w * ufrom))* v;
295 double jacobian7 = (- pow( ( w * dir), -2) * (w * o) *
fState[2][0] * ufrom * (w * vfrom)
296 - pow((w * dir), -2) * (w * o) *
fState[3][0] * vfrom * (w * vfrom)
297 + 1 / (w *
dir) * (w * o) * vfrom
298 - pow((w * dir), -2) * (w * o) * wfrom * (w * vfrom)
299 + pow((w * dir), -2) * (w * p1) *
fState[2][0] * ufrom * (w * vfrom)
300 + pow((w * dir), -2) * (w * p1) *
fState[3][0] * vfrom * (w * vfrom)
301 - 1 / (w *
dir) * (w * p1) * vfrom
302 + pow((w * dir), -2) * (w * p1) * wfrom * (w * vfrom))* v;
303 double jacobian8 = 0;
304 double jacobian9 = 0;
306 double jacobian10 = ((ufrom * u) / (w * dir)
307 - (dir * u) * pow((w * dir), -2) * (w * ufrom));
309 double jacobian11 = ((vfrom * u) / (w * dir)
310 - (dir * u) * pow((w * dir), -2) * (w * vfrom));
311 double jacobian12 = 0;
312 double jacobian13 = 0;
314 double jacobian14 = ((ufrom * v) / (w * dir)
315 - (dir * v) * pow((w * dir), -2) * (w * ufrom));
317 double jacobian15 = ((vfrom * v) / (w * dir)
318 - (dir * v) * pow((w * dir), -2) * (w * vfrom));
320 jacobian[0][0]=jacobian0; jacobian[0][1]=jacobian1; jacobian[0][2]=jacobian2; jacobian[0][3]=jacobian3;
321 jacobian[1][0]=jacobian4; jacobian[1][1]=jacobian5; jacobian[1][2]=jacobian6; jacobian[1][3]=jacobian7;
322 jacobian[2][0]=jacobian8; jacobian[2][1]=jacobian9; jacobian[2][2]=jacobian10; jacobian[2][3]=jacobian11;
323 jacobian[3][0]=jacobian12; jacobian[3][1]=jacobian13; jacobian[3][2]=jacobian14; jacobian[3][3]=jacobian15;
324 TMatrixT<double> jacobianT(jacobian);
326 covPred=jacobian*
fCov*(jacobianT);
328 statePred[0].Assign(state0);
329 statePred[1].Assign(state1);
330 statePred[2].Assign(state2);
331 statePred[3].Assign(state3);
336 TVector3& dirInPoca){
341 TVector3 wfrom=ufrom.Cross(vfrom);
342 TVector3 pfrom=ofrom +
fState[0][0]*ufrom +
fState[1][0]*vfrom;
347 double t = (dir * pos - dir * pfrom) / (dir * dir);
349 poca=pfrom + t *
dir;
350 dirInPoca=dir.Unit();
354 const TVector3& point2,
357 TVector3& poca_onwire){
362 TVector3 wfrom=ufrom.Cross(vfrom);
364 TVector3 pfrom=ofrom +
fState[0][0]*ufrom +
fState[1][0]*vfrom;
366 TVector3 lineDir=point2-point1;
368 TVector3 normal(dir.y()*lineDir.z()-dir.z()*lineDir.y(),
369 dir.x()*lineDir.z()-dir.z()*lineDir.x(),
370 dir.x()*lineDir.y()-dir.y()*lineDir.x());
372 normal=normal.Unit();
373 TVector3 planeNorm=dir.Cross(normal);
374 double t=(planeNorm*point2-planeNorm*pfrom)/(planeNorm*dir);
376 double t2=(lineDir*poca - lineDir*point1)/(lineDir*lineDir);
377 poca_onwire=point1+lineDir*
t2;
382 TMatrixT<double> statePred(
fState);
386 return pl.
getO()+(statePred[0][0]*pl.
getU())+(statePred[1][0]*pl.
getV());
389 TMatrixT<double> statePred(
fState);
393 TVector3 ret = pl.
getNormal()+(statePred[2][0]*pl.
getU())+(statePred[3][0]*pl.
getV());
virtual void getPosMom(const GFDetPlane &, TVector3 &pos, TVector3 &mom)
unsigned int fDimension
Dimensionality of track representation.
void extrapolateToPoint(const TVector3 &pos, TVector3 &poca, TVector3 &dirInPoca)
This method is to extrapolate the track to point of closest approach to a point in space...
Base Class for genfit track representations. Defines interface for track parameterizations.
TVector3 getNormal() const
TMatrixT< Double_t > fCov
The covariance matrix.
virtual double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
void setReferencePlane(const GFDetPlane &pl)
void extrapolateToLine(const TVector3 &point1, const TVector3 &point2, TVector3 &poca, TVector3 &dirInPoca, TVector3 &poca_onwire)
This method extrapolates to the point of closest approach to a line.
TMatrixT< Double_t > fState
The vector of track parameters.