Public Member Functions | Private Member Functions | Private Attributes | List of all members
trk::KalmanGeoHelper Class Reference

A collection of geometry functions used by the KalmanTrack tracking algorithm. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-04-12/TrackFit/KalmanGeoHelper.h"

Public Member Functions

 KalmanGeoHelper (bool badChan)
 
virtual ~KalmanGeoHelper ()
 
double BestTrackPoint (double z, double v, double halfz, double halfv, std::array< double, 2 > &bestPos, double z0, double v0, double dz, double dv)
 Estimate the position on a line within a box. More...
 
int CountMissedCellsOnLine (geo::View_t view, double vi, double zi, double vf, double zf, int maxCells)
 Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm. More...
 
int CountMissedCellsOnLine (geo::View_t view, double vi, double zi, double vf, double zf, std::vector< geo::OfflineChan > &hitsonline, std::vector< double > &distance)
 Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm. More...
 
double CalcMissFrac (double dist, double w, geo::View_t view)
 Estimation of the mip missed hit fraction as a function of the dist through a cell and cell depth w. More...
 
int MatchTrajectoryToPlane (TVector3 traj, int minplane, int maxplane)
 Finds the plane best matched to a trajectory point This method takes a trajectory point and finds the plane that closest matches this position in either view between minplane and maxplane. More...
 
int MatchTrajectoryToPlaneInView (TVector3 traj, geo::View_t view, int minplane, int maxplane)
 Variation of MatchTrajectoryToPlane that limits the matching plane to be in the input view. More...
 

Private Member Functions

void SetupHelper ()
 
double Efficiency (double w, double length)
 
void fillprivategeo ()
 

Private Attributes

bool fBadChannels
 
art::ServiceHandle< geo::Geometryfgeom
 
art::ServiceHandle< chaninfo::BadChanListfbadc
 
bool ndos = false
 
bool nd = false
 
bool fd = false
 
bool tb = false
 
int fnmaxplane
 
double fclipperLength
 
double fwlow
 
double fwhigh
 
double flengthmax
 
TF1 effpar0
 
TF1 effpar1
 
TF1 effpar2
 
TF1 effpar3
 

Detailed Description

A collection of geometry functions used by the KalmanTrack tracking algorithm.

Definition at line 21 of file KalmanGeoHelper.h.

Constructor & Destructor Documentation

trk::KalmanGeoHelper::KalmanGeoHelper ( bool  badChan)

Definition at line 50 of file KalmanGeoHelper.cxx.

References fBadChannels, and SetupHelper().

51  {
52  fBadChannels = badChan;
53  SetupHelper();
54  }
trk::KalmanGeoHelper::~KalmanGeoHelper ( )
virtual

Definition at line 57 of file KalmanGeoHelper.cxx.

58  {
59  }

Member Function Documentation

double trk::KalmanGeoHelper::BestTrackPoint ( double  z,
double  v,
double  halfz,
double  halfv,
std::array< double, 2 > &  bestPos,
double  z0,
double  v0,
double  dz,
double  dv 
)

Estimate the position on a line within a box.

This method takes a line, defined by its start point z0,v0 and direction dz and dv, and 2d box, defined by its center z,v and half lengths halfz, halfv and fills the bestPos with the best estimation of the closest point in that box t the line. For lines that do not intersect the box, this is defined as the corner of the box closest to the line. For lines that do intersect the box, the best position is the the point on the line bisecting the intersection points of the box.

Parameters
bestPosA vector that is filled with the best position of the line in the box with the first entry corresponding to z and the second to v
Returns
Total intersection length of the line with the box

Definition at line 130 of file KalmanGeoHelper.cxx.

References dz, registry_explorer::v, and test::z.

Referenced by CountMissedCellsOnLine(), trk::KalmanTrack::FilterTracker(), and trk::KalmanTrack::FindTracks().

