Public Member Functions | Private Attributes | Static Private Attributes | List of all members
genf::GFSpacepointHitPolicy Class Reference

#include <GFSpacepointHitPolicy.h>

Inheritance diagram for genf::GFSpacepointHitPolicy:

Public Member Functions

 GFSpacepointHitPolicy ()
 
const GFDetPlane & detPlane (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.

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 383 of file GFSpacepointHitPolicy.cxx.

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

Definition at line 81 of file GFSpacepointHitPolicy.h.

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

Hit coordinates in detector plane.

Definition at line 30 of file GFSpacepointHitPolicy.cxx.

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

Definition at line 60 of file GFSpacepointHitPolicy.cxx.

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

Hit covariances in detector plane.

Definition at line 90 of file GFSpacepointHitPolicy.cxx.

91 {
92  TVector3 _U;
93  TVector3 _V;
94 
95  _U = plane.getU();
96  _V = plane.getV();
97 
98  TMatrixT<Double_t> rawCov = hit->getRawHitCov();
99  //rawCov[0][0] = 12.e12; // Force this. EC, 3-Apr-2012.
100 
101  TMatrixT<Double_t> jac(3,2);
102 
103  // jac = dF_i/dx_j = s_unitvec * t_untivec, with s=u,v and t=x,y,z
104  jac[0][0] = _U[0];
105  jac[1][0] = _U[1];
106  jac[2][0] = _U[2];
107  jac[0][1] = _V[0];
108  jac[1][1] = _V[1];
109  jac[2][1] = _V[2];
110 
111  TMatrixT<Double_t> jac_t = jac.T();
112  TMatrixT<Double_t> jac_orig = jac;
113 
114  TMatrixT<Double_t> result=jac_t * (rawCov * jac_orig);
115  //std::cout << "hitCov="<<std::endl;
116  //result.Print();
117  return result;
118 }
result
static QCString result
Definition: fortranscanner.cpp:56614
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
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 121 of file GFSpacepointHitPolicy.cxx.

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

Member Data Documentation

GFDetPlane genf::GFSpacepointHitPolicy::fPlane
private

Definition at line 87 of file GFSpacepointHitPolicy.h.

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

Definition at line 84 of file GFSpacepointHitPolicy.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11