LArSoft  v06_85_00
Liquid Argon Software toolkit - http://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 > geom
 

Detailed Description

Definition at line 29 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 29 of file VertexFitAlg.cxx.

References 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().

30  {
31  // Minuit function for fitting the vertex position and vertex track directions
32 
33  fval = 0;
34  double vWire = 0, DirX, DirY, DirZ, DirU, dX, dU, arg;
35  unsigned short ipl, lastpl, indx;
36 
37  for(unsigned short itk = 0; itk < fVtxFitMinStr.HitX.size(); ++itk) {
38  lastpl = 4;
39  // index of the track Y direction vector. Z direction is the next one
40  indx = 3 + 2 * itk;
41  for(unsigned short iht = 0; iht < fVtxFitMinStr.HitX[itk].size(); ++iht) {
42  ipl = fVtxFitMinStr.Plane[itk][iht];
43  if(ipl != lastpl) {
44  // get the vertex position in this plane
45  // vertex wire number in the Detector coordinate system (equivalent to WireCoordinate)
46  //vtx wir = vtx Y * OrthY + vtx Z * OrthZ - wire offset
47  vWire = par[1] * fVtxFitMinStr.OrthY[ipl] + par[2] * fVtxFitMinStr.OrthZ[ipl] - fVtxFitMinStr.FirstWire[ipl];
48 // if(flag == 1) mf::LogVerbatim("VF")<<"fcn vtx "<<par[0]<<" "<<par[1]<<" "<<par[2]<<" vWire "<<vWire<<" OrthY "<<fVtxFitMinStr.OrthY[ipl]<<" OrthZ "<<fVtxFitMinStr.OrthZ[ipl];
49  lastpl = ipl;
50  } // ipl != lastpl
51  DirY = par[indx];
52  DirZ = par[indx + 1];
53  // rotate the track direction DirY, DirZ into the wire coordinate of this plane. The OrthVectors in ChannelMapStandardAlg
54  // are divided by the wire pitch so we need to correct for that here
55  DirU = fVtxFitMinStr.WirePitch * (DirY * fVtxFitMinStr.OrthY[ipl] + DirZ * fVtxFitMinStr.OrthZ[ipl]);
56  // distance (cm) between the wire and the vertex in the wire coordinate system (U)
57  dU = fVtxFitMinStr.WirePitch * (fVtxFitMinStr.Wire[itk][iht] - vWire);
58  if(std::abs(DirU) < 1E-3 || std::abs(dU) < 1E-3) {
59  // vertex is on the wire
60  dX = par[0] - fVtxFitMinStr.HitX[itk][iht];
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 // if(flag == 1) mf::LogVerbatim("VF")<<"fcn itk "<<itk<<" iht "<<iht<<" ipl "<<ipl<<" DirX "<<DirX<<" DirY "<<DirY<<" DirZ "<<DirZ
73 // <<" DirU "<<DirU<<" W "<<fVtxFitMinStr.Wire[itk][iht]<<" X "<<fVtxFitMinStr.HitX[itk][iht]<<" dU "<<dU<<" dX "<<dX<<" arg "<<arg;
74  fval += arg * arg;
75  } // iht
76  } //itk
77 
78  fval /= fVtxFitMinStr.DoF;
79 
80  // save the final chisq/dof in the struct on the last call
81  if(flag == 3) fVtxFitMinStr.ChiDoF = fval;
82 
83  } // fcnVtxPos
VertexFitMinuitStruct::Wire
std::vector< std::vector< unsigned short > > Wire
Definition: VertexFitMinuitStruct.h:28
E
Float_t E
Definition: plot.C:23
VertexFitMinuitStruct::HitX
std::vector< std::vector< double > > HitX
Definition: VertexFitMinuitStruct.h:24
VertexFitMinuitStruct::Plane
std::vector< std::vector< unsigned short > > Plane
Definition: VertexFitMinuitStruct.h:27
VertexFitMinuitStruct::OrthY
std::array< double, 3 > OrthY
Definition: VertexFitMinuitStruct.h:19
VertexFitMinuitStruct::FirstWire
std::array< double, 3 > FirstWire
Definition: VertexFitMinuitStruct.h:21
VertexFitMinuitStruct::HitXErr
std::vector< std::vector< double > > HitXErr
Definition: VertexFitMinuitStruct.h:26
VertexFitMinuitStruct::OrthZ
std::array< double, 3 > OrthZ
Definition: VertexFitMinuitStruct.h:20
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 87 of file VertexFitAlg.cxx.

References VertexFitMinuitStruct::ChiDoF, geo::GeometryCore::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, VertexFitMinuitStruct::TPC, geo::CryostatGeo::TPC(), VertexFitMinuitStruct::VtxPos, VertexFitMinuitStruct::Wire, geo::GeometryCore::WireCoordinate(), VertexFitMinuitStruct::WirePitch, and geo::GeometryCore::WirePitch().

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, itk;
111  for(itk = 0; itk < ntrks; ++itk) npts += hitX[itk].size();
112 
113  if(npts < ntrks) return;
114 
115  // Get the cryostat and tpc from the first hit
116  unsigned int cstat, tpc, nplanes, ipl, iht;
117  cstat = hitWID[0][0].Cryostat;
118  tpc = hitWID[0][0].TPC;
119  nplanes = geom->Cryostat(cstat).TPC(tpc).Nplanes();
120 
121  fVtxFitMinStr.Cstat = cstat;
122  fVtxFitMinStr.TPC = tpc;
123  fVtxFitMinStr.NPlanes = nplanes;
124  fVtxFitMinStr.WirePitch = geom->WirePitch(hitWID[0][0].Plane, tpc, cstat);
125 
126  // Put geometry conversion factors into the struct
127  for(ipl = 0; ipl < nplanes; ++ipl) {
128  fVtxFitMinStr.FirstWire[ipl] = -geom->WireCoordinate(0, 0, ipl, tpc, cstat);
129  fVtxFitMinStr.OrthY[ipl] = geom->WireCoordinate(1, 0, ipl, tpc, cstat) + fVtxFitMinStr.FirstWire[ipl];
130  fVtxFitMinStr.OrthZ[ipl] = geom->WireCoordinate(0, 1, ipl, tpc, cstat) + fVtxFitMinStr.FirstWire[ipl];
131  }
132  // and the vertex starting position
133  fVtxFitMinStr.VtxPos = VtxPos;
134 
135  // and the track direction and hits
136  fVtxFitMinStr.HitX = hitX;
137  fVtxFitMinStr.HitXErr = hitXErr;
138  fVtxFitMinStr.Plane.resize(ntrks);
139  fVtxFitMinStr.Wire.resize(ntrks);
140  for(itk = 0; itk < ntrks; ++itk) {
141  fVtxFitMinStr.Plane[itk].resize(hitX[itk].size());
142  fVtxFitMinStr.Wire[itk].resize(hitX[itk].size());
143  for(iht = 0; iht < hitWID[itk].size(); ++iht) {
144  fVtxFitMinStr.Plane[itk][iht] = hitWID[itk][iht].Plane;
145  fVtxFitMinStr.Wire[itk][iht] = hitWID[itk][iht].Wire;
146  }
147  } // itk
148  fVtxFitMinStr.Dir = TrkDir;
149 
150  fVtxFitMinStr.DoF = npts - npars;
151 
152  // define the starting parameters
153  std::vector<double> par(npars);
154  std::vector<double> stp(npars);
155  std::vector<double> parerr(npars);
156 
157  TMinuit *gMin = new TMinuit(npars);
158  gMin->SetFCN(VertexFitAlg::fcnVtxPos);
159  int errFlag = 0;
160  double arglist[10];
161 
162  // print level (-1 = none, 1 = yes)
163  arglist[0] = -1;
164  //
165  // ***** remember to delete gMin before returning on an error
166  //
167  gMin->mnexcm("SET PRINT", arglist, 1, errFlag);
168 
169  // the vertex position
170  unsigned short ipar;
171  for(ipar = 0; ipar < 3; ++ipar) {
172  par[ipar] = fVtxFitMinStr.VtxPos[ipar]; // in cm
173  stp[ipar] = 0.1; // 1 mm initial step
174  gMin->mnparm(ipar,"", par[ipar], stp[ipar], -1E6, 1E6, errFlag);
175  }
176  // use Y, Z track directions. There is no constraint that the direction vector is unit-normalized
177  // since we are only passing two of the components. Minuit could violate this requirement when
178  // fitting. fcnVtxPos prevents non-physical values and Minuit may throw error messages as a result.
179  for(itk = 0; itk < ntrks; ++itk) {
180  ipar = 3 + 2 * itk;
181  par[ipar] = fVtxFitMinStr.Dir[itk](1);
182  stp[ipar] = 0.03;
183  gMin->mnparm(ipar,"", par[ipar], stp[ipar], -1.05, 1.05, errFlag);
184  ++ipar;
185  par[ipar] = fVtxFitMinStr.Dir[itk](2);
186  stp[ipar] = 0.03;
187  gMin->mnparm(ipar,"", par[ipar], stp[ipar], -1.05, 1.05, errFlag);
188  } // itk
189 
190 /*
191  // Single call to fcnVtxPos for debugging it
192  std::cout<<"Starting: par ";
193  for(unsigned short ip = 0; ip < par.size(); ++ip) std::cout<<" "<<std::fixed<<std::setprecision(2)<<par[ip];
194  std::cout<<"\n";
195  // call fcn with starting parameters
196  arglist[0] = 1;
197  gMin->mnexcm("CALL", arglist, 1, errFlag);
198 */
199  // set strategy 0 for faster Minuit fitting
200  arglist[0] = 0.;
201  gMin->mnexcm("SET STRATEGY", arglist, 1, errFlag);
202 
203  // execute Minuit command: Migrad, max calls, tolerance
204  arglist[0] = 500; // max calls
205  arglist[1] = 1.; // tolerance on fval in fcn
206  gMin->mnexcm("MIGRAD", arglist, 2, errFlag);
207 
208  // call fcn to get final fit values
209  arglist[0] = 3;
210  gMin->mnexcm("CALL", arglist, 1, errFlag);
211  ChiDOF = fVtxFitMinStr.ChiDoF;
212 
213  // get the parameters
214  for(unsigned short ip = 0; ip < par.size(); ++ip) {
215  gMin->GetParameter(ip, par[ip], parerr[ip]);
216  }
217 
218  // return the vertex position and errors
219  for(ipar = 0; ipar < 3; ++ipar) {
220  VtxPos[ipar] = par[ipar];
221  VtxPosErr[ipar] = parerr[ipar];
222  }
223  // return the track directions and the direction errors if applicable
224  bool returnTrkDirErrs = (TrkDirErr.size() == TrkDir.size());
225  for(itk = 0; itk < ntrks; ++itk) {
226  ipar = 3 + 2 * itk;
227  double arg = 1 - par[ipar] * par[ipar] - par[ipar + 1] * par[ipar + 1];
228  if(arg < 0) arg = 0;
229  TrkDir[itk](0) = sqrt(arg);
230  TrkDir[itk](1) = par[ipar];
231  TrkDir[itk](2) = par[ipar + 1];
232  if(returnTrkDirErrs) {
233  // TODO This needs to be checked if there is a need for trajectory errors
234  double errY = parerr[ipar] / par[ipar];
235  double errZ = parerr[ipar + 1] / par[ipar + 1];
236  TrkDirErr[itk](0) = sqrt(arg * (errY * errY + errZ * errZ));
237  TrkDirErr[itk](1) = parerr[ipar];
238  TrkDirErr[itk](2) = parerr[ipar + 1];
239  }
240  } // itk
241 
242  delete gMin;
243 
244  } // VertexFit()
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:1066
VertexFitMinuitStruct::Wire
std::vector< std::vector< unsigned short > > Wire
Definition: VertexFitMinuitStruct.h:28
trkf::VertexFitAlg::geom
art::ServiceHandle< geo::Geometry > geom
Definition: VertexFitAlg.h:46
geo::TPCGeo::Nplanes
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:145
VertexFitMinuitStruct::HitX
std::vector< std::vector< double > > HitX
Definition: VertexFitMinuitStruct.h:24
VertexFitMinuitStruct::Plane
std::vector< std::vector< unsigned short > > Plane
Definition: VertexFitMinuitStruct.h:27
VertexFitMinuitStruct::OrthY
std::array< double, 3 > OrthY
Definition: VertexFitMinuitStruct.h:19
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:751
VertexFitMinuitStruct::Dir
std::vector< TVector3 > Dir
Definition: VertexFitMinuitStruct.h:30
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:280
VertexFitMinuitStruct::FirstWire
std::array< double, 3 > FirstWire
Definition: VertexFitMinuitStruct.h:21
trkf::VertexFitAlg::fcnVtxPos
static void fcnVtxPos(Int_t &, Double_t *, Double_t &fval, double *par, Int_t flag)
Definition: VertexFitAlg.cxx:29
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:174
VertexFitMinuitStruct::HitXErr
std::vector< std::vector< double > > HitXErr
Definition: VertexFitMinuitStruct.h:26
VertexFitMinuitStruct::OrthZ
std::array< double, 3 > OrthZ
Definition: VertexFitMinuitStruct.h:20
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
trkf::VertexFitAlg::fVtxFitMinStr
static VertexFitMinuitStruct fVtxFitMinStr
Definition: VertexFitAlg.h:40

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> trkf::VertexFitAlg::geom
private

Definition at line 46 of file VertexFitAlg.h.

Referenced by VertexFit().

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