132  {
133  // bestPos is a vector holding the best possible track hit position in the cell: bestPos[0] = z bestPos[1] = x/y.
134  // output vector the minimum distance between the projected track and the best track hit position in the cell.
135  // This function will take the initial state and z position and try to find the most compatible track hit point with the given hit
136 
137  double slope = dv/dz;
138  // Get the projected track at the cell boundaries.
139  double zlo = z - halfz;
140  double zhi = z + halfz;
141  double vlo = v - halfv;
142  double vhi = v + halfv;
143  double v1 = slope*(zlo-z0)+v0;
144  double v2 = slope*(zhi-z0)+v0;
145  double z1 = (vlo-v0)/slope+z0;
146  double z2 = (vhi-v0)/slope+z0;
147  double deltaz = 0;
148  double deltav = 0;
149  // Make a decision as to how to look for best track point as to avoid dividing by zero
150  if(dz == 0){
151  // either the z is in the z range or else it never will be, either case best position is z0, v
152  bestPos[0] = z0;
153  bestPos[1] = v;
154  // Increment deltav accordingly if actually in the cell
155  if(z >= zlo && z <= zhi){ deltav = vhi - vlo; }
156  } // end of if dz == 0
157  else{
158  // Go through the various cases to see if the track projection intersects the cell hit or not
159  if(slope>0){
160  // check if track enters cell from left side
161  if(v1 >= vlo && v1 <= vhi){
162  // check if track exits on the right side of the cell
163  if(v2 <= vhi){
164  bestPos[0] = z;
165  bestPos[1] = 0.5*v1+0.5*v2;
166  deltaz = zhi-zlo;
167  deltav = v2-v1;
168  }
169  else{
170  // track exits through the top of the cell
171  bestPos[0] = 0.5*zlo+0.5*z2;
172  bestPos[1] = 0.5*v1 + 0.5*vhi;
173  deltaz = z2-zlo;
174  deltav = vhi-v1;
175  }
176  }
177  // check if track enters cell from bottom side
178  else if(v1 < vlo && v2 >= vlo){
179  // check if track exits on right side of cell
180  if(v2 <= vhi){
181  bestPos[0] = 0.5*z1+0.5*zhi;
182  bestPos[1] = 0.5*vlo+0.5*v2;
183  deltaz = zhi-z1;
184  deltav = v2-vlo;
185  }
186  else{
187  // track exits through top of the cell
188  bestPos[0] = 0.5*z1 + 0.5*z2;
189  bestPos[1] = v;
190  deltaz = z2-z1;
191  deltav = vhi-vlo;
192  }
193  }
194  else{
195  // track does not go through cell
196  // check if track is above or below cell
197  if(v1 < vlo){
198  // closest point is bottom right corner of the cell
199  bestPos[0] = zhi;
200  bestPos[1] = vlo;
201  }
202  else{
203  // closest point is top left corner of the cell
204  bestPos[0] = zlo;
205  bestPos[1] = vhi;
206  }
207  }
208  }
209  else{
210  // check if track enters cell from left side
211  if(v1 >= vlo && v1<= vhi){
212  // check if the track exits on the right side of the cell
213  if(v2 >= vlo){
214  bestPos[0] = z;
215  bestPos[1] = 0.5*v1+0.5*v2;
216  deltaz = zhi-zlo;
217  deltav = v2-v1;
218  }
219  else{
220  // track exits through the bottom of the cell
221  bestPos[0] = 0.5*zlo+0.5*z1;
222  bestPos[1] = 0.5*v1+0.5*vlo;
223  deltaz = zlo-z1;
224  deltav = v1-vlo;
225  }
226  }
227  // check if the track enters the cell from the top side
228  else if(v1 >= vhi && v2 <= vhi){
229  // check if the track exits on the right side of the cell
230  if(v2 >= vlo){
231  bestPos[0] = 0.5*z2+0.5*zhi;
232  bestPos[1] = 0.5*vhi+0.5*v2;
233  deltaz = z2-zhi;
234  deltav = vhi-v2;
235  }
236  else{
237  // track exits through the bottom of the cell
238  bestPos[0] = 0.5*z2+0.5*z1;
239  bestPos[1] = v;
240  deltaz = z2-z1;
241  deltav = vhi-vlo;
242  }
243  }
244  else{
245  // track does not go through the cell
246  // check if the track is above or below the cell
247  if(v1 > vhi){
248  // closest point is top right corner of the cell
249  bestPos[0] = zhi;
250  bestPos[1] = vhi;
251  }
252  else{
253  // closest point is the bottom left corner of the cell
254  bestPos[0] = zlo;
255  bestPos[1] = vlo;
256  }
257  // if slope is truly zero than pick midpoint of cell instead of corner points for best z position
258  if(slope == 0){ bestPos[0] = z; }
259  }
260  }
261  } // end of else for if dz == 0
262  return (deltaz*deltaz+deltav*deltav);
263  }
double dz[NP][NC]
z
Definition: test.py:28
double trk::KalmanGeoHelper::CalcMissFrac ( double  dist,
double  w,
geo::View_t  view 
)

Estimation of the mip missed hit fraction as a function of the dist through a cell and cell depth w.

Definition at line 569 of file KalmanGeoHelper.cxx.

References dist, Efficiency(), fclipperLength, flengthmax, fwhigh, fwlow, demo0::length, and w.

Referenced by trk::KalmanTrackMerge::CanJoinTracks(), and trk::KalmanTrack::CheckTrack().

570  {
571  // if distance in cell is less clipper length, clipper hit return 1.0, unreliable hit efficiency
572  if(dist < fclipperLength){ return 1.0; }
573 
574  // define the length in cell to use
575  double length = dist;
576  if(length > flengthmax){ length = flengthmax; }
577  // define what cell depth to use
578  double depth = w;
579  if(depth < fwlow){ depth = fwlow; }
580  if(depth > fwhigh){ depth = fwhigh; }
581 
582  double effhit = Efficiency(depth,length);
583 
584  if(effhit > 1){ effhit = 1.0; }
585  else if(effhit<0){ effhit = 0.0; }
586  return 1.0-effhit;
587 
588  } // end of CalcMissFrac
double Efficiency(double w, double length)
double dist
Definition: runWimpSim.h:113
length
Definition: demo0.py:21
Float_t w
Definition: plot.C:20
int trk::KalmanGeoHelper::CountMissedCellsOnLine ( geo::View_t  view,
double  vi,
double  zi,
double  vf,
double  zf,
int  maxCells 
)

Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm.

Definition at line 267 of file KalmanGeoHelper.cxx.

References getBrightness::cell, geo::PlaneGeo::Cell(), geo::ClampRayToDetectorVolume(), geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geo::GeometryBase::DetLength(), dz, fbadc, fBadChannels, fgeom, fillprivategeo(), geo::CellGeo::HalfW(), makeTrainCVSamples::int, chaninfo::BadChanList::IsBad(), chaninfo::BadChanList::IsLowEfficiency(), geo::kY, std::max(), std::min(), geo::PlaneGeo::Ncells(), nd, ndos, geo::GeometryBase::NPlanes(), std_candles::pl, geo::GeometryBase::Plane(), NDAPDHVSetting::plane, elec2geo::pos, nd_projection_maker::ps, std::swap(), POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by trk::KalmanTrackMerge::CanJoinTracks(), trk::KalmanTrack::CheckTrack(), trk::KalmanTrack::FilterTracker(), and trk::KalmanTrack::SegmentFinder().

