LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Attributes | List of all members
trkf::VertexFitAlg Class Reference

#include "VertexFitAlg.h"

Public Member Functions

void VertexFit (std::vector< std::vector< geo::WireID >> const &hitWID, std::vector< std::vector< double >> const &hitX, std::vector< std::vector< double >> const &hitXErr, TVector3 &VtxPos, TVector3 &VtxPosErr, std::vector< TVector3 > &TrkDir, std::vector< TVector3 > &TrkDirErr, float &ChiDOF) const
 

Static Public Member Functions

static void fcnVtxPos (Int_t &, Double_t *, Double_t &fval, double *par, Int_t flag)
 

Static Public Attributes

static VertexFitMinuitStruct fVtxFitMinStr
 

Private Attributes

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

Detailed Description

Definition at line 28 of file VertexFitAlg.h.

Member Function Documentation

void trkf::VertexFitAlg::fcnVtxPos ( Int_t &  ,
Double_t *  ,
Double_t &  fval,
double *  par,
Int_t  flag 
)
static

Definition at line 27 of file VertexFitAlg.cxx.

References util::abs(), VertexFitMinuitStruct::ChiDoF, VertexFitMinuitStruct::DoF, E, VertexFitMinuitStruct::FirstWire, fVtxFitMinStr, VertexFitMinuitStruct::HitX, VertexFitMinuitStruct::HitXErr, VertexFitMinuitStruct::OrthY, VertexFitMinuitStruct::OrthZ, VertexFitMinuitStruct::Plane, VertexFitMinuitStruct::Wire, and VertexFitMinuitStruct::WirePitch.

Referenced by VertexFit().

28  {
29  // Minuit function for fitting the vertex position and vertex track directions
30 
31  fval = 0;
32  double vWire = 0, DirX, DirY, DirZ, DirU, dX, dU, arg;
33  unsigned short ipl, lastpl, indx;
34 
35  for (unsigned short itk = 0; itk < fVtxFitMinStr.HitX.size(); ++itk) {
36  lastpl = 4;
37  // index of the track Y direction vector. Z direction is the next one
38  indx = 3 + 2 * itk;
39  for (unsigned short iht = 0; iht < fVtxFitMinStr.HitX[itk].size(); ++iht) {
40  ipl = fVtxFitMinStr.Plane[itk][iht];
41  if (ipl != lastpl) {
42  // get the vertex position in this plane
43  // vertex wire number in the Detector coordinate system (equivalent to WireCoordinate)
44  //vtx wir = vtx Y * OrthY + vtx Z * OrthZ - wire offset
45  vWire = par[1] * fVtxFitMinStr.OrthY[ipl] + par[2] * fVtxFitMinStr.OrthZ[ipl] -
46  fVtxFitMinStr.FirstWire[ipl];
47  lastpl = ipl;
48  } // ipl != lastpl
49  DirY = par[indx];
50  DirZ = par[indx + 1];
51  // rotate the track direction DirY, DirZ into the wire coordinate of this plane. The OrthVectors in ChannelMapStandardAlg
52  // are divided by the wire pitch so we need to correct for that here
53  DirU = fVtxFitMinStr.WirePitch *
54  (DirY * fVtxFitMinStr.OrthY[ipl] + DirZ * fVtxFitMinStr.OrthZ[ipl]);
55  // distance (cm) between the wire and the vertex in the wire coordinate system (U)
56  dU = fVtxFitMinStr.WirePitch * (fVtxFitMinStr.Wire[itk][iht] - vWire);
57  if (std::abs(DirU) < 1E-3 || std::abs(dU) < 1E-3) {
58  // vertex is on the wire
59  dX = par[0] - fVtxFitMinStr.HitX[itk][iht];
60  }
61  else {
62  // project from vertex to the wire. We need to find dX/dU so first find DirX
63  DirX = 1 - DirY * DirY - DirZ * DirZ;
64  // DirX should be > 0 but the bounds on DirY and DirZ are +/- 1 so it is possible for a non-physical result.
65  if (DirX < 0) DirX = 0;
66  DirX = sqrt(DirX);
67  // Get the DirX sign from the relative X position of the hit and the vertex
68  if (fVtxFitMinStr.HitX[itk][iht] < par[0]) DirX = -DirX;
69  dX = par[0] + (dU * DirX / DirU) - fVtxFitMinStr.HitX[itk][iht];
70  }
71  arg = dX / fVtxFitMinStr.HitXErr[itk][iht];
72  fval += arg * arg;
73  } // iht
74  } //itk
75 
76  fval /= fVtxFitMinStr.DoF;
77 
78  // save the final chisq/dof in the struct on the last call
79  if (flag == 3) fVtxFitMinStr.ChiDoF = fval;
80 
81  } // fcnVtxPos
VertexFitMinuitStruct::Wire
std::vector< std::vector< unsigned short > > Wire
Definition: VertexFitMinuitStruct.h:33
VertexFitMinuitStruct::HitX
std::vector< std::vector< double > > HitX
Definition: VertexFitMinuitStruct.h:29
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
VertexFitMinuitStruct::Plane
std::vector< std::vector< unsigned short > > Plane
Definition: VertexFitMinuitStruct.h:32
VertexFitMinuitStruct::OrthY
std::array< double, 3 > OrthY
Definition: VertexFitMinuitStruct.h:24
VertexFitMinuitStruct::HitXErr
std::vector< std::vector< double > > HitXErr
Definition: VertexFitMinuitStruct.h:31
E
Float_t E
Definition: plot.C:20
VertexFitMinuitStruct::FirstWire
std::array< double, 3 > FirstWire
Definition: VertexFitMinuitStruct.h:26
VertexFitMinuitStruct::OrthZ
std::array< double, 3 > OrthZ
Definition: VertexFitMinuitStruct.h:25
trkf::VertexFitAlg::fVtxFitMinStr
static VertexFitMinuitStruct fVtxFitMinStr
Definition: VertexFitAlg.h:40
void trkf::VertexFitAlg::VertexFit ( std::vector< std::vector< geo::WireID >> const &  hitWID,
std::vector< std::vector< double >> const &  hitX,
std::vector< std::vector< double >> const &  hitXErr,
TVector3 &  VtxPos,
TVector3 &  VtxPosErr,
std::vector< TVector3 > &  TrkDir,
std::vector< TVector3 > &  TrkDirErr,
float &  ChiDOF 
) const

