Public Member Functions | Private Attributes | List of all members
trkf::TrackLineFitAlg Class Reference

#include <TrackLineFitAlg.h>

Public Member Functions

void TrkLineFit (std::vector< geo::WireID > &hitWID, std::vector< double > &hitX, std::vector< double > &hitXErr, double XOrigin, TVector3 &Pos, TVector3 &Dir, float &ChiDOF) const
 

Private Attributes

art::ServiceHandle< geo::Geometry const > geom
 

Detailed Description

Definition at line 25 of file TrackLineFitAlg.h.

Member Function Documentation

void trkf::TrackLineFitAlg::TrkLineFit ( std::vector< geo::WireID > &  hitWID,
std::vector< double > &  hitX,
std::vector< double > &  hitXErr,
double  XOrigin,
TVector3 &  Pos,
TVector3 &  Dir,
float &  ChiDOF 
) const

Definition at line 28 of file TrackLineFitAlg.cxx.

30  {
31  // Linear fit using X as the independent variable. Hits to be fitted
32  // are passed in the hits vector in a pair form (X, WireID). The
33  // fitted track position at XOrigin is returned in the Pos vector.
34  // The direction cosines are returned in the Dir vector.
35  //
36  // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
37  // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
38  // a matrix to calculate a track projected to a point at X, and v is
39  // a vector (Yo, Zo, dY/dX, dZ/dX).
40 
41  // assume failure
42  ChiDOF = -1;
43 
44  if(hitX.size() < 4) return;
45  if(hitX.size() != hitWID.size()) return;
46  if(hitX.size() != hitXErr.size()) return;
47 
48  const unsigned int nvars = 4;
49  unsigned int npts = hitX.size();
50 
51  TMatrixD A(npts, nvars);
52  // vector holding the Wire number
53  TVectorD w(npts);
54  unsigned short ninpl[3] = {0};
55  unsigned short nok = 0;
56  unsigned short iht;
57  unsigned int ipl, tpc, cstat;
58  double x, cw, sw, off, wght;
59  for(iht = 0; iht < hitX.size(); ++iht) {
60  cstat = hitWID[iht].Cryostat;
61  tpc = hitWID[iht].TPC;
62  ipl = hitWID[iht].Plane;
63  // get the wire plane offset
64  off = geom->WireCoordinate(0, 0, ipl, tpc, cstat);
65  // get the "cosine-like" component
66  cw = geom->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
67  // the "sine-like" component
68  sw = geom->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
69  x = hitX[iht] - XOrigin;
70  if(hitXErr[iht] > 0) {
71  wght = 1 / hitXErr[iht];
72  } else {
73  wght = 1;
74  }
75  A[iht][0] = wght * cw;
76  A[iht][1] = wght * sw;
77  A[iht][2] = wght * cw * x;
78  A[iht][3] = wght * sw * x;
79  w[iht] = wght * (hitWID[iht].Wire - off);
80  ++ninpl[ipl];
81  // need at least two points in a plane
82  if(ninpl[ipl] == 2) ++nok;
83  }
84 
85  // need at least 2 planes with at least two points
86  if(nok < 2) return;
87 
88  TDecompSVD svd(A);
89  bool ok;
90  TVectorD tVec = svd.Solve(w, ok);
91 
92  ChiDOF = 0;
93 
94  // not enough points to calculate Chisq
95  if(hitX.size() == 4) return;
96 
97  double ypr, zpr, diff;
98  for(iht = 0; iht < hitX.size(); ++iht) {
99  cstat = hitWID[iht].Cryostat;
100  tpc = hitWID[iht].TPC;
101  ipl = hitWID[iht].Plane;
102  off = geom->WireCoordinate(0, 0, ipl, tpc, cstat);
103  cw = geom->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
104  sw = geom->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
105  x = hitX[iht] - XOrigin;
106  ypr = tVec[0] + tVec[2] * x;
107  zpr = tVec[1] + tVec[3] * x;
108  if(hitXErr[iht] > 0) {
109  wght = 1 / hitXErr[iht];
110  } else {
111  wght = 1;
112  }
113  if(wght <= 0) wght = 1;
114  diff = (ypr * cw + zpr * sw - (hitWID[iht].Wire - off)) / wght;
115  ChiDOF += diff * diff;
116  }
117 
118  double werr2 = geom->WirePitch(0, tpc, cstat);
119  werr2 *= werr2;
120  ChiDOF /= werr2;
121  ChiDOF /= (double)(npts - 4);
122 
123  double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
124  Dir[0] = 1 / norm;
125  Dir[1] = tVec[2] / norm;
126  Dir[2] = tVec[3] / norm;
127 
128  Pos[0] = XOrigin;
129  Pos[1] = tVec[0];
130  Pos[2] = tVec[1];
131 /*
132  // covariance matrix
133  TMatrixD fV = svd.GetV();
134  PosCov(1, 1) = fV(0, 0);
135  PosCov(2, 1) = fV(1, 0);
136  PosCov(1, 2) = fV(0, 1);
137  PosCov(2, 2) = fV(1, 1);
138  // A conservative fake for Pos[0]
139  PosCov(0, 0) = PosCov(1,1);
140 */
141 
142  } // TrkLineFit()
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1128
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
trkf::TrackLineFitAlg::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackLineFitAlg.h:33
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:997
A
#define A
Definition: memgrp.cpp:38
train.x
list x
Definition: train.py:276

Member Data Documentation

art::ServiceHandle<geo::Geometry const> trkf::TrackLineFitAlg::geom
private

Definition at line 33 of file TrackLineFitAlg.h.

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