268  {
269  static double fDetHalfWidth = fgeom->DetHalfWidth();
270  static double fDetHalfHeight = fgeom->DetHalfHeight();
271  static double fDetLength = fgeom->DetLength();
272  static int fplaneMax = fgeom->NPlanes();
273 
274  fillprivategeo();
275 
276  // if we are in the nd or ndos nead to worry about muon catcher which has different spacing between planes and distance
277  if(nd || ndos){
278  // take the ratio of number of planes to distance from the first 100 planes otherwise will get screwed up by the muon catcher
279  fplaneMax = 100;
280  fDetLength = planez[100];
281  }
282  // Otherwise we would count when the ray was outside the detector in the other view.
283  TVector3 r1(vi,0,zi);
284  TVector3 r2(vf,0,zf);
285  if(view == geo::kY){
286  r1.SetXYZ(0,vi,zi);
287  r2.SetXYZ(0,vf,zf);
288  }
289  geo::ClampRayToDetectorVolume(&r1, &r2, &(*fgeom));
290  int nummiss = 0;
291  if(r1.Z() > r2.Z()) std::swap(r1, r2);
292 
293  double z1 = r1.Z();
294  double z2 = r2.Z();
295 
296  const unsigned int planeMax = fgeom->NPlanes();
297 
298  // find the approximate start plane
299  unsigned int iplane = 0;
300  int pl = (int) z1/(fDetLength/((double)fplaneMax));
301  if(pl < 0){iplane = 0; }
302  else if (pl >= (int)planeMax){iplane = planeMax-1;}
303  else iplane = pl;
304  // check that our initial plane guess is correct
305  // ask if the current position is in this plane
306  if(z1<(planez[iplane]-planehalfd[iplane])){
307  bool foundInitPlane = false;
308  // Did not get the right plane location. Check below if position is lower than the estimated plane center position
309  // only need to worry about geusses that are too high in plane number otherwise accept the estimated guess as the start plane
310  int currTestPlane = iplane-1;
311  while(currTestPlane >= 0 && !foundInitPlane){
312  if(z1>(planez[currTestPlane]-planehalfd[currTestPlane])){
313  iplane = currTestPlane;
314  foundInitPlane = true;
315  }
316  else{
317  if(currTestPlane == 0){iplane = 0;}
318  --currTestPlane;
319  }
320  }
321  }
322 
323  const double v1 = r1[view];
324  const double v2 = r2[view];
325 
326  const double denom = z2-z1;
327  const double dvdz = denom ? (v2-v1)/denom : 0;
328 
329  if(denom == 0 && r2[view] < r1[view]){
330  std::swap(r1, r2);
331  z1 = r1.Z();
332  z2 = r2.Z();
333  }
334 
335  for(unsigned int plane = iplane; plane < planeMax; ++plane){
336  const geo::PlaneGeo* geoplane = fgeom->Plane(plane);
337 
338  // Only look at the planes in the specified view
339  if(geoplane->View() != view) continue;
340  // Assuming here perfect alignment and identical cells across plane
341  const double zc = planez[plane];
342  const double dz = planehalfd[plane];
343 
344  if(z1 > zc+dz) continue;
345  // Plane positions increase monotonically
346  if(z2 < zc-dz) break;
347 
348  // Position entering the cell or starting
349  const double zin = std::max(zc-dz, z1);
350  double vlo = v1+(zin-z1)*dvdz;
351  // Position exiting or ending
352  const double zout = std::min(zc+dz, z2);
353  double vhi = v1+(zout-z1)*dvdz;
354 
355  // if slope is infinite set vlo and vhi as appropriately
356  if(denom == 0){
357  vlo = v1;
358  vhi = v2;
359  }
360 
361  // Ensure they're from low to high instead
362  if(vlo > vhi) std::swap(vlo, vhi);
363 
364  //create an approximate map between position and cell number for cells in the xz view.
365  unsigned int icell = 0;
366  int cell = (int) (vlo+fDetHalfWidth)/(2*fDetHalfWidth/((double) geoplane->Ncells()));
367  int viewind = 0;
368  if(view == geo::kY){
369  cell = (int) (vlo+fDetHalfHeight)/(2*fDetHalfHeight/((double) geoplane->Ncells()));
370  viewind = 1;
371  }
372  if(cell < 0){icell = 0; }
373  else if (cell >= (int) geoplane->Ncells()){icell = geoplane->Ncells()-1;}
374  else icell = cell;
375  // check that our initial cell guess is correct
376  const std::array<double, 3> & pos = cellxyz[plane][icell];
377  double dp = geoplane->Cell(icell)->HalfW();
378  if(vlo<(pos[viewind]-dp)){
379  bool foundInitCell = false;
380  // Did not get the right cell location. Check below if position is lower than the estimated cell center position
381  int currTestCell = icell-1;
382  while(currTestCell >= 0 && !foundInitCell){
383  const std::array<double, 3> & ps = cellxyz[plane][currTestCell];
384  dp = geoplane->Cell(currTestCell)->HalfW();
385  // ask if the current position is in this cell
386  if( vlo>(ps[viewind]-dp)){
387  icell = currTestCell;
388  foundInitCell = true;
389  }
390  else{
391  if(currTestCell == 0){icell = 0;}
392  --currTestCell;
393  }
394  }
395  }
396 
397  for(unsigned int jcell = icell; jcell < geoplane->Ncells();++jcell){
398  const std::array<double, 3> & xyz = cellxyz[plane][jcell];
399  const double vc = xyz[view];
400  const double dv = geoplane->Cell(jcell)->HalfW();
401 
402  if(vlo > vc+dv) continue;
403  // Cell positions increase monotonically
404  if(vhi < vc-dv) break;
405  if(fBadChannels){
406  // check if this is a bad channel
407  if(!fbadc->IsBad(plane,jcell) && !fbadc->IsLowEfficiency(plane,jcell)){ ++nummiss; }
408  }
409  else{ ++nummiss; }
410  if(nummiss > maxCells) return nummiss;
411 
412  } // end for jcell
413  } // end for plane
414  return nummiss;
415  } // end of CountMissedCellsOnLine function
This is the version 2, works for multi-track event.
T max(const caf::Proxy< T > &a, T b)
static std::vector< double > planez
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
static std::vector< double > planehalfd
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double HalfW() const
Definition: CellGeo.cxx:191
double DetLength() const
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
const PlaneGeo * Plane(unsigned int i) const
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
double dz[NP][NC]
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
double DetHalfHeight() const
art::ServiceHandle< chaninfo::BadChanList > fbadc
double DetHalfWidth() const
art::ServiceHandle< geo::Geometry > fgeom
static std::vector< std::vector< std::array< double, 3 > > > cellxyz
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
bool IsBad(int plane, int cell)
bool IsLowEfficiency(int plane, int cell)
int trk::KalmanGeoHelper::CountMissedCellsOnLine ( geo::View_t  view,
double  vi,
double  zi,
double  vf,
double  zf,
std::vector< geo::OfflineChan > &  hitsonline,
std::vector< double > &  distance 
)

Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm.