Definition at line 85 of file VertexFitAlg.cxx.

References VertexFitMinuitStruct::ChiDoF, geo::CryostatID::Cryostat, VertexFitMinuitStruct::Cstat, VertexFitMinuitStruct::Dir, VertexFitMinuitStruct::DoF, fcnVtxPos(), VertexFitMinuitStruct::FirstWire, fVtxFitMinStr, geom, VertexFitMinuitStruct::HitX, VertexFitMinuitStruct::HitXErr, VertexFitMinuitStruct::NPlanes, geo::TPCGeo::Nplanes(), VertexFitMinuitStruct::OrthY, VertexFitMinuitStruct::OrthZ, VertexFitMinuitStruct::Plane, util::size(), VertexFitMinuitStruct::TPC, geo::TPCID::TPC, geo::GeometryCore::TPC(), VertexFitMinuitStruct::VtxPos, VertexFitMinuitStruct::Wire, geo::GeometryCore::WireCoordinate(), VertexFitMinuitStruct::WirePitch, and geo::GeometryCore::WirePitch().

Referenced by trkf::CCTrackMaker::FitVertices().

93  {
94  // The passed set of hit WireIDs, X positions and X errors associated with a Track
95  // are fitted to a vertex position VtxPos. The fitted track direction vectors trkDir, TrkDirErr
96  // and ChiDOF are returned to the calling routine
97 
98  // assume failure
99  ChiDOF = 9999;
100 
101  // need at least hits for two tracks
102  if (hitX.size() < 2) return;
103  if (hitX.size() != hitWID.size()) return;
104  if (hitX.size() != hitXErr.size()) return;
105  if (hitX.size() != TrkDir.size()) return;
106 
107  // number of variables = 3 for the vertex position + 2 * number of track directions
108  const unsigned int ntrks = hitX.size();
109  const unsigned int npars = 3 + 2 * ntrks;
110  unsigned int npts = 0;
111  for (unsigned int itk = 0; itk < ntrks; ++itk)
112  npts += hitX[itk].size();
113 
114  if (npts < ntrks) return;
115 
116  // Get the cryostat and tpc from the first hit
117  geo::TPCID const& tpcid = hitWID[0][0];
118  unsigned int const nplanes = geom->TPC(tpcid).Nplanes();
119 
120  fVtxFitMinStr.Cstat = tpcid.Cryostat;
121  fVtxFitMinStr.TPC = tpcid.TPC;
122  fVtxFitMinStr.NPlanes = nplanes;
123  fVtxFitMinStr.WirePitch = geom->WirePitch(hitWID[0][0]);
124 
125  // Put geometry conversion factors into the struct
126  for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
127  geo::PlaneID const planeID{tpcid, ipl};
128  fVtxFitMinStr.FirstWire[ipl] = -geom->WireCoordinate(geo::Point_t{0, 0, 0}, planeID);
129  fVtxFitMinStr.OrthY[ipl] =
130  geom->WireCoordinate(geo::Point_t{0, 1, 0}, planeID) + fVtxFitMinStr.FirstWire[ipl];
131  fVtxFitMinStr.OrthZ[ipl] =
132  geom->WireCoordinate(geo::Point_t{0, 0, 1}, planeID) + fVtxFitMinStr.FirstWire[ipl];
133  }
134  // and the vertex starting position
135  fVtxFitMinStr.VtxPos = VtxPos;
136 
137  // and the track direction and hits
138  fVtxFitMinStr.HitX = hitX;
139  fVtxFitMinStr.HitXErr = hitXErr;
140  fVtxFitMinStr.Plane.resize(ntrks);
141  fVtxFitMinStr.Wire.resize(ntrks);
142  for (unsigned int itk = 0; itk < ntrks; ++itk) {
143  fVtxFitMinStr.Plane[itk].resize(hitX[itk].size());
144  fVtxFitMinStr.Wire[itk].resize(hitX[itk].size());
145  for (std::size_t iht = 0; iht < hitWID[itk].size(); ++iht) {
146  fVtxFitMinStr.Plane[itk][iht] = hitWID[itk][iht].Plane;
147  fVtxFitMinStr.Wire[itk][iht] = hitWID[itk][iht].Wire;
148  }
149  } // itk
150  fVtxFitMinStr.Dir = TrkDir;
151 
152  fVtxFitMinStr.DoF = npts - npars;
153 
154  // define the starting parameters
155  std::vector<double> par(npars);
156  std::vector<double> stp(npars);
157  std::vector<double> parerr(npars);
158 
159  TMinuit* gMin = new TMinuit(npars);
160  gMin->SetFCN(VertexFitAlg::fcnVtxPos);
161  int errFlag = 0;
162  double arglist[10];
163 
164  // print level (-1 = none, 1 = yes)
165  arglist[0] = -1;
166  //
167  // ***** remember to delete gMin before returning on an error
168  //
169  gMin->mnexcm("SET PRINT", arglist, 1, errFlag);
170 
171  // the vertex position
172  for (unsigned int ipar = 0; ipar < 3u; ++ipar) {
173  par[ipar] = fVtxFitMinStr.VtxPos[ipar]; // in cm
174  stp[ipar] = 0.1; // 1 mm initial step
175  gMin->mnparm(ipar, "", par[ipar], stp[ipar], -1E6, 1E6, errFlag);
176  }
177  // use Y, Z track directions. There is no constraint that the direction vector is unit-normalized
178  // since we are only passing two of the components. Minuit could violate this requirement when
179  // fitting. fcnVtxPos prevents non-physical values and Minuit may throw error messages as a result.
180  for (unsigned int itk = 0; itk < ntrks; ++itk) {
181  unsigned int ipar = 3 + 2 * itk;
182  par[ipar] = fVtxFitMinStr.Dir[itk](1);
183  stp[ipar] = 0.03;
184  gMin->mnparm(ipar, "", par[ipar], stp[ipar], -1.05, 1.05, errFlag);
185  ++ipar;
186  par[ipar] = fVtxFitMinStr.Dir[itk](2);
187  stp[ipar] = 0.03;
188  gMin->mnparm(ipar, "", par[ipar], stp[ipar], -1.05, 1.05, errFlag);
189  } // itk
190 
191  // set strategy 0 for faster Minuit fitting
192  arglist[0] = 0.;
193  gMin->mnexcm("SET STRATEGY", arglist, 1, errFlag);
194 
195  // execute Minuit command: Migrad, max calls, tolerance
196  arglist[0] = 500; // max calls
197  arglist[1] = 1.; // tolerance on fval in fcn
198  gMin->mnexcm("MIGRAD", arglist, 2, errFlag);
199 
200  // call fcn to get final fit values
201  arglist[0] = 3;
202  gMin->mnexcm("CALL", arglist, 1, errFlag);
203  ChiDOF = fVtxFitMinStr.ChiDoF;
204 
205  // get the parameters
206  for (unsigned short ip = 0; ip < par.size(); ++ip) {
207  gMin->GetParameter(ip, par[ip], parerr[ip]);
208  }
209 
210  // return the vertex position and errors
211  for (unsigned int ipar = 0; ipar < 3u; ++ipar) {
212  VtxPos[ipar] = par[ipar];
213  VtxPosErr[ipar] = parerr[ipar];
214  }
215  // return the track directions and the direction errors if applicable
216  bool returnTrkDirErrs = (TrkDirErr.size() == TrkDir.size());
217  for (unsigned int itk = 0; itk < ntrks; ++itk) {
218  unsigned int const ipar = 3 + 2 * itk;
219  double arg = 1 - par[ipar] * par[ipar] - par[ipar + 1] * par[ipar + 1];
220  if (arg < 0) arg = 0;
221  TrkDir[itk](0) = sqrt(arg);
222  TrkDir[itk](1) = par[ipar];
223  TrkDir[itk](2) = par[ipar + 1];
224  if (returnTrkDirErrs) {
225  // TODO This needs to be checked if there is a need for trajectory errors
226  double errY = parerr[ipar] / par[ipar];
227  double errZ = parerr[ipar + 1] / par[ipar + 1];
228  TrkDirErr[itk](0) = sqrt(arg * (errY * errY + errZ * errZ));
229  TrkDirErr[itk](1) = parerr[ipar];
230  TrkDirErr[itk](2) = parerr[ipar + 1];
231  }
232  } // itk
233 
234  delete gMin;
235 
236  } // VertexFit()
geo::GeometryCore::WireCoordinate
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1093
VertexFitMinuitStruct::Wire
std::vector< std::vector< unsigned short > > Wire
Definition: VertexFitMinuitStruct.h:33
geo::TPCGeo::Nplanes
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
VertexFitMinuitStruct::HitX
std::vector< std::vector< double > > HitX
Definition: VertexFitMinuitStruct.h:29
VertexFitMinuitStruct::Plane
std::vector< std::vector< unsigned short > > Plane
Definition: VertexFitMinuitStruct.h:32
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
VertexFitMinuitStruct::OrthY
std::array< double, 3 > OrthY
Definition: VertexFitMinuitStruct.h:24
geo::GeometryCore::TPC
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
VertexFitMinuitStruct::Dir
std::vector< TVector3 > Dir
Definition: VertexFitMinuitStruct.h:35
VertexFitMinuitStruct::HitXErr
std::vector< std::vector< double > > HitXErr
Definition: VertexFitMinuitStruct.h:31
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
VertexFitMinuitStruct::FirstWire
std::array< double, 3 > FirstWire
Definition: VertexFitMinuitStruct.h:26
geo::Point_t
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
trkf::VertexFitAlg::fcnVtxPos
static void fcnVtxPos(Int_t &, Double_t *, Double_t &fval, double *par, Int_t flag)
Definition: VertexFitAlg.cxx:27
trkf::VertexFitAlg::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: VertexFitAlg.h:45
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
VertexFitMinuitStruct::OrthZ
std::array< double, 3 > OrthZ
Definition: VertexFitMinuitStruct.h:25
trkf::VertexFitAlg::fVtxFitMinStr
static VertexFitMinuitStruct fVtxFitMinStr
Definition: VertexFitAlg.h:40
geo::GeometryCore::WirePitch
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:763

Member Data Documentation

VertexFitMinuitStruct trkf::VertexFitAlg::fVtxFitMinStr
static

Definition at line 40 of file VertexFitAlg.h.

Referenced by fcnVtxPos(), and VertexFit().

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

Definition at line 45 of file VertexFitAlg.h.

Referenced by VertexFit().

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