LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
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 27 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.

References geom, norm, sw, w, geo::GeometryCore::WireCoordinate(), geo::GeometryCore::WirePitch(), and x.

Referenced by trkf::TrackTrajectoryAlg::TrackTrajectory().

35  {
36  // Linear fit using X as the independent variable. Hits to be fitted
37  // are passed in the hits vector in a pair form (X, WireID). The
38  // fitted track position at XOrigin is returned in the Pos vector.
39  // The direction cosines are returned in the Dir vector.
40  //
41  // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
42  // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
43  // a matrix to calculate a track projected to a point at X, and v is
44  // a vector (Yo, Zo, dY/dX, dZ/dX).
45 
46  // assume failure
47  ChiDOF = -1;
48 
49  if (hitX.size() < 4) return;
50  if (hitX.size() != hitWID.size()) return;
51  if (hitX.size() != hitXErr.size()) return;
52 
53  const unsigned int nvars = 4;
54  unsigned int npts = hitX.size();
55 
56  TMatrixD A(npts, nvars);
57  // vector holding the Wire number
58  TVectorD w(npts);
59  unsigned short ninpl[3] = {0};
60  unsigned short nok = 0;
61  for (std::size_t iht = 0; iht < hitX.size(); ++iht) {
62  auto const& wid = hitWID[iht];
63  // get the wire plane offset
64  double const off = geom->WireCoordinate(geo::Point_t{0, 0, 0}, wid);
65  // get the "cosine-like" component
66  double const cw = geom->WireCoordinate(geo::Point_t{0, 1, 0}, wid) - off;
67  // the "sine-like" component
68  double const sw = geom->WireCoordinate(geo::Point_t{0, 0, 1}, wid) - off;
69  double const x = hitX[iht] - XOrigin;
70  double wght{1.};
71  if (hitXErr[iht] > 0) { wght = 1 / hitXErr[iht]; }
72  A[iht][0] = wght * cw;
73  A[iht][1] = wght * sw;
74  A[iht][2] = wght * cw * x;
75  A[iht][3] = wght * sw * x;
76  w[iht] = wght * (hitWID[iht].Wire - off);
77  unsigned int ipl = wid.Plane;
78  ++ninpl[ipl];
79  // need at least two points in a plane
80  if (ninpl[ipl] == 2) ++nok;
81  }
82 
83  // need at least 2 planes with at least two points
84  if (nok < 2) return;
85 
86  TDecompSVD svd(A);
87  bool ok;
88  TVectorD tVec = svd.Solve(w, ok);
89 
90  ChiDOF = 0;
91 
92  // not enough points to calculate Chisq
93  if (hitX.size() == 4) return;
94 
95  // FIXME: The 'tpc' and 'cstat' variables are suspect as they are
96  // updated for each iteration below, but only the last
97  // value is used in the geom->WirePitch(...) calculation
98  // below.
99  unsigned int tpc{-1u}, cstat{-1u};
100  for (std::size_t iht = 0; iht < hitX.size(); ++iht) {
101  auto const& wid = hitWID[iht];
102  tpc = wid.TPC;
103  cstat = wid.Cryostat;
104  double const off = geom->WireCoordinate(geo::Point_t{0, 0, 0}, wid);
105  double const cw = geom->WireCoordinate(geo::Point_t{0, 1, 0}, wid) - off;
106  double const sw = geom->WireCoordinate(geo::Point_t{0, 0, 1}, wid) - off;
107  double const x = hitX[iht] - XOrigin;
108  double const ypr = tVec[0] + tVec[2] * x;
109  double const zpr = tVec[1] + tVec[3] * x;
110  double wght{1.};
111  if (hitXErr[iht] > 0) { wght = 1 / hitXErr[iht]; }
112  assert(wght > 0.);
113  double const diff = (ypr * cw + zpr * sw - (hitWID[iht].Wire - off)) / wght;
114  ChiDOF += diff * diff;
115  }
116 
117  double werr2 = geom->WirePitch(geo::PlaneID{cstat, tpc, 0});
118  werr2 *= werr2;
119  ChiDOF /= werr2;
120  ChiDOF /= (double)(npts - 4);
121 
122  double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
123  Dir[0] = 1 / norm;
124  Dir[1] = tVec[2] / norm;
125  Dir[2] = tVec[3] / norm;
126 
127  Pos[0] = XOrigin;
128  Pos[1] = tVec[0];
129  Pos[2] = tVec[1];
130  } // TrkLineFit()
Float_t x
Definition: compare.C:6
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
art::ServiceHandle< geo::Geometry const > geom
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
Float_t sw
Definition: plot.C:20
Float_t norm
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.

Member Data Documentation

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

Definition at line 38 of file TrackLineFitAlg.h.

Referenced by TrkLineFit().


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