Definition at line 419 of file KalmanGeoHelper.cxx.

References BestTrackPoint(), getBrightness::cell, geo::PlaneGeo::Cell(), geo::ClampRayToDetectorVolume(), geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geo::GeometryBase::DetLength(), dz, fgeom, fillprivategeo(), geo::CellGeo::HalfW(), makeTrainCVSamples::int, geo::kY, std::max(), std::min(), geo::PlaneGeo::Ncells(), nd, ndos, geo::GeometryBase::NPlanes(), std_candles::pl, geo::GeometryBase::Plane(), NDAPDHVSetting::plane, elec2geo::pos, nd_projection_maker::ps, std::sqrt(), std::swap(), POTSpillRate::view, and geo::PlaneGeo::View().

420  {
421  static double fDetHalfWidth = fgeom->DetHalfWidth();
422  static double fDetHalfHeight = fgeom->DetHalfHeight();
423  static double fDetLength = fgeom->DetLength();
424  static int fplaneMax = fgeom->NPlanes();
425 
426  fillprivategeo();
427 
428  // if we are in the nd or ndos nead to worry about muon catcher which has different spacing between planes and distance
429  if(nd || ndos){
430  // take the ratio of number of planes to distance from the first 100 planes otherwise muon catcher will screw this up
431  fplaneMax = 100;
432  fDetLength = planez[100];
433  }
434  // Otherwise we would count when the ray was outside the detector in the other view.
435  TVector3 r1(vi,0,zi);
436  TVector3 r2(vf,0,zf);
437  if(view == geo::kY){
438  r1.SetXYZ(0,vi,zi);
439  r2.SetXYZ(0,vf,zf);
440  }
441  geo::ClampRayToDetectorVolume(&r1, &r2, &(*fgeom));
442  int nummiss = 0;
443  if(r1.Z() > r2.Z()) std::swap(r1, r2);
444 
445  double z1 = r1.Z();
446  double z2 = r2.Z();
447 
448  static const unsigned int planeMax = fgeom->NPlanes();
449 
450  // find the approximate start plane
451  unsigned int iplane = 0;
452  int pl = (int) z1/(fDetLength/((double)fplaneMax));
453  if(pl < 0){iplane = 0; }
454  else if (pl >= (int)planeMax){iplane = planeMax-1;}
455  else iplane = pl;
456  // check that our initial plane guess is correct
457  // ask if the current position is in this plane
458  if(z1<(planez[iplane]-planehalfd[iplane])){
459  bool foundInitPlane = false;
460  // Did not get the right plane location. Check below if position is lower than the estimated plane center position
461  // only need to worry about geusses that are too high in plane number otherwise accept the estimated guess as the start plane
462  int currTestPlane = iplane-1;
463  while(currTestPlane >= 0 && !foundInitPlane){
464  if(z1>(planez[currTestPlane]-planehalfd[currTestPlane])){
465  iplane = currTestPlane;
466  foundInitPlane = true;
467  }
468  else{
469  if(currTestPlane == 0){iplane = 0;}
470  --currTestPlane;
471  }
472  }
473  }
474 
475  const double v1 = r1[view];
476  const double v2 = r2[view];
477 
478  const double denom = z2-z1;
479  const double dvdz = denom ? (v2-v1)/denom : 0;
480 
481  if(denom == 0 && r2[view] < r1[view]){
482  std::swap(r1, r2);
483  z1 = r1.Z();
484  z2 = r2.Z();
485  }
486 
487  for(unsigned int plane = iplane; plane < planeMax; ++plane){
488  const geo::PlaneGeo* geoplane = fgeom->Plane(plane);
489 
490  // Only look at the planes in the specified view
491  if(geoplane->View() != view) continue;
492  // Assuming here perfect alignment and identical cells across plane
493  const double zc = planez[plane];
494  const double dz = planehalfd[plane];
495 
496  if(z1 > zc+dz) continue;
497  // Plane positions increase monotonically
498  if(z2 < zc-dz) break;
499 
500  // Position entering the cell or starting
501  const double zin = std::max(zc-dz, z1);
502  double vlo = v1+(zin-z1)*dvdz;
503  // Position exiting or ending
504  const double zout = std::min(zc+dz, z2);
505  double vhi = v1+(zout-z1)*dvdz;
506 
507  // if slope is infinite set vlo and vhi as appropriately
508  if(denom == 0){
509  vlo = v1;
510  vhi = v2;
511  }
512 
513  // Ensure they're from low to high instead
514  if(vlo > vhi) std::swap(vlo, vhi);
515 
516  //create an approximate map between position and cell number for cells in the xz view.
517  unsigned int icell = 0;
518  int cell = (int) (vlo+fDetHalfWidth)/(2*fDetHalfWidth/((double) geoplane->Ncells()));
519  int viewind = 0;
520  if(view == geo::kY){
521  cell = (int) (vlo+fDetHalfHeight)/(2*fDetHalfHeight/((double) geoplane->Ncells()));
522  viewind = 1;
523  }
524  if(cell < 0){icell = 0; }
525  else if(cell >= (int) geoplane->Ncells()){icell = geoplane->Ncells()-1;}
526  else{ icell = cell; }
527  // check that our initial cell guess is correct
528  const std::array<double, 3> & pos = cellxyz[plane][icell];
529  double dp = geoplane->Cell(icell)->HalfW();
530  if(vlo<(pos[viewind]-dp)){
531  bool foundInitCell = false;
532  // Did not get the right cell location. Check below if position is lower than the estimated cell center position
533  int currTestCell = icell-1;
534  while(currTestCell >= 0 && !foundInitCell){
535  const std::array<double, 3> & ps = cellxyz[plane][currTestCell];
536  dp = geoplane->Cell(currTestCell)->HalfW();
537  // ask if the current position is in this cell
538  if( vlo>(ps[viewind]-dp)){
539  icell = currTestCell;
540  foundInitCell = true;
541  }
542  else{
543  if(currTestCell == 0){icell = 0;}
544  --currTestCell;
545  }
546  }
547  }
548 
549  for(unsigned int jcell = icell; jcell < geoplane->Ncells();++jcell){
550  const std::array<double, 3> & xyz = cellxyz[plane][jcell];
551  const double vc = xyz[view];
552  const double dv = geoplane->Cell(jcell)->HalfW();
553 
554  if(vlo > vc+dv) continue;
555  // Cell positions increase monotonically
556  if(vhi < vc-dv) break;
557  // get the distance in the cell
558  std::array<double, 2> bestPos;
559  double dist2 = BestTrackPoint(xyz[2],vc,dz,dv,bestPos,z1,v1,z2-z1,v2-v1);
560  ++nummiss;
561  hitsonline.push_back(geo::OfflineChan(plane,jcell));
562  distance.push_back(sqrt(dist2));
563  } // end for jcell
564  } // end for plane
565  return nummiss;
566  } // end of CountMissedCellsOnLine function
This is the version 2, works for multi-track event.
T max(const caf::Proxy< T > &a, T b)
static std::vector< double > planez
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
static std::vector< double > planehalfd
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double HalfW() const
Definition: CellGeo.cxx:191
double DetLength() const
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
const PlaneGeo * Plane(unsigned int i) const
unsigned distance(const T &t1, const T &t2)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
double dz[NP][NC]
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
double DetHalfHeight() const
double DetHalfWidth() const
A (plane, cell) pair.
Definition: OfflineChan.h:17
art::ServiceHandle< geo::Geometry > fgeom
static std::vector< std::vector< std::array< double, 3 > > > cellxyz
double BestTrackPoint(double z, double v, double halfz, double halfv, std::array< double, 2 > &bestPos, double z0, double v0, double dz, double dv)
Estimate the position on a line within a box.
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
double trk::KalmanGeoHelper::Efficiency ( double  w,
double  length 
)
private

