LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
genf::GFSpacepointHitPolicy Class Reference

#include "GFSpacepointHitPolicy.h"

Inheritance diagram for genf::GFSpacepointHitPolicy:

Public Member Functions

 GFSpacepointHitPolicy ()
 
const GFDetPlanedetPlane (GFAbsRecoHit *, GFAbsTrackRep *)
 Get detector plane perpendicular to track. More...
 
TMatrixT< double > hitCoord (GFAbsRecoHit *, const GFDetPlane &)
 Hit coordinates in detector plane. More...
 
TMatrixT< double > hitCoord (GFAbsRecoHit *, const GFDetPlane &, const GFDetPlane &)
 
TMatrixT< double > hitCov (GFAbsRecoHit *, const GFDetPlane &)
 Hit covariances in detector plane. More...
 
TMatrixT< double > hitCov (GFAbsRecoHit *hit, const GFDetPlane &plane, const GFDetPlane &planePrev, const TMatrixT< Double_t > &state, const Double_t &mass)
 
virtual ~GFSpacepointHitPolicy ()
 
const std::string & getName ()
 

Private Attributes

GFDetPlane fPlane
 

Static Private Attributes

static const std::string fPolicyName = "GFSpacepointHitPolicy"
 

Detailed Description

Definition at line 55 of file GFSpacepointHitPolicy.h.

Constructor & Destructor Documentation

genf::GFSpacepointHitPolicy::GFSpacepointHitPolicy ( )
inline

Definition at line 59 of file GFSpacepointHitPolicy.h.

References detPlane(), hitCoord(), and hitCov().

59 {;}
virtual genf::GFSpacepointHitPolicy::~GFSpacepointHitPolicy ( )
inlinevirtual

Definition at line 79 of file GFSpacepointHitPolicy.h.

79 {;}

Member Function Documentation

const genf::GFDetPlane & genf::GFSpacepointHitPolicy::detPlane ( GFAbsRecoHit hit,
GFAbsTrackRep rep 
)

Get detector plane perpendicular to track.

The detector plane is contructed from the position of the hit and the track representation. For this the track is extrapolated to the point of closest approach to the hit.

Definition at line 382 of file GFSpacepointHitPolicy.cxx.

References genf::GFAbsTrackRep::extrapolateToPoint(), fPlane, genf::GFAbsRecoHit::getRawHitCoord(), genf::GFDetPlane::setNormal(), and genf::GFDetPlane::setO().

Referenced by GFSpacepointHitPolicy().

383 {
384  TMatrixT<Double_t> rawcoord = hit->getRawHitCoord();
385  TVector3 point(rawcoord[0][0],rawcoord[1][0],rawcoord[2][0]);
386 
387  TVector3 poca,dirInPoca;
388  rep->extrapolateToPoint(point,poca,dirInPoca);
389 
390  fPlane.setO(point);
391  fPlane.setNormal(dirInPoca);
392 
393  return fPlane;
394 }
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:94
Detector simulation of raw signals on wires.
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:158
const std::string& genf::GFSpacepointHitPolicy::getName ( )
inline

Definition at line 81 of file GFSpacepointHitPolicy.h.

References fPolicyName.

81 {return fPolicyName;}
static const std::string fPolicyName
TMatrixT< Double_t > genf::GFSpacepointHitPolicy::hitCoord ( GFAbsRecoHit hit,
const GFDetPlane plane 
)

Hit coordinates in detector plane.

Definition at line 29 of file GFSpacepointHitPolicy.cxx.

References genf::GFDetPlane::getO(), genf::GFAbsRecoHit::getRawHitCoord(), genf::GFDetPlane::getU(), and genf::GFDetPlane::getV().

Referenced by GFSpacepointHitPolicy().

