SlTrackRep.cxx
Go to the documentation of this file.
2 #include <cmath>
3 
5  :GFAbsTrackRep(4),_backw(0)
6 {
7 
8 }
9 
10 genf::SlTrackRep::SlTrackRep(const TMatrixT<double>& _state,
11  const TMatrixT<double>& sigma)
12  :GFAbsTrackRep(4),_backw(0)
13 {
14  fState=_state;
15  for(int i=0;i<sigma.GetNrows();++i){
16  fCov[i][i]=sigma[i][0];
17  }
18  setReferencePlane(GFDetPlane(TVector3(0,0,0),TVector3(1,0,0),TVector3(0,1,0)));
19 }
20 
21 genf::SlTrackRep::SlTrackRep(const TVector3& pos, const TVector3& dir) :GFAbsTrackRep(4),_backw(0)
22 {
23  GFDetPlane d(pos,dir);
25  fState[0][0]=0.;
26  fState[1][0]=0.;
27  fState[2][0]=0.;
28  fState[3][0]=0.;
29  fCov[0][0]=1.;
30  fCov[1][1]=1.;
31  fCov[2][2]=1.;
32  fCov[3][3]=1.;
33 }
34 
36  const TMatrixT<double>& _state,
37  const TMatrixT<double>& sigma)
38  :GFAbsTrackRep(4),_backw(0)
39 {
40  fState=_state;
41  for(int i=0;i<sigma.GetNrows();++i){
42  fCov[i][i]=sigma[i][0];
43  }
45 }
46 
47 genf::SlTrackRep::SlTrackRep(const TMatrixT<double>& _state,
48  const TMatrixT<double>& sigma,
49  const double z)
50  :GFAbsTrackRep(4),_backw(0)
51 {
52 
53  fState=_state;
54  for(int i=0;i<sigma.GetNrows();++i){
55  fCov[i][i]=sigma[i][0];
56  }
57  setReferencePlane(GFDetPlane(TVector3(0,0,z),TVector3(1,0,0),TVector3(0,1,0)));
58 }
59 
60 
62 {
63 
64 }
66  TMatrixT<double>& statePred){
67  statePred.ResizeTo(fDimension,1);
68  TVector3 o=pl.getO();
69  TVector3 u=pl.getU();
70  TVector3 v=pl.getV();
71  TVector3 w=u.Cross(v);
72  TVector3 ofrom=fRefPlane.getO();
73  TVector3 ufrom=fRefPlane.getU();
74  TVector3 vfrom=fRefPlane.getV();
75  TVector3 wfrom=ufrom.Cross(vfrom);
76  TVector3 p1=ofrom+fState[0][0]*ufrom+fState[1][0]*vfrom;
77 
78  TVector3 dir=(fState[2][0] * ufrom + fState[3][0] * vfrom + wfrom);//in global coordinates
79 
80  if(w*dir==0){
81  std::cerr<<"track paralell to detector plane"<<std::endl
82  <<"extrapolation impossible, aborting"<<std::endl;
83  throw;
84  }
85 
86  double t= ((( w * o) - (w * p1)) /
87  (w * dir));
88 
89  double dist=t*dir.Mag();
90 
91  TVector3 p2=p1+t*dir;
92 
93  double state0=(p2-o)*u;
94  double state1=(p2-o)*v;
95  double state2=(dir*u)/(dir*w);
96  double state3=(dir*v)/(dir*w);
97 
98  statePred[0].Assign(state0);
99  statePred[1].Assign(state1);
100  statePred[2].Assign(state2);
101  statePred[3].Assign(state3);
102  return dist;
103 }
105  TMatrixT<double>& statePred,
106  TMatrixT<double>& covPred)
107 {
108  //std::cout<<std::endl<<std::endl<<"extrapolate to plane with state and cov"<<std::endl<<std::endl;
109  statePred.ResizeTo(fDimension,1);
110  covPred.ResizeTo(fDimension,fDimension);
111  TVector3 o=pl.getO();
112  TVector3 u=pl.getU();
113  TVector3 v=pl.getV();
114  TVector3 w=u.Cross(v);
115  TVector3 ofrom=fRefPlane.getO();
116  TVector3 ufrom=fRefPlane.getU();
117  TVector3 vfrom=fRefPlane.getV();
118  TVector3 wfrom=ufrom.Cross(vfrom);
119  /*
120  std::cout<<"ufrom ";
121  ufrom.Print();
122  std::cout<<"vfrom ";
123  vfrom.Print();
124  std::cout<<"wfrom ";
125  wfrom.Print();
126  */
127  TVector3 p1=ofrom+fState[0][0]*ufrom+fState[1][0]*vfrom;
128  /*
129  std::cout<<"p1 ";
130  p1.Print();
131  */
132  TVector3 dir=(fState[2][0] * ufrom + fState[3][0] * vfrom + wfrom);//in global coordinates
133  /*
134  std::cout<<"dir ";
135  dir.Print();
136  */
137  if(w*dir==0){
138  std::cerr<<"track paralell to detector plane"<<std::endl
139  <<"extrapolation impossible, aborting"<<std::endl;
140  throw;
141  }
142 
143  double t= ((( w * o) - (w * p1)) /
144  (w * dir));
145 
146  double dist=t*dir.Mag();
147 
148  TVector3 p2=p1+t*dir;
149  /*
150  std::cout<<"xtrapolate point p2 ";
151  p2.Print();
152  */
153  double state0=(p2-o)*u;
154  double state1=(p2-o)*v;
155  double state2=(dir*u)/(dir*w);
156  double state3=(dir*v)/(dir*w);
157  /*
158  std::cout<<"start state:"<<std::endl;
159  std::cout<<fState[0][0]<<std::endl
160  <<fState[1][0]<<std::endl
161  <<fState[2][0]<<std::endl
162  <<fState[3][0]<<std::endl;
163  std::cout<<"end state:"<<std::endl;
164  std::cout<<fState0<<std::endl
165  <<fState1<<std::endl
166  <<fState2<<std::endl
167  <<fState3<<std::endl<<std::endl;
168  */
169  TMatrixT<double> jacobian(4,4);
170  /*
171  cg3 = .(ufrom
172  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ufrom) * c * ufrom
173  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ufrom) * d * vfrom
174  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ufrom) *wfrom, u);
175 
176  cg4 = .(vfrom
177  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, vfrom) * c * ufrom
178  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, vfrom) * d * vfrom
179  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, vfrom) *wfrom, u);
180 
181  cg5 = .(- pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * c * ufrom * .(w, ufrom)
182  + 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, o) * ufrom
183  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * d * vfrom * .(w, ufrom)
184  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) *wfrom * .(w, ufrom)
185  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * c * ufrom * .(w, ufrom)
186  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ofrom + a * ufrom + b * vfrom) * ufrom
187  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * d * vfrom * .(w, ufrom)
188  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) *wfrom * .(w, ufrom), u);
189 
190  cg6 = .(-pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * c * ufrom * .(w, vfrom)
191  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * d * vfrom * .(w, vfrom)
192  + 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, o) * vfrom
193  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) *wfrom * .(w, vfrom)
194  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * c * ufrom * .(w, vfrom)
195  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * d * vfrom * .(w, vfrom)
196  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ofrom + a * ufrom + b * vfrom) * vfrom
197  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) *wfrom * .(w, vfrom), u);
198 
199  cg7 = .(ufrom
200  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ufrom) * c * ufrom
201  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ufrom) * d * vfrom
202  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ufrom) *wfrom, v);
203 
204  cg8 = .(vfrom
205  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, vfrom) * c * ufrom
206  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, vfrom) * d * vfrom
207  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, vfrom) *wfrom, v);
208 
209  cg9 = .(-pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * c * ufrom * .(w, ufrom)
210  + 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, o) * ufrom
211  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * d * vfrom * .(w, ufrom)
212  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) *wfrom * .(w, ufrom)
213  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * c * ufrom * .(w, ufrom)
214  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ofrom + a * ufrom + b * vfrom) * ufrom
215  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * d * vfrom * .(w, ufrom)
216  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) *wfrom * .(w, ufrom), v);
217 
218  cg10 = .(-pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * c * ufrom * .(w, vfrom)
219  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) * d * vfrom * .(w, vfrom)
220  + 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, o) * vfrom
221  - pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, o) *wfrom * .(w, vfrom)
222  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * c * ufrom * .(w, vfrom)
223  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) * d * vfrom * .(w, vfrom)
224  - 1 / .(w, c * ufrom + d * vfrom +wfrom) * .(w, ofrom + a * ufrom + b * vfrom) * vfrom
225  + pow(.(w, c * ufrom + d * vfrom +wfrom), -2) * .(w, ofrom + a * ufrom + b * vfrom) *wfrom * .(w, vfrom), v);
226 
227  cg11 = 0;
228  cg12 = 0;
229 
230  cg13 = .(ufrom, u) / .(c * ufrom + d * vfrom +wfrom, w)
231  - .(c * ufrom + d * vfrom +wfrom, u) * pow(.(c * ufrom + d * vfrom +wfrom, w), -2) * .(ufrom, w);
232 
233  cg14 = .(vfrom, u) / .(c * ufrom + d * vfrom +wfrom, w)
234  - .(c * ufrom + d * vfrom +wfrom, u) * pow(.(c * ufrom + d * vfrom +wfrom, w), -2) * .(vfrom, w);
235  cg15 = 0;
236  cg16 = 0;
237 
238  cg17 = .(ufrom, v) / .(c * ufrom + d * vfrom +wfrom, w)
239  - .(c * ufrom + d * vfrom +wfrom, v) * pow(.(c * ufrom + d * vfrom +wfrom, w), -2) * .(ufrom, w);
240 
241  cg18 = .(vfrom, v) / .(c * ufrom + d * vfrom +wfrom, w)
242  - .(c * ufrom + d * vfrom +wfrom, v) * pow(.(c * ufrom + d * vfrom +wfrom, w), -2) * .(vfrom, w);
243 
244  */
245 
246 
247  double jacobian0 = (ufrom
248  - 1 / (w * dir) * (w * ufrom) * fState[2][0] * ufrom
249  - 1 / (w * dir) * (w * ufrom) * fState[3][0] * vfrom
250  - 1 / (w * dir) * (w * ufrom) * wfrom) * u;
251 
252  double jacobian1 = (vfrom
253  - 1 / (w * dir) * (w * vfrom) * fState[2][0] * ufrom
254  - 1 / (w * dir) * (w * vfrom) * fState[3][0] * vfrom
255  - 1 / (w * dir) * (w * vfrom) * wfrom)* u;
256 
257  double jacobian2 = (- pow( (w * dir), -2) * (w * o) * fState[2][0] * ufrom * (w * ufrom)
258  + 1 / (w * dir) * (w * o) * ufrom
259  - pow((w * dir), -2) * (w * o) * fState[3][0] * vfrom * (w * ufrom)
260  - pow((w * dir), -2) * (w * o) * wfrom * (w * ufrom)
261  + pow((w * dir), -2) * (w * p1) * fState[2][0] * ufrom * (w * ufrom)
262  - 1 / (w * dir) * (w * p1) * ufrom
263  + pow((w * dir), -2) * (w * p1) * fState[3][0] * vfrom * (w * ufrom)
264  + pow((w * dir), -2) * (w * p1) * wfrom * (w * ufrom))* u;
265 
266  double jacobian3 = (-pow( ( w * dir), -2) * (w * o) * fState[2][0] * ufrom * (w * vfrom)
267  - pow((w * dir), -2) * (w * o) * fState[3][0] * vfrom * (w * vfrom)
268  + 1 / (w * dir) * (w * o) * vfrom
269  - pow((w * dir), -2) * (w * o) * wfrom * (w * vfrom)
270  + pow((w * dir), -2) * (w * p1) * fState[2][0] * ufrom * (w * vfrom)
271  + pow((w * dir), -2) * (w * p1) * fState[3][0] * vfrom * (w * vfrom)
272  - 1 / (w * dir) * (w * p1) * vfrom
273  + pow((w * dir), -2) * (w * p1) * wfrom * (w * vfrom))* u;
274 
275  double jacobian4 = (ufrom
276  - 1 / (w * dir) * (w * ufrom) * fState[2][0] * ufrom
277  - 1 / (w * dir) * (w * ufrom) * fState[3][0] * vfrom
278  - 1 / (w * dir) * (w * ufrom) * wfrom)* v;
279 
280  double jacobian5 = (vfrom
281  - 1 / (w * dir) * (w * vfrom) * fState[2][0] * ufrom
282  - 1 / (w * dir) * (w * vfrom) * fState[3][0] * vfrom
283  - 1 / (w * dir) * (w * vfrom) * wfrom)* v;
284 
285  double jacobian6 = (- pow( ( w * dir), -2) * (w * o) * fState[2][0] * ufrom * (w * ufrom)
286  + 1 / (w * dir) * (w * o) * ufrom
287  - pow((w * dir), -2) * (w * o) * fState[3][0] * vfrom * (w * ufrom)
288  - pow((w * dir), -2) * (w * o) * wfrom * (w * ufrom)
289  + pow((w * dir), -2) * (w * p1) * fState[2][0] * ufrom * (w * ufrom)
290  - 1 / (w * dir) * (w * p1) * ufrom
291  + pow((w * dir), -2) * (w * p1) * fState[3][0] * vfrom * (w * ufrom)
292  + pow((w * dir), -2) * (w * p1) * wfrom * (w * ufrom))* v;
293 
294  double jacobian7 = (- pow( ( w * dir), -2) * (w * o) * fState[2][0] * ufrom * (w * vfrom)
295  - pow((w * dir), -2) * (w * o) * fState[3][0] * vfrom * (w * vfrom)
296  + 1 / (w * dir) * (w * o) * vfrom
297  - pow((w * dir), -2) * (w * o) * wfrom * (w * vfrom)
298  + pow((w * dir), -2) * (w * p1) * fState[2][0] * ufrom * (w * vfrom)
299  + pow((w * dir), -2) * (w * p1) * fState[3][0] * vfrom * (w * vfrom)
300  - 1 / (w * dir) * (w * p1) * vfrom
301  + pow((w * dir), -2) * (w * p1) * wfrom * (w * vfrom))* v;
302  double jacobian8 = 0;
303  double jacobian9 = 0;
304 
305  double jacobian10 = ((ufrom * u) / (w * dir)
306  - (dir * u) * pow((w * dir), -2) * (w * ufrom));
307 
308  double jacobian11 = ((vfrom * u) / (w * dir)
309  - (dir * u) * pow((w * dir), -2) * (w * vfrom));
310  double jacobian12 = 0;
311  double jacobian13 = 0;
312 
313  double jacobian14 = ((ufrom * v) / (w * dir)
314  - (dir * v) * pow((w * dir), -2) * (w * ufrom));
315 
316  double jacobian15 = ((vfrom * v) / (w * dir)
317  - (dir * v) * pow((w * dir), -2) * (w * vfrom));
318 
319  jacobian[0][0]=jacobian0; jacobian[0][1]=jacobian1; jacobian[0][2]=jacobian2; jacobian[0][3]=jacobian3;
320  jacobian[1][0]=jacobian4; jacobian[1][1]=jacobian5; jacobian[1][2]=jacobian6; jacobian[1][3]=jacobian7;
321  jacobian[2][0]=jacobian8; jacobian[2][1]=jacobian9; jacobian[2][2]=jacobian10; jacobian[2][3]=jacobian11;
322  jacobian[3][0]=jacobian12; jacobian[3][1]=jacobian13; jacobian[3][2]=jacobian14; jacobian[3][3]=jacobian15;
323  TMatrixT<double> jacobianT(jacobian);
324  jacobianT.T();
325  covPred=jacobian*fCov*(jacobianT);
326  // std::cout<<"covariance[0][0]"<<covPred[0][0]<<std::endl;
327  statePred[0].Assign(state0);
328  statePred[1].Assign(state1);
329  statePred[2].Assign(state2);
330  statePred[3].Assign(state3);
331  return dist;
332 }
334  TVector3& poca,
335  TVector3& dirInPoca){
336 
337  TVector3 ofrom=fRefPlane.getO();
338  TVector3 ufrom=fRefPlane.getU();
339  TVector3 vfrom=fRefPlane.getV();
340  TVector3 wfrom=ufrom.Cross(vfrom);
341  TVector3 pfrom=ofrom + fState[0][0]*ufrom + fState[1][0]*vfrom;
342  TVector3 dir=(fState[2][0] * ufrom + fState[3][0] * vfrom + wfrom);//in global coordinates
343  // pfrom.Print();
344  dir=dir.Unit();
345  //dir.Print();
346  double t = (dir * pos - dir * pfrom) / (dir * dir);
347 
348  poca=pfrom + t * dir;
349  dirInPoca=dir.Unit();
350 
351 }
352 void genf::SlTrackRep::extrapolateToLine(const TVector3& point1,
353  const TVector3& point2,
354  TVector3& poca,
355  TVector3& dirInPoca,
356  TVector3& poca_onwire){
357 
358  TVector3 ofrom=fRefPlane.getO();
359  TVector3 ufrom=fRefPlane.getU();
360  TVector3 vfrom=fRefPlane.getV();
361  TVector3 wfrom=ufrom.Cross(vfrom);
362 
363  TVector3 pfrom=ofrom + fState[0][0]*ufrom + fState[1][0]*vfrom;
364  TVector3 dir=(fState[2][0] * ufrom + fState[3][0] * vfrom + wfrom);//in global coordinates
365  TVector3 lineDir=point2-point1;
366  //normal normal to both lineDir and Dir
367  TVector3 normal(dir.y()*lineDir.z()-dir.z()*lineDir.y(),
368  dir.x()*lineDir.z()-dir.z()*lineDir.x(),
369  dir.x()*lineDir.y()-dir.y()*lineDir.x());
370 
371  normal=normal.Unit();
372  TVector3 planeNorm=dir.Cross(normal);
373  double t=(planeNorm*point2-planeNorm*pfrom)/(planeNorm*dir);
374  poca=pfrom+t*dir;
375  double t2=(lineDir*poca - lineDir*point1)/(lineDir*lineDir);
376  poca_onwire=point1+lineDir*t2;
377  dirInPoca=dir;
378 
379 }
381  TMatrixT<double> statePred(fState);
382  if(pl!=fRefPlane){
383  extrapolate(pl,statePred);
384  }
385  return pl.getO()+(statePred[0][0]*pl.getU())+(statePred[1][0]*pl.getV());
386 }
388  TMatrixT<double> statePred(fState);
389  if(pl!=fRefPlane){
390  extrapolate(pl,statePred);
391  }
392  TVector3 ret = pl.getNormal()+(statePred[2][0]*pl.getU())+(statePred[3][0]*pl.getV());
393  ret.SetMag(1.);
394  return ret;
395 }
396 void genf::SlTrackRep::getPosMom(const GFDetPlane& pl,TVector3& pos,TVector3& mom){
397  pos=getPos(pl);
398  mom=getMom(pl);
399 }
virtual void getPosMom(const GFDetPlane &, TVector3 &pos, TVector3 &mom)
Definition: SlTrackRep.cxx:396
constexpr T pow(T x)
Definition: pow.h:72
unsigned int fDimension
Dimensionality of track representation.
Definition: GFAbsTrackRep.h:88
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...
Definition: SlTrackRep.cxx:333
string dir
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:83
TVector3 getNormal() const
Definition: GFDetPlane.cxx:145
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:94
TVector3 getO() const
Definition: GFDetPlane.h:81
virtual double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
Definition: SlTrackRep.cxx:104
void setReferencePlane(const GFDetPlane &pl)
Definition: SlTrackRep.h:24
TVector3 getU() const
Definition: GFDetPlane.h:82
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
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.
Definition: SlTrackRep.cxx:352
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
virtual ~SlTrackRep()
Definition: SlTrackRep.cxx:61
TVector3 getV() const
Definition: GFDetPlane.h:83
QTextStream & endl(QTextStream &s)