Definition at line 592 of file KalmanGeoHelper.cxx.

References nd_projection_maker::eff, effpar0, effpar1, effpar2, effpar3, fd, nd, and tb.

Referenced by CalcMissFrac().

593  {
594  double eff = 0.0;
595 
596  if(fd){
597  double llow = 3.4;
598  double lhigh = 3.75;
599 
600  if( length < llow){
601  eff = effpar0.Eval(w)*length + effpar1.Eval(w);
602  }
603  else if( length < lhigh){
604  eff = effpar2.Eval(w)*length + (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
605  }
606  else{
607  eff = effpar3.Eval(w)*length + (effpar2.Eval(w)-effpar3.Eval(w))*lhigh +
608  (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
609  }
610  }
611  else if(nd){
612  double llow = 3.0;
613  double lhigh = 3.75;
614 
615  if( length < llow){
616  eff = effpar0.Eval(w)*length + effpar1.Eval(w);
617  }
618  else if( length < lhigh){
619  eff = effpar2.Eval(w)*length + (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
620  }
621  else{
622  eff = effpar3.Eval(w)*length + (effpar2.Eval(w)-effpar3.Eval(w))*lhigh +
623  (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
624  }
625  }
626  else if(tb){
627  double llow = 3.0;
628  double lhigh = 3.75;
629 
630  if( length < llow){
631  eff = effpar0.Eval(w)*length + effpar1.Eval(w);
632  }
633  else if( length < lhigh){
634  eff = effpar2.Eval(w)*length + (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
635  }
636  else{
637  eff = effpar3.Eval(w)*length + (effpar2.Eval(w)-effpar3.Eval(w))*lhigh +
638  (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
639  }
640  }
641 
642  else{ eff = 1.0; }
643  // unless clipper, ndos shouldn't miss hits (need to check for bad channel definition of efficiency/badchannels outside of this function)
644 
645  return eff;
646  }
length
Definition: demo0.py:21
Float_t w
Definition: plot.C:20
void trk::KalmanGeoHelper::fillprivategeo ( )
private

Definition at line 21 of file KalmanGeoHelper.cxx.

References plot_validation_datamc::c, geo::PlaneGeo::Cell(), fgeom, geo::CellGeo::GetCenter(), geo::CellGeo::HalfD(), MECModelEnuComparisons::i, geo::PlaneGeo::Ncells(), geo::GeometryBase::NPlanes(), and geo::GeometryBase::Plane().

Referenced by CountMissedCellsOnLine(), and MatchTrajectoryToPlaneInView().

22  {
23  if(!planez.empty()) return;
24 
25  double posp[3];
26  for(unsigned int i = 0; i < fgeom->NPlanes(); i++){
27  fgeom->Plane(i)->Cell(0)->GetCenter(posp);
28  planez.push_back(posp[2]);
29  planehalfd.push_back(fgeom->Plane(i)->Cell(0)->HalfD());
30  }
31 
32  cellxyz.resize(fgeom->NPlanes());
33 
34  for(unsigned int p = 0; p < fgeom->NPlanes(); p++){
35  const geo::PlaneGeo* geoplane = fgeom->Plane(p);
36  for(unsigned int c = 0; c < geoplane->Ncells(); c++){
37  double xyz[3];
38  fgeom->Plane(p)->Cell(c)->GetCenter(xyz);
39  std::array<double, 3> vxyz;
40  vxyz[0] = xyz[0];
41  vxyz[1] = xyz[1];
42  vxyz[2] = xyz[2];
43  cellxyz[p].push_back(vxyz);
44  }
45  }
46 
47  }
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
static std::vector< double > planez
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
static std::vector< double > planehalfd
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const PlaneGeo * Plane(unsigned int i) const
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
art::ServiceHandle< geo::Geometry > fgeom
static std::vector< std::vector< std::array< double, 3 > > > cellxyz
unsigned int NPlanes() const
int trk::KalmanGeoHelper::MatchTrajectoryToPlane ( TVector3  traj,
int  minplane,
int  maxplane 
)

Finds the plane best matched to a trajectory point This method takes a trajectory point and finds the plane that closest matches this position in either view between minplane and maxplane.

Parameters
trajTrajectory point to be matched to a plane
minplaneMinimum plane to search for a match to the trajectory point
maxplaneMaximum plane to search for a match to the trajectory point
Returns
Plane that matches the trajectory point between minplane and maxplane

Definition at line 650 of file KalmanGeoHelper.cxx.

References geo::PlaneGeo::Cell(), geo::GeometryBase::CellId(), PandAna.Demos.tute_pid_validation::cid, geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), fgeom, fnmaxplane, geo::CellGeo::HalfD(), geo::CellGeo::HalfW(), geo::GeometryBase::IdToPlane(), makeTrainCVSamples::int, ip, geo::kY, geo::PlaneGeo::Ncells(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, registry_explorer::v, and geo::PlaneGeo::View().

Referenced by MatchTrajectoryToPlaneInView(), and trk::KalmanTrackMerge::ViewMergeTracks().

651  {
652  int plane = minplane;
653  int mintrialplane = minplane;
654  int maxtrialplane = maxplane;
655  if(minplane < 0){ mintrialplane = 0; }
656  if(maxplane > fnmaxplane-1){ maxtrialplane = fnmaxplane-1; }
657 
658  // Try finding what plane/cell this point is in
659  const geo::CellUniqueId cid = fgeom->CellId(traj.X(),traj.Y(),traj.Z(),1,1,0.1,2);
660 
661  // if found the cell unique id use it; otherwise backout plane from individual plane width checks
662  if(cid){ fgeom->IdToPlane(cid, &plane); }
663  else{
664  for(int ip = maxtrialplane; ip >= mintrialplane; --ip){
665  int trialPlane = ip;
666 
667  const geo::PlaneGeo* geoplane = fgeom->Plane(trialPlane);
668  // get the estimated cell in this plane
669  int mintrialcell = 0;
670  int maxtrialcell = geoplane->Ncells();
671  double v = traj.X();
672  int trialCell = (int) (v+fgeom->DetHalfWidth())/(2*fgeom->DetHalfWidth()/((double) (geoplane->Ncells()-1)));
673  int viewind = 0;
674  if(geoplane->View() == geo::kY){
675  v = traj.Y();
676  trialCell = (int) (v+fgeom->DetHalfHeight())/(2*fgeom->DetHalfHeight()/((double) (geoplane->Ncells()-1)));
677  viewind = 1;
678  }
679 
680  if(trialCell < 0){trialCell = 0; }
681  else if (trialCell >= (int) geoplane->Ncells()){trialCell = geoplane->Ncells()-1;}
682 
683  // check that our initial cell guess is correct
684  const std::array<double, 3> & trialCenter = cellxyz[ip][trialCell];
685  double dp = geoplane->Cell(trialCell)->HalfW();
686 
687  if(v < (trialCenter[viewind]-dp)){
688  // trial cell is too high, try iterating down to cell 0 until find right cell
689  maxtrialcell = trialCell;
690  for(int ic = maxtrialcell; ic >= mintrialcell; --ic){
691  const std::array<double, 3> & tc = cellxyz[trialPlane][ic];
692  dp = geoplane->Cell(ic)->HalfW();
693  trialCell = ic;
694  if(v > (tc[viewind]-dp)){ break; }
695  }
696  }
697  else if(v > (trialCenter[viewind]+dp)){
698  // trial cell is too low, try iterating up to Ncells in plane until find right cell
699  mintrialcell = trialCell;
700  for(int ic = mintrialcell; ic <= maxtrialcell; ++ic){
701  const std::array<double, 3> & tc = cellxyz[trialPlane][ic];
702  dp = geoplane->Cell(ic)->HalfW();
703  trialCell = ic;
704  if(v > (tc[viewind]-dp)){ break; }
705  }
706  }
707 
708  const std::array<double, 3> & tc = cellxyz[trialPlane][trialCell];
709  double trialPlaneWidth = fgeom->Plane(trialPlane)->Cell(trialCell)->HalfD();
710  // starting from max plane first so if the z is larger than the smallest z
711  // on the plane then this is the correct plane
712  if(traj.Z() > (tc[2]-trialPlaneWidth)){
713  plane = trialPlane;
714  break;
715  }
716  }
717  }
718 
719  return plane;
720  } // end of MatchTrajectoryToPlane
double HalfD() const
Definition: CellGeo.cxx:205
const PlaneGeo * IdToPlane(const CellUniqueId &id, int *iplane) const
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double HalfW() const
Definition: CellGeo.cxx:191
TString ip
Definition: loadincs.C:5
const PlaneGeo * Plane(unsigned int i) const
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
double DetHalfHeight() const
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
double DetHalfWidth() const
art::ServiceHandle< geo::Geometry > fgeom
static std::vector< std::vector< std::array< double, 3 > > > cellxyz
int trk::KalmanGeoHelper::MatchTrajectoryToPlaneInView ( TVector3  traj,
geo::View_t  view,
int  minplane,
int  maxplane 
)

Variation of MatchTrajectoryToPlane that limits the matching plane to be in the input view.

Definition at line 723 of file KalmanGeoHelper.cxx.

References ana::assert(), stan::math::fabs(), fgeom, fillprivategeo(), geo::kPLANE_NOT_FOUND, geo::kXorY, MatchTrajectoryToPlane(), geo::GeometryBase::NextPlaneOtherView(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by trk::KalmanTrackMerge::CalcMatchScore(), and trk::KalmanTrackMerge::JoinTracks().

724  {
725 
726  assert(fgeom->Plane(minplane)->View() == view);
727  assert(fgeom->Plane(maxplane)->View() == view);
728  int plane = MatchTrajectoryToPlane(traj,minplane,maxplane);
729  geo::View_t trialview = fgeom->Plane(plane)->View();
730 
731  fillprivategeo();
732 
733  if(view == geo::kXorY || trialview == view){ return plane; }
734  else{
735  int lowtrialplane = fgeom->NextPlaneOtherView(plane,-1);
736  int hightrialplane = fgeom->NextPlaneOtherView(plane,+1);
737  if(lowtrialplane == geo::kPLANE_NOT_FOUND){
738  if(hightrialplane == geo::kPLANE_NOT_FOUND){ return minplane; }
739  else{ return hightrialplane; }
740  }
741  else if(hightrialplane == geo::kPLANE_NOT_FOUND){
742  if(lowtrialplane >= minplane){ return lowtrialplane; }
743  else{ return minplane; }
744  }
745  else{
746  // check whether the high or low plane is the closest in z
747  if(fabs(planez[lowtrialplane]+planehalfd[lowtrialplane] - traj.Z()) <=
748  fabs(planez[hightrialplane]-planehalfd[hightrialplane] - traj.Z())){ return lowtrialplane; }
749  else{ return hightrialplane; }
750  }
751  }
752  return minplane;
753 
754  } // end of MatchTrajectoryToPlaneInView
static std::vector< double > planez
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
static std::vector< double > planehalfd
const PlaneGeo * Plane(unsigned int i) const
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
art::ServiceHandle< geo::Geometry > fgeom
assert(nhit_max >=nhit_nbins)
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const
int MatchTrajectoryToPlane(TVector3 traj, int minplane, int maxplane)
Finds the plane best matched to a trajectory point This method takes a trajectory point and finds the...
void trk::KalmanGeoHelper::SetupHelper ( )
private

Definition at line 62 of file KalmanGeoHelper.cxx.

References om::cerr, geo::GeometryBase::DetId(), E, effpar0, effpar1, effpar2, effpar3, allTimeWatchdog::endl, fclipperLength, fd, fgeom, flengthmax, fnmaxplane, fwhigh, fwlow, novadaq::cnv::kFARDET, novadaq::cnv::kNDOS, novadaq::cnv::kNEARDET, novadaq::cnv::kTESTBEAM, nd, ndos, geo::GeometryBase::NPlanes(), and tb.

Referenced by KalmanGeoHelper().

63  {
65  if(fgeom->DetId() == novadaq::cnv::kFARDET) fd = true;
66  else if(fgeom->DetId() == novadaq::cnv::kNEARDET) nd = true;
67  else if(fgeom->DetId() == novadaq::cnv::kNDOS) ndos = true;
68  else if(fgeom->DetId() == novadaq::cnv::kTESTBEAM) tb = true;
69  else{
70  std::cerr<<"This geometry is not ndos, tb, nd, or fd. This object doesn't know what to do"<<std::endl;
71  }
72  if(fd){
73  fclipperLength = 1.5;
74  fwlow = -700;
75  fwhigh = 700;
76  flengthmax = 6.0;
77  effpar0 = TF1("par0","pol4",fwlow,fwhigh);
78  effpar1 = TF1("par1","pol4",fwlow,fwhigh);
79  effpar2 = TF1("par2","pol5",fwlow,fwhigh);
80  effpar3 = TF1("par3","pol4",fwlow,fwhigh);
81 
82  // low threshold tuning
83  effpar0.SetParameters(0.204074, -2.45217E-6, -3.61656E-7, 4.38289E-11, 4.3173E-13);
84  effpar1.SetParameters(0.126228, 3.77528E-4, 7.85744E-7, -4.65944E-10, -1.0345E-12);
85  effpar2.SetParameters(0.327222, -5.80372E-4, 7.42265E-7, 1.27296E-9, -1.13353E-12, -1.13839E-15);
86  effpar3.SetParameters(1.35419E-2, -8.28103E-5, 8.92098E-8, 5.64673E-11, -6.05079E-14);
87 
88  }
89  else if(nd){
90  fclipperLength = 1.5;
91  fwlow = -200;
92  fwhigh = 200;
93  flengthmax = 6.0;
94  effpar0 = TF1("par0","pol4",fwlow,fwhigh);
95  effpar1 = TF1("par1","pol4",fwlow,fwhigh);
96  effpar2 = TF1("par2","pol5",fwlow,fwhigh);
97  effpar3 = TF1("par3","pol4",fwlow,fwhigh);
98 
99  // low threshold tuning
100  effpar0.SetParameters(0.139861, -5.05093E-5, 2.66785E-7, -3.65295E-10, -1.46979E-11);
101  effpar1.SetParameters(0.525903, 2.00583E-4, -4.07681E-7, 5.28495E-10, 2.01945E-11);
102  effpar2.SetParameters(0.0239177, -6.04344E-5, -1.08969E-7, -3.69116E-10, 2.01116E-11, 5.79678E-14);
103  effpar3.SetParameters(-1.03056E-3, 1.7984E-5, -3.04568E-7, -7.85076E-10, 3.55217E-12);
104 
105  }
106  else if(ndos){
107  fclipperLength = 1.5;
108  }
109  // not sure how the effpars are defined so for test beam just use nd setup
110  else if(tb){
111  fclipperLength = 1.5;
112  fwlow = -130;
113  fwhigh = 130;
114  flengthmax = 6.0;
115  effpar0 = TF1("par0","pol4",fwlow,fwhigh);
116  effpar1 = TF1("par1","pol4",fwlow,fwhigh);
117  effpar2 = TF1("par2","pol5",fwlow,fwhigh);
118  effpar3 = TF1("par3","pol4",fwlow,fwhigh);
119 
120  // low threshold tuning this produces a dist with relatively flat efficiencies as a fn of W
121  effpar0.SetParameters(0.139861, -5.05093E-5, 2.66785E-7, -3.65295E-10, -1.46979E-11);
122  effpar1.SetParameters(0.525903, 2.00583E-4, -4.07681E-7, 5.28495E-10, 2.01945E-11);
123  effpar2.SetParameters(0.0239177, -6.04344E-5, -1.08969E-7, -3.69116E-10, 2.01116E-11, 5.79678E-14);
124  effpar3.SetParameters(-1.03056E-3, 1.7984E-5, -3.04568E-7, -7.85076E-10, 3.55217E-12);
125  }
126 
127  }
OStream cerr
Definition: OStream.cxx:7
Far Detector at Ash River, MN.
Float_t E
Definition: plot.C:20
Prototype Near Detector on the surface at FNAL.
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
art::ServiceHandle< geo::Geometry > fgeom
unsigned int NPlanes() const

Member Data Documentation

TF1 trk::KalmanGeoHelper::effpar0
private

Definition at line 71 of file KalmanGeoHelper.h.

Referenced by Efficiency(), and SetupHelper().

TF1 trk::KalmanGeoHelper::effpar1
private

Definition at line 72 of file KalmanGeoHelper.h.

Referenced by Efficiency(), and SetupHelper().

TF1 trk::KalmanGeoHelper::effpar2
private

Definition at line 73 of file KalmanGeoHelper.h.

Referenced by Efficiency(), and SetupHelper().

TF1 trk::KalmanGeoHelper::effpar3
private

Definition at line 74 of file KalmanGeoHelper.h.

Referenced by Efficiency(), and SetupHelper().

art::ServiceHandle<chaninfo::BadChanList> trk::KalmanGeoHelper::fbadc
private

Definition at line 61 of file KalmanGeoHelper.h.

Referenced by CountMissedCellsOnLine().

bool trk::KalmanGeoHelper::fBadChannels
private

Definition at line 59 of file KalmanGeoHelper.h.

Referenced by CountMissedCellsOnLine(), and KalmanGeoHelper().

double trk::KalmanGeoHelper::fclipperLength
private

Definition at line 67 of file KalmanGeoHelper.h.

Referenced by CalcMissFrac(), and SetupHelper().

bool trk::KalmanGeoHelper::fd = false
private

Definition at line 64 of file KalmanGeoHelper.h.

Referenced by Efficiency(), and SetupHelper().

art::ServiceHandle<geo::Geometry> trk::KalmanGeoHelper::fgeom
private
double trk::KalmanGeoHelper::flengthmax
private

Definition at line 70 of file KalmanGeoHelper.h.

Referenced by CalcMissFrac(), and SetupHelper().

int trk::KalmanGeoHelper::fnmaxplane
private

Definition at line 66 of file KalmanGeoHelper.h.

Referenced by MatchTrajectoryToPlane(), and SetupHelper().

double trk::KalmanGeoHelper::fwhigh
private

Definition at line 69 of file KalmanGeoHelper.h.

Referenced by CalcMissFrac(), and SetupHelper().

double trk::KalmanGeoHelper::fwlow
private

Definition at line 68 of file KalmanGeoHelper.h.

Referenced by CalcMissFrac(), and SetupHelper().

bool trk::KalmanGeoHelper::nd = false
private

Definition at line 63 of file KalmanGeoHelper.h.

Referenced by CountMissedCellsOnLine(), Efficiency(), and SetupHelper().

bool trk::KalmanGeoHelper::ndos = false
private

Definition at line 62 of file KalmanGeoHelper.h.

Referenced by CountMissedCellsOnLine(), and SetupHelper().

bool trk::KalmanGeoHelper::tb = false
private

Definition at line 65 of file KalmanGeoHelper.h.

Referenced by Efficiency(), and SetupHelper().


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