LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GFSpacepointHitPolicy.cxx
Go to the documentation of this file.
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 #include "larreco/Genfit/GFSpacepointHitPolicy.h"
20 #include "larreco/Genfit/GFAbsTrackRep.h"
21 
22 #include "TMath.h"
23 
24 #include "larreco/Genfit/GFAbsRecoHit.h"
25 
26 const std::string genf::GFSpacepointHitPolicy::fPolicyName = "GFSpacepointHitPolicy";
27 
28 TMatrixT<Double_t> genf::GFSpacepointHitPolicy::hitCoord(GFAbsRecoHit* hit, const GFDetPlane& plane)
29 {
30  TMatrixT<Double_t> returnMat(2, 1);
31 
32  TMatrixT<Double_t> _D(3, 1);
33  TVector3 _U;
34  TVector3 _V;
35 
36  _D[0][0] = (plane.getO())[0];
37  _D[1][0] = (plane.getO())[1];
38  _D[2][0] = (plane.getO())[2];
39 
40  _D *= -1.;
41  _D += hit->getRawHitCoord();
42  //now the vector _D points from the origin of the plane to the hit point
43 
44  _U = plane.getU();
45  _V = plane.getV();
46 
47  returnMat[0][0] = _D[0][0] * _U[0] + _D[1][0] * _U[1] + _D[2][0] * _U[2];
48  returnMat[1][0] = _D[0][0] * _V[0] + _D[1][0] * _V[1] + _D[2][0] * _V[2];
49  //std::cout << "hitCoord="<<std::endl;
50  //returnMat.Print();
51  return returnMat;
52 }
53 
54 TMatrixT<Double_t> genf::GFSpacepointHitPolicy::hitCoord(GFAbsRecoHit* hit,
55  const GFDetPlane& plane,
56  const GFDetPlane& /* planePrev */)
57 {
58  TMatrixT<Double_t> returnMat(
59  5, 1); // Just return last 2 elements. Will calculate the rest in GFKalman.cxx.
60 
61  TMatrixT<Double_t> _D(3, 1);
62  TVector3 _U;
63  TVector3 _V;
64 
65  _D[0][0] = (plane.getO())[0];
66  _D[1][0] = (plane.getO())[1];
67  _D[2][0] = (plane.getO())[2];
68 
69  _D *= -1.;
70  _D += hit->getRawHitCoord();
71  //now the vector _D points from the origin of the plane to the hit point
72 
73  _U = plane.getU();
74  _V = plane.getV();
75 
76  returnMat[0][0] = _D[0][0] * _U[0] + _D[1][0] * _U[1] + _D[2][0] * _U[2];
77  returnMat[1][0] = _D[0][0] * _V[0] + _D[1][0] * _V[1] + _D[2][0] * _V[2];
78  //std::cout << "hitCoord="<<std::endl;
79  //returnMat.Print();
80  return returnMat;
81 }
82 
83 TMatrixT<Double_t> genf::GFSpacepointHitPolicy::hitCov(GFAbsRecoHit* hit, const GFDetPlane& plane)
84 {
85  TVector3 _U;
86  TVector3 _V;
87 
88  _U = plane.getU();
89  _V = plane.getV();
90 
91  TMatrixT<Double_t> rawCov = hit->getRawHitCov();
92  //rawCov[0][0] = 12.e12; // Force this. EC, 3-Apr-2012.
93 
94  TMatrixT<Double_t> jac(3, 2);
95 
96  // jac = dF_i/dx_j = s_unitvec * t_untivec, with s=u,v and t=x,y,z
97  jac[0][0] = _U[0];
98  jac[1][0] = _U[1];
99  jac[2][0] = _U[2];
100  jac[0][1] = _V[0];
101  jac[1][1] = _V[1];
102  jac[2][1] = _V[2];
103 
104  TMatrixT<Double_t> jac_t = jac.T();
105  TMatrixT<Double_t> jac_orig = jac;
106 
107  TMatrixT<Double_t> result = jac_t * (rawCov * jac_orig);
108  //std::cout << "hitCov="<<std::endl;
109  //result.Print();
110  return result;
111 }
112 
113 TMatrixT<Double_t> genf::GFSpacepointHitPolicy::hitCov(GFAbsRecoHit* hit,
114  const GFDetPlane& plane,
115  const GFDetPlane& planePrev,
116  const TMatrixT<Double_t>& state,
117  const Double_t& mass)
118 {
119  TVector3 _U;
120  TVector3 _V;
121 
122  _U = plane.getU();
123  _V = plane.getV();
124 
125  // Double_t eps(1.0e-6);
126  TMatrixT<Double_t> rawCov = hit->getRawHitCov();
127 
128  //rawCov[0][0] = 12.e12; // Force this. EC, 3-Apr-2012.
129  std::vector<double> tmpRawCov;
130  tmpRawCov.push_back(rawCov[0][0]);
131  tmpRawCov.push_back(rawCov[1][1]);
132  tmpRawCov.push_back(rawCov[2][2]);
133  tmpRawCov.push_back(rawCov[2][1]); // y-z correlated cov element.
134 
135  static std::vector<double> oldRawCov(tmpRawCov);
136  static std::vector<double> oldOldRawCov(tmpRawCov);
137  static GFDetPlane planePrevPrev(planePrev);
138  rawCov.ResizeTo(7, 7); // this becomes used now for something else entirely:
139  // the 7x7 raw errors on x,y,z,px,py,pz,th, which sandwiched between 5x7
140  // Jacobian converts it to the cov matrix for the 5x5 state space.
141  rawCov[0][0] = tmpRawCov[0]; // x
142  rawCov[1][1] = tmpRawCov[1]; // y
143  rawCov[2][2] = tmpRawCov[2]; // z
144  rawCov[1][2] = tmpRawCov[3]; // yz
145  rawCov[2][1] = tmpRawCov[3]; // yz
146 
147  // This Jacobian concerns the transfrom from planar to detector coords.
148  // TMatrixT<Double_t> jac(3,2);
149  TMatrixT<Double_t> jac(7, 5); // X,Y,Z,UX,UY,UZ,Theta in detector coords
150 
151  // 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
152 
153  // More correctly :
154  Double_t dist(0.3); // place holder!!
155  Double_t C;
156  //Double_t p (2.5); // place holder!!!
157  //p = abs(1/state[0][0]);
158 
159  Double_t mom = fabs(1.0 / state[0][0]);
160  Double_t beta = mom / sqrt(mass * mass + mom * mom);
161  dist = (plane.getO() - planePrev.getO()).Mag();
162  C = 0.0136 / beta * sqrt(dist / 14.0) * (1 + 0.038 * log(dist / 14.0));
163  TVector3 _D(plane.getNormal());
164 
165  // px^ = (x1-x2)/d;
166  // (sigpx)^2 = (sig_x1^2+sig_x2^2)*[(y1-y2)^2 + (z1-z2)^2]/d^4 +
167  // (x1-x2)^2*(y1-y2)^2*sigma_y1^2/d^6 + ...
168  // (x1-x2)^2*(z1-z2)^2*sigma_z2^2/d^6 +
169 
170  rawCov[3][3] =
171  (oldRawCov[0] + tmpRawCov[0]) / pow(dist, 4.) *
172  (pow((plane.getO() - planePrev.getO()).Y(), 2.) +
173  pow((plane.getO() - planePrev.getO()).Z(), 2.)) +
174  pow((plane.getO() - planePrev.getO()).X(), 2.) *
175  (pow((plane.getO() - planePrev.getO()).Y(), 2.) * (oldRawCov[1] + tmpRawCov[1]) +
176  pow((plane.getO() - planePrev.getO()).Z(), 2.) * (oldRawCov[2] + tmpRawCov[2])) /
177  pow(dist, 6.)
178  // y-z cross-term
179  + 3.0 * (plane.getO() - planePrev.getO()).X() * (plane.getO() - planePrev.getO()).Y() *
180  (plane.getO() - planePrev.getO()).Z() / pow(dist, 5.) * (tmpRawCov[3] + oldRawCov[3]);
181 
182  rawCov[4][4] =
183  (oldRawCov[1] + tmpRawCov[1]) / pow(dist, 4.) *
184  (pow((plane.getO() - planePrev.getO()).X(), 2.) +
185  pow((plane.getO() - planePrev.getO()).Z(), 2.)) +
186  pow((plane.getO() - planePrev.getO()).Y(), 2.) *
187  (pow((plane.getO() - planePrev.getO()).X(), 2.) * (oldRawCov[0] + tmpRawCov[0]) +
188  pow((plane.getO() - planePrev.getO()).Z(), 2.) * (oldRawCov[2] + tmpRawCov[2])) /
189  pow(dist, 6.)
190  // y-z cross-term
191  + 3.0 *
192  ((plane.getO() - planePrev.getO()).X() * (plane.getO() - planePrev.getO()).Y() *
193  (plane.getO() - planePrev.getO()).Z() / pow(dist, 5.) +
194  (plane.getO() - planePrev.getO()).Z() / pow(dist, 3.)) *
195  (tmpRawCov[3] + oldRawCov[3]);
196 
197  rawCov[5][5] =
198  (oldRawCov[2] + tmpRawCov[2]) / pow(dist, 4.) *
199  (pow((plane.getO() - planePrev.getO()).X(), 2.) +
200  pow((plane.getO() - planePrev.getO()).Y(), 2.)) +
201  pow((plane.getO() - planePrev.getO()).Z(), 2.) *
202  (pow((plane.getO() - planePrev.getO()).Y(), 2.) * (oldRawCov[0] + tmpRawCov[0]) +
203  pow((plane.getO() - planePrev.getO()).X(), 2.) * (oldRawCov[1] + tmpRawCov[1])) /
204  pow(dist, 6.) +
205  // y-z cross-term
206  3.0 *
207  ((plane.getO() - planePrev.getO()).X() * (plane.getO() - planePrev.getO()).Y() *
208  (plane.getO() - planePrev.getO()).Z() / pow(dist, 5.) +
209  (plane.getO() - planePrev.getO()).Y() / pow(dist, 3.)) *
210  (tmpRawCov[3] + oldRawCov[3]);
211 
212  Double_t num = (plane.getO() - planePrev.getO()) * (planePrev.getO() - planePrevPrev.getO());
213  Double_t d1 = (plane.getO() - planePrev.getO()).Mag(); // same as dist above
214  Double_t d2 = (planePrev.getO() - planePrevPrev.getO()).Mag();
215 
216  // This is the error on cosTh^2. There are 9 terms for diagonal errors, one for each
217  // x,y,z coordinate at 3 points. Below the first 3 terms are at pt 1.
218  // Second 3 are at 3. Third 3 are at 2, which is slightly more complicated.
219  // There are three more terms for non-diagonal erros.
220  double dcosTh =
221  pow(((planePrev.getO() - planePrevPrev.getO()).X() * d1 * d2 -
222  num * (plane.getO() - planePrev.getO()).X() * d2 / d1) /
223  pow(d1, 2.0) / pow(d2, 2.0),
224  2.0) *
225  tmpRawCov[0] +
226  pow(((planePrev.getO() - planePrevPrev.getO()).Y() * d1 * d2 +
227  num * (plane.getO() - planePrev.getO()).Y() * d2 / d1) /
228  pow(d1, 2.0) / pow(d2, 2.0),
229  2.0) *
230  tmpRawCov[1] +
231  pow(((planePrev.getO() - planePrevPrev.getO()).Z() * d1 * d2 +
232  num * (plane.getO() - planePrev.getO()).Z() * d2 / d1) /
233  pow(d1, 2.0) / pow(d2, 2.0),
234  2.0) *
235  tmpRawCov[2] +
236 
237  pow(((plane.getO() - planePrev.getO()).X() * d1 * d2 -
238  num * (planePrev.getO() - planePrevPrev.getO()).X() * d1 / d2) /
239  pow(d1, 2.0) / pow(d2, 2.0),
240  2.0) *
241  oldOldRawCov[0] +
242  pow(((plane.getO() - planePrev.getO()).Y() * d1 * d2 -
243  num * (planePrev.getO() - planePrevPrev.getO()).Y() * d1 / d2) /
244  pow(d1, 2.0) / pow(d2, 2.0),
245  2.0) *
246  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),
250  2.0) *
251  oldOldRawCov[2] +
252 
253  pow(((plane.getO() - planePrevPrev.getO()).X() * d1 * d2 -
254  num * (plane.getO() - planePrev.getO()).X() * d2 / d1 +
255  num * (planePrev.getO() - planePrevPrev.getO()).X() * d1 / d2) /
256  pow(d1, 2.0) / pow(d2, 2.0),
257  2.0) *
258  oldRawCov[0] +
259  pow(((plane.getO() - planePrevPrev.getO()).Y() * d1 * d2 -
260  num * (plane.getO() - planePrev.getO()).Y() * d2 / d1 +
261  num * (planePrev.getO() - planePrevPrev.getO()).Y() * d1 / d2) /
262  pow(d1, 2.0) / pow(d2, 2.0),
263  2.0) *
264  oldRawCov[1] +
265  pow(((plane.getO() - planePrevPrev.getO()).Z() * d1 * d2 -
266  num * (plane.getO() - planePrev.getO()).Z() * d2 / d1 +
267  num * (planePrev.getO() - planePrevPrev.getO()).Z() * d1 / d2) /
268  pow(d1, 2.0) / pow(d2, 2.0),
269  2.0) *
270  oldRawCov[2]
271  // And now the off-diagonal terms. This is a mess unto itself.
272  // First, d^2(costh)/d(z1)d(y1)
273  + ((plane.getO() - planePrev.getO()).Z() * (planePrev.getO() - planePrevPrev.getO()).Y() /
274  pow(d1, 3.) / d2 -
275  (plane.getO() - planePrev.getO()).Y() * (planePrev.getO() - planePrevPrev.getO()).Z() /
276  pow(d1, 3.) / d2 +
277  3. * (plane.getO() - planePrev.getO()).Y() * (plane.getO() - planePrev.getO()).Z() * num /
278  pow(d1, 5.) / d2) *
279  tmpRawCov[3]
280  // Next, d^2(costh)/d(z3)d(y3). Above w z1<->z3, y1<->y3
281  + ((planePrevPrev.getO() - planePrev.getO()).Z() * (planePrev.getO() - plane.getO()).Y() /
282  pow(d1, 3.) / d2 -
283  (planePrevPrev.getO() - planePrev.getO()).Y() * (planePrev.getO() - plane.getO()).Z() /
284  pow(d1, 3.) / d2 +
285  3. * (planePrevPrev.getO() - planePrev.getO()).Y() *
286  (planePrevPrev.getO() - planePrev.getO()).Z() * num / pow(d1, 5.) / d2) *
287  oldOldRawCov[3]
288  // Last, d^2(costh)/d(z2)d(y2). This is an even bigger mess.
289  + (((plane.getO() - planePrev.getO()).Y() - (planePrev.getO() - planePrevPrev.getO()).Y()) *
290  (-(plane.getO() - planePrev.getO()).Z() / pow(d1, 3.) / d2 +
291  (planePrev.getO() - planePrevPrev.getO()).Z() / pow(d2, 3.) / d1) +
292  (plane.getO() - planePrev.getO()).Y() *
293  ((-(planePrev.getO() - planePrevPrev.getO()).Z() + (plane.getO() - planePrev.getO()).Z()) /
294  pow(d1, 3.) / d2 -
295  num * (-3. * (plane.getO() - planePrev.getO()).Z() * d1 * d2 +
296  (planePrev.getO() - planePrevPrev.getO()).Z() * pow(d1, 3.) / d2)) /
297  pow(d1, 6.) / pow(d2, 2.) -
298  (planePrev.getO() - planePrevPrev.getO()).Y() *
299  ((-(planePrev.getO() - planePrevPrev.getO()).Z() + (plane.getO() - planePrev.getO()).Z()) /
300  pow(d2, 3.) / d1 -
301  num * (3. * (planePrev.getO() - planePrevPrev.getO()).Z() * d1 * d2 -
302  (plane.getO() - planePrev.getO()).Z() * pow(d2, 3.) / d1)) /
303  pow(d1, 2.) / pow(d2, 6.)) *
304  oldRawCov[3];
305 
306  // That was delta(cos(theta))^2. I want delta(theta)^2. dTh^2 = d(cosTh)^2/sinTh^2
307  Double_t theta(TMath::ACos((plane.getO() - planePrev.getO()).Unit() *
308  (planePrev.getO() - planePrevPrev.getO()).Unit()));
309  rawCov[6][6] =
310  TMath::Min(pow(dcosTh, 2.) / pow(TMath::Sin(theta), 2.), pow(TMath::Pi() / 2.0, 2.));
311 
312  // This means I'm too close to the endpoints to have histories that allow
313  // proper calculation of above rawCovs. Use below defaults instead.
314  if (d1 == 0 || d2 == 0) {
315  rawCov[3][3] = pow(0.2, 2.0); // Unit Px
316  rawCov[4][4] = pow(0.2, 2.0); // Unit Py
317  rawCov[5][5] = pow(0.2, 2.0); // Unit Pz
318  rawCov[6][6] = pow(0.1, 2.0); // theta. 0.3/3mm, say.
319  dist = 0.3;
320  C = 0.0136 / beta * sqrt(dist / 14.0) * (1 + 0.038 * log(dist / 14.0));
321  }
322 
323  // This forces a huge/tiny theta error, which effectively freezes/makes-us-sensitive theta, as is.
324  //rawCov[6][6] = /*9.99e9*/ 0.1;
325 
326  TVector3 u = plane.getU();
327  TVector3 v = plane.getV();
328  TVector3 w = u.Cross(v);
329 
330  // Below line has been done, with code change in GFKalman that has updated
331  // the plane orientation by now.
332  // TVector3 pTilde = 1.0 * (w + state[1][0] * u + state[2][0] * v);
333  TVector3 pTilde = w;
334  double pTildeMag = pTilde.Mag();
335 
336  jac.Zero();
337  jac[6][0] = 1. / C; // Should be 1/C?; Had 1 until ... 12-Feb-2013
338 
339  jac[0][3] = _U[0];
340  jac[1][3] = _U[1];
341  jac[2][3] = _U[2];
342 
343  jac[0][4] = _V[0];
344  jac[1][4] = _V[1];
345  jac[2][4] = _V[2];
346 
347  // cnp'd from RKTrackRep.cxx, line ~496
348  // da{x,y,z}/du'
349  jac[3][1] = 1.0 / pTildeMag * (u.X() - pTilde.X() / (pTildeMag * pTildeMag) * u * pTilde);
350  jac[4][1] = 1.0 / pTildeMag * (u.Y() - pTilde.Y() / (pTildeMag * pTildeMag) * u * pTilde);
351  jac[5][1] = 1.0 / pTildeMag * (u.Z() - pTilde.Z() / (pTildeMag * pTildeMag) * u * pTilde);
352  // da{x,y,z}/dv'
353  jac[3][2] = 1.0 / pTildeMag * (v.X() - pTilde.X() / (pTildeMag * pTildeMag) * v * pTilde);
354  jac[4][2] = 1.0 / pTildeMag * (v.Y() - pTilde.Y() / (pTildeMag * pTildeMag) * v * pTilde);
355  jac[5][2] = 1.0 / pTildeMag * (v.Z() - pTilde.Z() / (pTildeMag * pTildeMag) * v * pTilde);
356 
357  TMatrixT<Double_t> jac_orig = jac;
358  TMatrixT<Double_t> jac_t = jac.T();
359 
360  TMatrixT<Double_t> result = jac_t * (rawCov * jac_orig);
361  //std::cout << "hitCov="<<std::endl;
362  //result.Print();
363 
364  oldOldRawCov = oldRawCov;
365  oldRawCov = tmpRawCov;
366  planePrevPrev = planePrev;
367 
368  return result;
369 }
370 
371 const genf::GFDetPlane& genf::GFSpacepointHitPolicy::detPlane(GFAbsRecoHit* hit, GFAbsTrackRep* rep)
372 {
373  TMatrixT<Double_t> rawcoord = hit->getRawHitCoord();
374  TVector3 point(rawcoord[0][0], rawcoord[1][0], rawcoord[2][0]);
375 
376  TVector3 poca, dirInPoca;
377  rep->extrapolateToPoint(point, poca, dirInPoca);
378 
379  fPlane.setO(point);
380  fPlane.setNormal(dirInPoca);
381 
382  return fPlane;
383 }
384 
385 //ClassImp(GFSpacepointHitPolicy)
Y
Float_t Y
Definition: plot.C:37
genf::GFDetPlane::setO
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:94
genf::GFAbsRecoHit::getRawHitCoord
TMatrixT< Double_t > getRawHitCoord() const
Get raw hit coordinates.
Definition: GFAbsRecoHit.h:185
genf::GFSpacepointHitPolicy::hitCoord
TMatrixT< double > hitCoord(GFAbsRecoHit *, const GFDetPlane &)
Hit coordinates in detector plane.
Definition: GFSpacepointHitPolicy.cxx:28
genf::GFAbsTrackRep
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:81
Z
Float_t Z
Definition: plot.C:37
genf::GFDetPlane::getNormal
TVector3 getNormal() const
Definition: GFDetPlane.cxx:133
genf::GFSpacepointHitPolicy::hitCov
TMatrixT< double > hitCov(GFAbsRecoHit *, const GFDetPlane &)
Hit covariances in detector plane.
Definition: GFSpacepointHitPolicy.cxx:83
genf::GFDetPlane::getO
TVector3 getO() const
Definition: GFDetPlane.h:77
genf::GFSpacepointHitPolicy::fPolicyName
static const std::string fPolicyName
Definition: GFSpacepointHitPolicy.h:85
genf::GFAbsTrackRep::extrapolateToPoint
virtual void extrapolateToPoint(const TVector3 &point, TVector3 &poca, TVector3 &normVec)
This method is to extrapolate the track to point of closest approach to a point in space...
Definition: GFAbsTrackRep.cxx:80
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
genf::GFSpacepointHitPolicy::detPlane
const GFDetPlane & detPlane(GFAbsRecoHit *, GFAbsTrackRep *)
Get detector plane perpendicular to track.
Definition: GFSpacepointHitPolicy.cxx:371
genf::GFDetPlane::getU
TVector3 getU() const
Definition: GFDetPlane.h:78
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:537
genf::GFDetPlane::setNormal
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:151
genf::GFDetPlane::getV
TVector3 getV() const
Definition: GFDetPlane.h:79
X
Float_t X
Definition: plot.C:37
genf::GFAbsRecoHit::getRawHitCov
TMatrixT< Double_t > getRawHitCov() const
Get raw hit covariances.
Definition: GFAbsRecoHit.h:180
w
Float_t w
Definition: plot.C:20