30 {
31  TMatrixT<Double_t> returnMat(2,1);
32 
33  TMatrixT<Double_t> _D(3,1);
34  TVector3 _U;
35  TVector3 _V;
36 
37  _D[0][0] = (plane.getO())[0];
38  _D[1][0] = (plane.getO())[1];
39  _D[2][0] = (plane.getO())[2];
40 
41  _D *= -1.;
42  _D += hit->getRawHitCoord();
43  //now the vector _D points from the origin of the plane to the hit point
44 
45 
46  _U = plane.getU();
47  _V = plane.getV();
48 
49 
50  returnMat[0][0] = _D[0][0] * _U[0] + _D[1][0] * _U[1] + _D[2][0] * _U[2];
51  returnMat[1][0] = _D[0][0] * _V[0] + _D[1][0] * _V[1] + _D[2][0] * _V[2];
52  //std::cout << "hitCoord="<<std::endl;
53  //returnMat.Print();
54  return returnMat;
55 }
Detector simulation of raw signals on wires.
TMatrixT< Double_t > genf::GFSpacepointHitPolicy::hitCoord ( GFAbsRecoHit hit,
const GFDetPlane plane,
const GFDetPlane  
)

Definition at line 59 of file GFSpacepointHitPolicy.cxx.

References genf::GFDetPlane::getO(), genf::GFAbsRecoHit::getRawHitCoord(), genf::GFDetPlane::getU(), and genf::GFDetPlane::getV().

60 {
61  TMatrixT<Double_t> returnMat(5,1); // Just return last 2 elements. Will calculate the rest in GFKalman.cxx.
62 
63  TMatrixT<Double_t> _D(3,1);
64  TVector3 _U;
65  TVector3 _V;
66 
67  _D[0][0] = (plane.getO())[0];
68  _D[1][0] = (plane.getO())[1];
69  _D[2][0] = (plane.getO())[2];
70 
71  _D *= -1.;
72  _D += hit->getRawHitCoord();
73  //now the vector _D points from the origin of the plane to the hit point
74 
75 
76  _U = plane.getU();
77  _V = plane.getV();
78 
79 
80  returnMat[0][0] = _D[0][0] * _U[0] + _D[1][0] * _U[1] + _D[2][0] * _U[2];
81  returnMat[1][0] = _D[0][0] * _V[0] + _D[1][0] * _V[1] + _D[2][0] * _V[2];
82  //std::cout << "hitCoord="<<std::endl;
83  //returnMat.Print();
84  return returnMat;
85 }
Detector simulation of raw signals on wires.
TMatrixT< Double_t > genf::GFSpacepointHitPolicy::hitCov ( GFAbsRecoHit hit,
const GFDetPlane plane 
)

Hit covariances in detector plane.

Definition at line 89 of file GFSpacepointHitPolicy.cxx.

References genf::GFAbsRecoHit::getRawHitCov(), genf::GFDetPlane::getU(), and genf::GFDetPlane::getV().

Referenced by GFSpacepointHitPolicy().

90 {
91  TVector3 _U;
92  TVector3 _V;
93 
94  _U = plane.getU();
95  _V = plane.getV();
96 
97  TMatrixT<Double_t> rawCov = hit->getRawHitCov();
98  //rawCov[0][0] = 12.e12; // Force this. EC, 3-Apr-2012.
99 
100  TMatrixT<Double_t> jac(3,2);
101 
102  // jac = dF_i/dx_j = s_unitvec * t_untivec, with s=u,v and t=x,y,z
103  jac[0][0] = _U[0];
104  jac[1][0] = _U[1];
105  jac[2][0] = _U[2];
106  jac[0][1] = _V[0];
107  jac[1][1] = _V[1];
108  jac[2][1] = _V[2];
109 
110  TMatrixT<Double_t> jac_t = jac.T();
111  TMatrixT<Double_t> jac_orig = jac;
112 
113  TMatrixT<Double_t> result=jac_t * (rawCov * jac_orig);
114  //std::cout << "hitCov="<<std::endl;
115  //result.Print();
116  return result;
117 }
Detector simulation of raw signals on wires.
TMatrixT< Double_t > genf::GFSpacepointHitPolicy::hitCov ( GFAbsRecoHit hit,
const GFDetPlane plane,
const GFDetPlane planePrev,
const TMatrixT< Double_t > &  state,
const Double_t &  mass 
)

Definition at line 120 of file GFSpacepointHitPolicy.cxx.

References beta, genf::GFDetPlane::getNormal(), genf::GFDetPlane::getO(), genf::GFAbsRecoHit::getRawHitCov(), genf::GFDetPlane::getU(), genf::GFDetPlane::getV(), w, X, Y, and Z.

