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