121 {
122  TVector3 _U;
123  TVector3 _V;
124 
125  _U = plane.getU();
126  _V = plane.getV();
127 
128 // Double_t eps(1.0e-6);
129  TMatrixT<Double_t> rawCov = hit->getRawHitCov();
130 
131  //rawCov[0][0] = 12.e12; // Force this. EC, 3-Apr-2012.
132  std::vector <double> tmpRawCov;
133  tmpRawCov.push_back(rawCov[0][0]);
134  tmpRawCov.push_back(rawCov[1][1]);
135  tmpRawCov.push_back(rawCov[2][2]);
136  tmpRawCov.push_back(rawCov[2][1]); // y-z correlated cov element.
137 
138  static std::vector <double> oldRawCov(tmpRawCov);
139  static std::vector <double> oldOldRawCov(tmpRawCov);
140  static GFDetPlane planePrevPrev(planePrev);
141  rawCov.ResizeTo(7,7); // this becomes used now for something else entirely:
142  // the 7x7 raw errors on x,y,z,px,py,pz,th, which sandwiched between 5x7
143  // Jacobian converts it to the cov matrix for the 5x5 state space.
144  rawCov[0][0] = tmpRawCov[0]; // x
145  rawCov[1][1] = tmpRawCov[1]; // y
146  rawCov[2][2] = tmpRawCov[2]; // z
147  rawCov[1][2] = tmpRawCov[3]; // yz
148  rawCov[2][1] = tmpRawCov[3]; // yz
149 
150  // This Jacobian concerns the transfrom from planar to detector coords.
151  // TMatrixT<Double_t> jac(3,2);
152  TMatrixT<Double_t> jac(7,5); // X,Y,Z,UX,UY,UZ,Theta in detector coords
153 
154  // jac = dF_i/dx_j = s_unitvec * t_unitvec, with s=|q|/p,du/dw, dv/dw, u,v and t=th, x,y,z
155 
156  // More correctly :
157  Double_t dist(0.3); // place holder!!
158  Double_t C;
159  //Double_t p (2.5); // place holder!!!
160  //p = abs(1/state[0][0]);
161 
162  Double_t mom = fabs(1.0/state[0][0]);
163  Double_t beta = mom/sqrt(mass*mass+mom*mom);
164  dist = (plane.getO()-planePrev.getO()).Mag();
165  C = 0.0136/beta*sqrt(dist/14.0)*(1+0.038*log(dist/14.0));
166  TVector3 _D (plane.getNormal());
167 
168  // px^ = (x1-x2)/d;
169  // (sigpx)^2 = (sig_x1^2+sig_x2^2)*[(y1-y2)^2 + (z1-z2)^2]/d^4 +
170  // (x1-x2)^2*(y1-y2)^2*sigma_y1^2/d^6 + ...
171  // (x1-x2)^2*(z1-z2)^2*sigma_z2^2/d^6 +
172 
173  rawCov[3][3] = (oldRawCov[0]+tmpRawCov[0])/pow(dist,4.) *
174  (
175  pow((plane.getO()-planePrev.getO()).Y(),2.) +
176  pow((plane.getO()-planePrev.getO()).Z(),2.)
177  )
178  +
179  pow((plane.getO()-planePrev.getO()).X(),2.) *
180  (
181  pow((plane.getO()-planePrev.getO()).Y(),2.) * (oldRawCov[1]+tmpRawCov[1]) +
182  pow((plane.getO()-planePrev.getO()).Z(),2.) * (oldRawCov[2]+tmpRawCov[2])
183  ) / pow(dist,6.)
184  // y-z cross-term
185  +
186  3.0* (plane.getO()-planePrev.getO()).X() * (plane.getO()-planePrev.getO()).Y() * (plane.getO()-planePrev.getO()).Z() / pow(dist,5.) * (tmpRawCov[3] + oldRawCov[3])
187  ;
188 
189  rawCov[4][4] = (oldRawCov[1]+tmpRawCov[1])/pow(dist,4.) *
190  (
191  pow((plane.getO()-planePrev.getO()).X(),2.) +
192  pow((plane.getO()-planePrev.getO()).Z(),2.)
193  )
194  +
195  pow((plane.getO()-planePrev.getO()).Y(),2.) *
196  (
197  pow((plane.getO()-planePrev.getO()).X(),2.) * (oldRawCov[0]+tmpRawCov[0]) +
198  pow((plane.getO()-planePrev.getO()).Z(),2.) * (oldRawCov[2]+tmpRawCov[2])
199  ) / pow(dist,6.)
200  // y-z cross-term
201  +
202  3.0* ( (plane.getO()-planePrev.getO()).X() * (plane.getO()-planePrev.getO()).Y() * (plane.getO()-planePrev.getO()).Z() / pow(dist,5.) + (plane.getO()-planePrev.getO()).Z() / pow(dist,3.) ) * ( tmpRawCov[3] + oldRawCov[3] )
203  ;
204 
205  rawCov[5][5] = (oldRawCov[2]+tmpRawCov[2])/pow(dist,4.) *
206  (
207  pow((plane.getO()-planePrev.getO()).X(),2.) +
208  pow((plane.getO()-planePrev.getO()).Y(),2.)
209  )
210  +
211  pow((plane.getO()-planePrev.getO()).Z(),2.) *
212  (
213  pow((plane.getO()-planePrev.getO()).Y(),2.) * (oldRawCov[0]+tmpRawCov[0])+
214  pow((plane.getO()-planePrev.getO()).X(),2.) * (oldRawCov[1]+tmpRawCov[1])
215  ) / pow(dist,6.)
216  +
217  // y-z cross-term
218  3.0* ( (plane.getO()-planePrev.getO()).X() * (plane.getO()-planePrev.getO()).Y() * (plane.getO()-planePrev.getO()).Z() / pow(dist,5.) + (plane.getO()-planePrev.getO()).Y() / pow(dist,3.) ) * ( tmpRawCov[3] + oldRawCov[3] )
219  ;
220 
221 
222  Double_t num = (plane.getO()-planePrev.getO())*(planePrev.getO()-planePrevPrev.getO());
223  Double_t d1 = (plane.getO()-planePrev.getO()).Mag(); // same as dist above
224  Double_t d2 = (planePrev.getO()-planePrevPrev.getO()).Mag();
225 
226  // This is the error on cosTh^2. There are 9 terms for diagonal errors, one for each
227  // x,y,z coordinate at 3 points. Below the first 3 terms are at pt 1.
228  // Second 3 are at 3. Third 3 are at 2, which is slightly more complicated.
229  // There are three more terms for non-diagonal erros.
230  double dcosTh =
231  pow(((planePrev.getO()-planePrevPrev.getO()).X()*d1*d2 -
232  num * (plane.getO()-planePrev.getO()).X() *d2/d1) /
233  pow(d1,2.0)/pow(d2,2.0),2.0) * tmpRawCov[0] +
234  pow(((planePrev.getO()-planePrevPrev.getO()).Y()*d1*d2 +
235  num * (plane.getO()-planePrev.getO()).Y() *d2/d1) /
236  pow(d1,2.0)/pow(d2,2.0),2.0) * tmpRawCov[1] +
237  pow(((planePrev.getO()-planePrevPrev.getO()).Z()*d1*d2 +
238  num * (plane.getO()-planePrev.getO()).Z() *d2/d1) /
239  pow(d1,2.0)/pow(d2,2.0),2.0) * tmpRawCov[2] +
240 
241  pow(((plane.getO()-planePrev.getO()).X()*d1*d2 -
242  num * (planePrev.getO()-planePrevPrev.getO()).X() *d1/d2) /
243  pow(d1,2.0)/pow(d2,2.0),2.0) * oldOldRawCov[0] +
244  pow(((plane.getO()-planePrev.getO()).Y()*d1*d2 -
245  num * (planePrev.getO()-planePrevPrev.getO()).Y() *d1/d2) /
246  pow(d1,2.0)/pow(d2,2.0),2.0) * oldOldRawCov[1] +
247  pow(((plane.getO()-planePrev.getO()).Z()*d1*d2 -
248  num * (planePrev.getO()-planePrevPrev.getO()).Z() *d1/d2) /
249  pow(d1,2.0)/pow(d2,2.0),2.0) * oldOldRawCov[2] +
250 
251  pow(((plane.getO()-planePrevPrev.getO()).X()*d1*d2 -
252  num * (plane.getO()-planePrev.getO()).X() *d2/d1 +
253  num * (planePrev.getO()-planePrevPrev.getO()).X() *d1/d2
254  ) /
255  pow(d1,2.0)/pow(d2,2.0),2.0) * oldRawCov[0] +
256  pow(((plane.getO()-planePrevPrev.getO()).Y()*d1*d2 -
257  num * (plane.getO()-planePrev.getO()).Y() *d2/d1 +
258  num * (planePrev.getO()-planePrevPrev.getO()).Y() *d1/d2
259  ) /
260  pow(d1,2.0)/pow(d2,2.0),2.0) * oldRawCov[1] +
261  pow(((plane.getO()-planePrevPrev.getO()).Z()*d1*d2 -
262  num * (plane.getO()-planePrev.getO()).Z() *d2/d1 +
263  num * (planePrev.getO()-planePrevPrev.getO()).Z() *d1/d2
264  ) /
265  pow(d1,2.0)/pow(d2,2.0),2.0) * oldRawCov[2]
266  // And now the off-diagonal terms. This is a mess unto itself.
267  // First, d^2(costh)/d(z1)d(y1)
268  +
269  (
270  (plane.getO()-planePrev.getO()).Z() * (planePrev.getO()-planePrevPrev.getO()).Y() / pow(d1,3.)/d2
271  -
272  (plane.getO()-planePrev.getO()).Y() * (planePrev.getO()-planePrevPrev.getO()).Z() / pow(d1,3.)/d2
273  +
274  3.*(plane.getO()-planePrev.getO()).Y() * (plane.getO()-planePrev.getO()).Z() * num / pow(d1,5.)/d2
275  ) * tmpRawCov[3]
276  // Next, d^2(costh)/d(z3)d(y3). Above w z1<->z3, y1<->y3
277  +
278  (
279  (planePrevPrev.getO()-planePrev.getO()).Z() * (planePrev.getO()-plane.getO()).Y() / pow(d1,3.)/d2
280  -
281  (planePrevPrev.getO()-planePrev.getO()).Y() * (planePrev.getO()-plane.getO()).Z() / pow(d1,3.)/d2
282  +
283  3.*(planePrevPrev.getO()-planePrev.getO()).Y() * (planePrevPrev.getO()-planePrev.getO()).Z() * num / pow(d1,5.)/d2
284  ) * oldOldRawCov[3]
285  // Last, d^2(costh)/d(z2)d(y2). This is an even bigger mess.
286  +
287  (
288  ( (plane.getO()-planePrev.getO()).Y() - (planePrev.getO()-planePrevPrev.getO()).Y() ) *
289  (-(plane.getO()-planePrev.getO()).Z()/pow(d1,3.)/d2 +
290  (planePrev.getO()-planePrevPrev.getO()).Z()/pow(d2,3.)/d1
291  )
292  + (plane.getO()-planePrev.getO()).Y() *
293  (
294  (-(planePrev.getO()-planePrevPrev.getO()).Z() +
295  (plane.getO()-planePrev.getO()).Z()
296  ) / pow(d1,3.)/d2
297  - num*
298  ( -3.* (plane.getO()-planePrev.getO()).Z()*d1*d2 +
299  (planePrev.getO()-planePrevPrev.getO()).Z()*pow(d1,3.)/d2
300  )
301  ) / pow(d1,6.) / pow(d2,2.)
302  -
303  (planePrev.getO()-planePrevPrev.getO()).Y()*
304  (
305  (-(planePrev.getO()-planePrevPrev.getO()).Z() +
306  (plane.getO()-planePrev.getO()).Z()
307  ) / pow(d2,3.)/d1
308  - num*
309  ( 3.* (planePrev.getO()-planePrevPrev.getO()).Z()*d1*d2 -
310  (plane.getO()-planePrev.getO()).Z()*pow(d2,3.)/d1
311  )
312  ) / pow(d1,2.) / pow(d2,6.)
313  ) * oldRawCov[3]
314  ;
315 
316  // That was delta(cos(theta))^2. I want delta(theta)^2. dTh^2 = d(cosTh)^2/sinTh^2
317  Double_t theta(TMath::ACos((plane.getO()-planePrev.getO()).Unit() * (planePrev.getO()-planePrevPrev.getO()).Unit()));
318  rawCov[6][6] = TMath::Min(pow(dcosTh,2.)/pow(TMath::Sin(theta),2.),pow(TMath::Pi()/2.0,2.));
319 
320  // This means I'm too close to the endpoints to have histories that allow
321  // proper calculation of above rawCovs. Use below defaults instead.
322  if (d1 == 0 || d2 == 0)
323  {
324  rawCov[3][3] = pow(0.2,2.0); // Unit Px
325  rawCov[4][4] = pow(0.2,2.0); // Unit Py
326  rawCov[5][5] = pow(0.2,2.0); // Unit Pz
327  rawCov[6][6] = pow(0.1,2.0); // theta. 0.3/3mm, say.
328  dist = 0.3;
329  C = 0.0136/beta*sqrt(dist/14.0)*(1+0.038*log(dist/14.0));
330  }
331 
332  // This forces a huge/tiny theta error, which effectively freezes/makes-us-sensitive theta, as is.
333  //rawCov[6][6] = /*9.99e9*/ 0.1;
334 
335  TVector3 u=plane.getU();
336  TVector3 v=plane.getV();
337  TVector3 w=u.Cross(v);
338 
339  // Below line has been done, with code change in GFKalman that has updated
340  // the plane orientation by now.
341  // TVector3 pTilde = 1.0 * (w + state[1][0] * u + state[2][0] * v);
342  TVector3 pTilde = w;
343  double pTildeMag = pTilde.Mag();
344 
345  jac.Zero();
346  jac[6][0] = 1./C; // Should be 1/C?; Had 1 until ... 12-Feb-2013
347 
348  jac[0][3] = _U[0];
349  jac[1][3] = _U[1];
350  jac[2][3] = _U[2];
351 
352  jac[0][4] = _V[0];
353  jac[1][4] = _V[1];
354  jac[2][4] = _V[2];
355 
356  // cnp'd from RKTrackRep.cxx, line ~496
357  // da{x,y,z}/du'
358  jac[3][1] = 1.0/pTildeMag*(u.X()-pTilde.X()/(pTildeMag*pTildeMag)*u*pTilde);
359  jac[4][1] = 1.0/pTildeMag*(u.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*u*pTilde);
360  jac[5][1] = 1.0/pTildeMag*(u.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*u*pTilde);
361  // da{x,y,z}/dv'
362  jac[3][2] = 1.0/pTildeMag*(v.X()-pTilde.X()/(pTildeMag*pTildeMag)*v*pTilde);
363  jac[4][2] = 1.0/pTildeMag*(v.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*v*pTilde);
364  jac[5][2] = 1.0/pTildeMag*(v.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*v*pTilde);
365 
366  TMatrixT<Double_t> jac_orig = jac;
367  TMatrixT<Double_t> jac_t = jac.T();
368 
369  TMatrixT<Double_t> result=jac_t * (rawCov * jac_orig);
370  //std::cout << "hitCov="<<std::endl;
371  //result.Print();
372 
373 
374  oldOldRawCov = oldRawCov;
375  oldRawCov = tmpRawCov;
376  planePrevPrev = planePrev;
377 
378  return result;
379 }
Float_t Y
Definition: plot.C:39
Double_t beta
Float_t Z
Definition: plot.C:39
Detector simulation of raw signals on wires.
Float_t X
Definition: plot.C:39
Float_t w
Definition: plot.C:23

Member Data Documentation

GFDetPlane genf::GFSpacepointHitPolicy::fPlane
private

Definition at line 87 of file GFSpacepointHitPolicy.h.

Referenced by detPlane().

const std::string genf::GFSpacepointHitPolicy::fPolicyName = "GFSpacepointHitPolicy"
staticprivate

Definition at line 84 of file GFSpacepointHitPolicy.h.

Referenced by getName().


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