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

An algorithm to perform cosmic ray track fitting. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-21/TrackFit/WindowTrackingAlg.h"

Inheritance diagram for trk::WindowTrackingAlg:
trk::TrackAlg

Public Member Functions

 WindowTrackingAlg (fhicl::ParameterSet const &pset)
 
virtual ~WindowTrackingAlg ()
 
std::vector< rb::TrackMakeTrack (const rb::Cluster *slice)
 

Private Member Functions

rb::Track Make3DTrack (std::vector< art::Ptr< rb::CellHit > > &allHits, std::vector< art::Ptr< rb::CellHit > > &removedHits, double &chiSqr, bool(*comp)(unsigned short int, unsigned short int))
 
rb::Track FitView (std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &hits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, std::vector< art::Ptr< rb::CellHit > > &removedHits, geo::View_t view, double &chiSqr)
 
rb::Track MakeViewTrack (std::vector< art::Ptr< rb::CellHit > > const &trkHits, geo::View_t const &view, double &chiSqr, bool const &dir, std::vector< art::Ptr< rb::CellHit > > &removedHits)
 
void LookForBremsstrahlungHits (rb::Track &track, std::vector< art::Ptr< rb::CellHit > > &removedHits)
 
unsigned int FitWindow (std::vector< art::Ptr< rb::CellHit > > const &prevWindow, std::vector< art::Ptr< rb::CellHit > > &currWindow, std::vector< art::Ptr< rb::CellHit > > &removedHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, geo::View_t const &view)
 
rb::Track ShortTrack (std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &xViewHits, std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &yViewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, double &chiSqr)
 
rb::Track ShortViewTrack (std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, bool const &twoPlanes, int const &firstPlane, int const &secondPlane, geo::View_t const &view)
 
void ShortTrackExtraPlane (std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, std::vector< std::pair< geo::View_t, TVector3 > > &trajPoints, rb::Track &viewTrack, unsigned short int const &removedPlane)
 
void ShortTrackRemoveExtraPlane (std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > &planeCoG, unsigned short int &removedPlane)
 
void SetTrackEndPoints (rb::Track &track)
 
void FindEndPoint (geo::View_t const &view, double *cellCenter, TVector3 const &endPoint, TVector3 const &endDir, bool const &start, rb::Track &track)
 
void CheckTrackDirectionInY (rb::Track &track, TVector3 &start, TVector3 &end)
 
void DetermineInitialDirection (rb::Track &track)
 
bool ValidTrack (std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &planeToHit)
 

Private Attributes

double fCellHalfDepth
 half of cell size in z More...
 
double fDHitGood
 
unsigned int fWindowSize
 Number of planes in the sliding window. More...
 
unsigned int fMinViewPlanes
 Minimum number of planes in a view to fit. More...
 
unsigned int fMinHitsTryAgain
 Minimum number of hits left to look for another track. More...
 
double fMinHitFraction
 Minimum value of hits/plane to look for a track. More...
 
unsigned int fNTrajectory
 number of trajectory points to use in determining incoming direction More...
 
int fMaxPlaneSeparation
 Maximum number of planes allowed to separate planes with hits in slice. More...
 
art::ServiceHandle< geo::GeometryfGeo
 handle to the geometry More...
 

Detailed Description

An algorithm to perform cosmic ray track fitting.

Definition at line 32 of file WindowTrackingAlg.h.

Constructor & Destructor Documentation

trk::WindowTrackingAlg::WindowTrackingAlg ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 55 of file WindowTrackingAlg.cxx.

References fDHitGood, fMaxPlaneSeparation, fMinHitFraction, fMinHitsTryAgain, fMinViewPlanes, fNTrajectory, fWindowSize, and fhicl::ParameterSet::get().

56  : TrackAlg(pset)
57  , fCellHalfDepth(-1.)
58  {
59  fDHitGood = pset.get<double >("DHitGood" );
60  fMinViewPlanes = pset.get<unsigned int>("MinViewPlanes" );
61  fWindowSize = pset.get<unsigned int>("WindowSize" );
62  fMinHitsTryAgain = pset.get<unsigned int>("MinHitsToTryAgain" );
63  fMinHitFraction = pset.get<double >("MinHitFraction" );
64  fNTrajectory = pset.get<unsigned int>("NTrajectoryForDir" );
65  fMaxPlaneSeparation = pset.get<int >("MaxPlaneSeparation");
66  }
unsigned int fMinHitsTryAgain
Minimum number of hits left to look for another track.
double fCellHalfDepth
half of cell size in z
unsigned int fMinViewPlanes
Minimum number of planes in a view to fit.
double fMinHitFraction
Minimum value of hits/plane to look for a track.
TrackAlg(fhicl::ParameterSet const &pset)
Definition: TrackAlg.cxx:15
unsigned int fWindowSize
Number of planes in the sliding window.
unsigned int fNTrajectory
number of trajectory points to use in determining incoming direction
int fMaxPlaneSeparation
Maximum number of planes allowed to separate planes with hits in slice.
trk::WindowTrackingAlg::~WindowTrackingAlg ( )
virtual

Definition at line 70 of file WindowTrackingAlg.cxx.

70 {}

Member Function Documentation

void trk::WindowTrackingAlg::CheckTrackDirectionInY ( rb::Track track,
TVector3 &  start,
TVector3 &  end 
)
private

Definition at line 1341 of file WindowTrackingAlg.cxx.

References rb::Cluster::AllCells(), test3::ctrack, febshutoff_auto::end, rb::Track::NTrajectoryPoints(), and rb::Track::TrajectoryPoint().

Referenced by Make3DTrack(), and ShortTrack().

1344  {
1345 
1346  if(start.Y() < end.Y()){
1347 
1348  TVector3 pt2 = track.TrajectoryPoint(track.NTrajectoryPoints()-2);
1349  rb::Track ctrack(track.AllCells(), end, pt2-end);
1350  for(int p = track.NTrajectoryPoints()-2; p >= 0; --p)
1351  ctrack.AppendTrajectoryPoint(track.TrajectoryPoint(p));
1352 
1353  track = ctrack;
1354 
1355  start = ctrack.Start();
1356  end = ctrack.Stop();
1357  }// end if we have to flip the direction on this track
1358 
1359  return;
1360  }
size_t NTrajectoryPoints() const
Definition: Track.h:83
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
ctrack
Definition: test3.py:11
const char * p
Definition: xmltok.h:285
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
void trk::WindowTrackingAlg::DetermineInitialDirection ( rb::Track track)
private

Definition at line 988 of file WindowTrackingAlg.cxx.

References febshutoff_auto::end, fNTrajectory, geo::LinFitMinDperp(), rb::Track::NTrajectoryPoints(), elec2geo::pos, rb::Prong::SetDir(), febshutoff_auto::start, std::swap(), confusionMatrixTree::t, rb::Track::TrajectoryPoint(), w, submit_syst::x, submit_syst::y, and test::z.

Referenced by Make3DTrack().

989  {
990 
991  // use at most the first N (default = 20) trajectory points on a track to get
992  // the incoming direction
993  size_t numTrajPts = (track.NTrajectoryPoints() > fNTrajectory) ? fNTrajectory : track.NTrajectoryPoints();
994 
995  // put the first n trajectory points into vectors
996  std::vector<double> x(numTrajPts, 0.);
997  std::vector<double> y(numTrajPts, 0.);
998  std::vector<double> z(numTrajPts, 0.);
999  std::vector<double> w(numTrajPts, 1.);
1000 
1001  TVector3 pos;
1002  for(size_t t = 0; t < numTrajPts; ++t){
1003  pos = track.TrajectoryPoint(t);
1004  x[t] = pos.X();
1005  y[t] = pos.Y();
1006  z[t] = pos.Z();
1007  }
1008 
1009  double xzStart[2] = {0.};
1010  double xzEnd[2] = {0.};
1011  double yzStart[2] = {0.};
1012  double yzEnd[2] = {0.};
1013 
1014  geo::LinFitMinDperp(z, x, w, &xzStart[0], &xzStart[1], &xzEnd[0], &xzEnd[1]);
1015  geo::LinFitMinDperp(z, y, w, &yzStart[0], &yzStart[1], &yzEnd[0], &yzEnd[1]);
1016 
1017  if(xzStart[0] < z.front()) std::swap(xzStart, xzEnd);
1018  if(yzStart[0] < z.front()) std::swap(yzStart, yzEnd);
1019 
1020  TVector3 start(xzStart[1], yzStart[1], z.front());
1021  TVector3 end (xzEnd[1], yzEnd[1], z.back());
1022 
1023  track.SetDir(end - start);
1024 
1025  return;
1026  }
size_t NTrajectoryPoints() const
Definition: Track.h:83
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
z
Definition: test.py:28
unsigned int fNTrajectory
number of trajectory points to use in determining incoming direction
Float_t w
Definition: plot.C:20
void trk::WindowTrackingAlg::FindEndPoint ( geo::View_t const &  view,
double *  cellCenter,
TVector3 const &  endPoint,
TVector3 const &  endDir,
bool const &  start,
rb::Track track 
)
private

Definition at line 1449 of file WindowTrackingAlg.cxx.

References std::abs(), rb::Track::AppendTrajectoryPoint(), geo::PlaneGeo::Cell(), fGeo, geo::CellGeo::HalfD(), std::isnan(), geo::kY, geo::GeometryBase::Plane(), rb::Track::PrependTrajectoryPoint(), gen_hdf5record::pt, febshutoff_auto::start, and registry_explorer::v.

Referenced by SetTrackEndPoints().

1455  {
1456  double deltaZ = fGeo->Plane(0)->Cell(0)->HalfD();
1457 
1458  double &v = cellCenter[0];
1459  double epV = endPoint.X();
1460  double epVDir = epDir.X();
1461  int signDirV = 1.;
1462 
1463 // std::string startEnd("start");
1464 // if(!start) startEnd = "end";
1465 
1466  if(view == geo::kY){
1467  v = cellCenter[1];
1468  epV = endPoint.Y();
1469  epVDir = epDir.Y();
1470  }
1471 
1472 // LOG_VERBATIM("WindowTrackingAlg") << startEnd << " point: "
1473 // << endPoint.X() << " " << endPoint.Y() << " " << endPoint.Z()
1474 // << " " << epDir.X() << " " << epDir.Y() << " " << epDir.Z()
1475 // << " cell at "
1476 // << v << " " << cellCenter[2] << " " << deltaZ;
1477 
1478  // bail if the direction in the view is 0
1479  // the current end point will have to be good enough
1480  if(epVDir == 0.) return;
1481 
1482  double deltaV = std::abs(v - epV);
1483  TVector3 dV(epDir.X()/epVDir, epDir.Y()/epVDir, epDir.Z()/epVDir);
1484  // point the track backwards
1485  if (epVDir > 0. && start) signDirV = -1;
1486  else if(epVDir < 0. && !start) signDirV = -1;
1487  TVector3 pt(endPoint.X() + deltaV*signDirV*dV.X(),
1488  endPoint.Y() + deltaV*signDirV*dV.Y(),
1489  endPoint.Z() + deltaV*signDirV*dV.Z());
1490 
1491  // only prepend if the point is different from the original starting point
1492  if(pt != endPoint){
1493 
1494  // if the new start point is more than half a cell width
1495  // from the z position of the first cell in the start plane
1496  // just carry the track to the boundary of the plane
1497  if(std::abs(pt.Z() - cellCenter[2]) > deltaZ){
1498  signDirV = 1.;
1499  if (epDir.Z() > 0. && start) signDirV = -1.;
1500  else if(epDir.Z() < 0. && !start) signDirV = -1.;
1501  pt.SetXYZ(endPoint.X() + deltaZ*signDirV*epDir.X()/epDir.Z(),
1502  endPoint.Y() + deltaZ*signDirV*epDir.Y()/epDir.Z(),
1503  endPoint.Z() + deltaZ*signDirV);
1504  }
1505 
1506  // check that we don't have nans for the new coordinates
1507  if( !std::isnan(pt.X()) &&
1508  !std::isnan(pt.Y()) &&
1509  !std::isnan(pt.Z())){
1510  if(start) track.PrependTrajectoryPoint(pt);
1511  else track.AppendTrajectoryPoint(pt);
1512 
1513  // LOG_VERBATIM("WindowTrackingAlg") << "new " << startEnd << " point is: "
1514  // << pt.X() << " "
1515  // << pt.Y() << " "
1516  // << pt.Z();
1517 
1518  }
1519 
1520  } // end if the new point and the previous endPoint are different
1521 
1522  return;
1523  }
double HalfD() const
Definition: CellGeo.cxx:205
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const PlaneGeo * Plane(unsigned int i) const
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
void PrependTrajectoryPoint(TVector3 pt)
Support constructing tracks backwards.
Definition: Track.cxx:134
rb::Track trk::WindowTrackingAlg::FitView ( std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &  hits,
std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &  planeCoG,
std::vector< art::Ptr< rb::CellHit > > &  removedHits,
geo::View_t  view,
double &  chiSqr 
)
private

Definition at line 459 of file WindowTrackingAlg.cxx.

References dir, FitWindow(), fWindowSize, make_syst_table_plots::h, hits(), and MakeViewTrack().

Referenced by Make3DTrack().

466  {
467  // start off with an N plane window, then step out one plane at a time to add to the window
468  size_t numPlanes = hits.size();
469 
470  // another vector that has all the hits for a given window
471  std::vector< art::Ptr<rb::CellHit> > currWindow;
472  std::vector< art::Ptr<rb::CellHit> > prevWindow;
473 
474  // get the iterator to the hits of the first plane in the window
475  auto startItr = hits.begin();
476 
477  // get the iterator the hits of the last plane in the window
478  auto endItr = hits.begin();
479 
480  if(numPlanes > fWindowSize) std::advance(endItr, fWindowSize);
481  else endItr = hits.end();
482 
483  bool dir = (startItr->first > endItr->first);
484 
485  // fit the initial window
486  unsigned short int startPlane = startItr->second.front()->Plane();
487  for(auto itr = startItr; itr != endItr; ++itr){
488  for(size_t h = 0; h < itr->second.size(); ++h){
489  currWindow.push_back(itr->second[h]);
490  // LOG_VERBATIM("WindowTrackAlg") << "adding to initial window: " << itr->second[h]->Plane() << " " << itr->second[h]->Cell();
491  }
492  }
493 
494  // this is the fit of the initial window, there is nothing in prevWindow
495  // if all the hits are removed from currWindow, return an empty track
496  this->FitWindow(prevWindow, currWindow, removedHits, planeCoG, view);
497  if(currWindow.size() < 1){
498  //LOG_DEBUG("WindowTrack") << "currWindow size is " << currWindow.size()
499  // << " returning empty track";
500  return rb::Track();
501  }
502 
503  // put the current window's hits into the vector we'll use
504  // to make the cluster
505  std::vector< art::Ptr<rb::CellHit> > clHits(currWindow);
506 
507  // if there are no more planes to look at, just make the track
508  // otherwise, endItr is already on the element past the initial
509  // window
510  if(endItr == hits.end() ){
511  // make the track from the window
512  return this->MakeViewTrack(clHits, view, chiSqr, dir, removedHits);
513  }
514 
515  prevWindow.swap(currWindow);
516 
517  //LOG_DEBUG("WindowTrackingAlg") << "start sliding window";
518 
519  // slide the window until you run out of planes
520  while(endItr != hits.end()){
521 
522  //LOG_DEBUG("WindowTrackingAlg") << "prev window size now " << prevWindow.size()
523  // << " curr window has size " << endItr->second.size()
524  // << " starts with plane " << endItr->second.front()->Plane();
525 
526  // remove the first plane of hits from the previousWindow, but only if
527  // we added hits from the current window the last time
528  auto itr = prevWindow.begin();
529  while( itr != prevWindow.end() && currWindow.size() > 0){
530  //LOG_DEBUG("WindowTrackingAlg") << (*itr)->Plane() << " " << (*itr)->Cell();
531  if((*itr)->Plane() == startPlane) itr = prevWindow.erase(itr);
532  else ++itr;
533  }
534  startPlane = prevWindow.front()->Plane();
535 
536  // put the current plane's hits to the the currWindow vector
537  currWindow.swap(endItr->second);
538 
539  this->FitWindow(prevWindow, currWindow, removedHits, planeCoG, view);
540 
541  // add the surviving currWindow hits to the prevWindow and to the clHits
542  for(size_t h = 0; h < currWindow.size(); ++h){
543  prevWindow.push_back(currWindow[h]);
544  clHits.push_back(currWindow[h]);
545  }
546 
547  ++endItr;
548  }// end loop over the windows
549 
550  // make the track using the hits that were found along the
551  // trajectory
552  return this->MakeViewTrack(clHits, view, chiSqr, dir, removedHits);
553  }
rb::Track MakeViewTrack(std::vector< art::Ptr< rb::CellHit > > const &trkHits, geo::View_t const &view, double &chiSqr, bool const &dir, std::vector< art::Ptr< rb::CellHit > > &removedHits)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
unsigned int FitWindow(std::vector< art::Ptr< rb::CellHit > > const &prevWindow, std::vector< art::Ptr< rb::CellHit > > &currWindow, std::vector< art::Ptr< rb::CellHit > > &removedHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, geo::View_t const &view)
void hits()
Definition: readHits.C:15
TDirectory * dir
Definition: macro.C:5
unsigned int fWindowSize
Number of planes in the sliding window.
unsigned int trk::WindowTrackingAlg::FitWindow ( std::vector< art::Ptr< rb::CellHit > > const &  prevWindow,
std::vector< art::Ptr< rb::CellHit > > &  currWindow,
std::vector< art::Ptr< rb::CellHit > > &  removedHits,
std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &  planeCoG,
geo::View_t const &  view 
)
private

Definition at line 649 of file WindowTrackingAlg.cxx.

References std::abs(), b, plot_validation_datamc::c, geo::PlaneGeo::Cell(), geo::DsqrToLine(), e, febshutoff_auto::end, fDHitGood, fGeo, geo::CellGeo::GetCenter(), geo::kX, geo::LinFitMinDperp(), getGoodRuns4SAM::n, genie::units::pb, geo::GeometryBase::Plane(), TB_WatchdogFx::pw, util::sqr(), febshutoff_auto::start, registry_explorer::v, w, submit_syst::x, and submit_syst::y.

Referenced by FitView().

655  {
656  size_t prevWinSize = prevWindow.size();
657 
658  double start[2] = {0.};
659  double end[2] = {0.};
660 
661  std::vector<double> x(prevWinSize + currWindow.size(), 0.);
662  std::vector<double> y(prevWinSize + currWindow.size(), 0.);
663  std::vector<double> w(prevWinSize + currWindow.size(), 0.);
664  std::vector<size_t> planeBoundaries;
665  double chiSqr = 0.;
666 
667  // fill the vectors, first with the previous window, then the current one
668  double xyz[3] = {0.};
669  double &v = (view == geo::kX) ? xyz[0] : xyz[1];
670  for(size_t c = 0; c < prevWinSize; ++c){
671  fGeo->Plane(prevWindow[c]->Plane())->Cell(prevWindow[c]->Cell())->GetCenter(xyz);
672  x[c] = xyz[2];
673  y[c] = v;
674  w[c] = 2.;//prevWindow[c]->ADC();
675  //LOG_VERBATIM("WindowTrackingAlg")
676  //<< "previous window points: "
677  //<< prevWindow[c]->Plane() << " "
678  //<< prevWindow[c]->Cell() << " "
679  //<< x[c] << " " << y[c] << " "
680  //<< w[c];
681 
682  }
683 
684  // fit the previous window to be able to figure out which cells
685  // in the current window should go in the initial fit
686  if(prevWinSize > 0)
687  chiSqr = geo::LinFitMinDperp(x, y, w,
688  &start[0], &start[1],
689  &end[0], &end[1]);
690 
691  unsigned int ndrop = 0;
692 
693  // variables to help identify outliers early on
694  unsigned short int currPlane = currWindow.front()->Plane();
695  unsigned short int prevPlane = currWindow.front()->Plane();
696  float cog = 0;
697  for(size_t c = 0; c < currWindow.size(); ++c){
698  currPlane = currWindow[c]->Plane();
699  fGeo->Plane(currPlane)->Cell(currWindow[c]->Cell())->GetCenter(xyz);
700  cog = planeCoG.find(currPlane)->second.first;
701  cog /= planeCoG.find(currPlane)->second.second;
702  x[c + prevWinSize] = xyz[2];
703  y[c + prevWinSize] = v;
704  w[c + prevWinSize] = 1.;//currWindow[c]->ADC();
705 
706  if(currPlane != prevPlane) planeBoundaries.push_back(c);
707 
708  prevPlane = currPlane;
709 
710  // First pass, drop hits that are clear outliers (more than 2x tolerance
711  // from anything else) in the current window. If it turns out they really
712  // are good hits (eg surrounded by bad channels) they will be added back
713  // in by the last step. But this helps prevent distant noise hits from
714  // throwing off the original fit.
715  if(prevWinSize > 0){
716  const double dSq = geo::DsqrToLine(x[c + prevWinSize],
717  y[c + prevWinSize],
718  start[0], start[1],
719  end[0], end[1]);
720 
721  if(dSq > util::sqr(fDHitGood*2)){
722  //LOG_VERBATIM("WindowTrackingAlg")
723  //<< "dropping hit at (z,v) = ("
724  //<< x[c + prevWinSize] << ","
725  //<< y[c + prevWinSize]
726  //<< ") for being too separated from others: "
727  //<< dSq
728  //<< " line start: " << start[0] << " " << start[1]
729  //<< " line end: " << end[0] << " " << end[1]
730  //<< " expected pos "
731  //<< x[c + prevWinSize]*(end[1]-start[1])/(end[0]-start[0]) + (start[1] - start[0]*(end[1]-start[1])/(end[0]-start[0]));
732 
733  w[c + prevWinSize] = 0;
734  ++ndrop;
735  }
736  }// end if we have a previous window to check this hit against
737  else{
738  // set the weight to 0 of cells from the same plane
739  if(std::abs(v - cog) > fDHitGood*4){
740  w[c + prevWinSize] = 0.;
741  ++ndrop;
742  }
743  }
744 
745  //LOG_VERBATIM("WindowTrackingAlg")
746  //<< "current window points: "
747  //<< currWindow[c]->Plane() << " "
748  //<< currWindow[c]->Cell() << " "
749  //<< x[c + prevWinSize] << " "
750  //<< y[c + prevWinSize] << " "
751  //<< w[c + prevWinSize] << " "
752  //<< cog << " "
753  //<< currWindow[c]->ADC();
754 
755  }
756 
757  // check if we have dropped all the hits from the current window
758  // and have none in the previous window. In that case, we are
759  // on the first window of the given slice but have no way to do
760  // a fit. Remove all the hits from the current window
761  if(ndrop == currWindow.size() && prevWinSize < 1){
762  currWindow.clear();
763  return ndrop;
764  }
765 
766  // sometimes you get a plane in the initial window that has
767  // its hits well away from the others and that throws off the
768  // fit. Run through the initial window and drop each plane's
769  // cells in turn. If a plane causes the chi^2 to spike drop it
770  if(prevWinSize < 1){
771 
772  // add the final boundary to the list
773  planeBoundaries.push_back(w.size());
774 
775  size_t prevPB = 0;
776  double bestChi2 = 1.e6;
777  double worstChi2 = 0.;
778  std::pair<size_t, size_t> bestBounds(0, 0);
779  for(auto const p : planeBoundaries){
780 
781  // set the weights for the plane to 0
782  std::vector< std::pair<size_t, double> > prevWeights;
783  for(auto pb = prevPB; pb < p; ++pb){
784  prevWeights.push_back(std::pair<size_t, double>(pb, w[pb]));
785  w[pb] = 0.;
786  }
787 
788  // for(size_t c = prevWinSize; c < x.size(); ++c)
789  // LOG_DEBUG("WindowTrackingAlg") << "look for outliers: " << p << " "
790  // << x[c + prevWinSize] << " "
791  // << y[c + prevWinSize] << " "
792  // << w[c + prevWinSize];
793 
794  try{
795  chiSqr = geo::LinFitMinDperp(x, y, w,
796  &start[0], &start[1],
797  &end[0], &end[1]);
798  }
799  catch(cet::exception &e){
800 
801  //LOG_DEBUG("WindowTrackingAlg") << "caught exception " << e
802  // << " and moving on to next plane";
803 
804  // there were likely no non-zero weights, so just move on to the next plane
805  continue;
806  }
807 
808  //LOG_DEBUG("WindowTrackingAlg") << "chisq = " << chiSqr;
809 
810  if(chiSqr < bestChi2){
811  bestChi2 = chiSqr;
812  bestBounds.first = prevPB;
813  bestBounds.second = p;
814  //LOG_DEBUG("WindowTrackingAlg") << "best chisqr, bounds: " << bestBounds.first << " " << bestBounds.second;
815  }
816  if(chiSqr > worstChi2) worstChi2 = chiSqr;
817 
818  // set the weights back to 1 to evaluate the remaining planes
819  for(auto pw : prevWeights) w[pw.first] = pw.second;
820 
821  // set the previous boundary to the current value
822  prevPB = p;
823 
824  }// end loop over plane boundaries
825 
826  // set the worst plane contributions to 0 weight for now
827  // but only do that if it adds 2 units of chi^2 or more
828  // for each cell in the window
829  if(worstChi2 - bestChi2 > 2.*w.size())
830  for(auto b = bestBounds.first; b < bestBounds.second; ++b) w[b] = 0.;
831 
832  }// end if on the first window for a slice in the chosen view
833 
834  bool droppedAny = true;
835  while(droppedAny){
836 
837  try{
838  chiSqr = geo::LinFitMinDperp(x, y, w,
839  &start[0], &start[1],
840  &end[0], &end[1]);
841  }
842  catch(cet::exception &e){
843 
844  // there were likely no non-zero weights, so stop trying to do a fit
845  break;
846  }
847 
848  //LOG_DEBUG("WindowTracking") << "end points from fit "
849  // << start[0] << " "
850  // << start[1] << " "
851  // << end[0] << " "
852  // << end[1] << " chi^2 = "
853  // << chiSqr;
854 
855  // loop over the hits and for the ones that are more than fDHitGood cm
856  // off the line set the weight to be 0 so that they do not contribute to
857  // future fits. Only drop the worst hit, and then iterate to get the next
858  // worst and so on. Otherwise we can end up dropping lots of hits we
859  // shouldn't due to one outlier.
860  double dSqMax = 0;
861  int dSqMaxIdx = -1;
862  for(size_t n = prevWinSize; n < x.size(); ++n){
863  // Already dropped
864  if(w[n] == 0) continue;
865 
866  const double dSq = geo::DsqrToLine(x[n], y[n],
867  start[0], start[1],
868  end[0], end[1]);
869  if(dSq > dSqMax){
870  dSqMax = dSq;
871  dSqMaxIdx = n;
872  //LOG_DEBUG("WindowTracking") << "highest dSq is at: "
873  // << x[n] << " " << y[n]
874  // << " " << dSq;
875  }
876  } // end for n
877 
878  droppedAny = false;
879  if(dSqMax > util::sqr(fDHitGood)){
880  //LOG_VERBATIM("WindowTracking")
881  //<< "dropping hit at ("
882  //<< x[dSqMaxIdx] << ", "
883  //<< y[dSqMaxIdx] << ")"
884  //<< " dSq = " << dSqMax << " expected "
885  //<< x[dSqMaxIdx]*(end[1]-start[1])/(end[0]-start[0]) + (start[1] - start[0]*(end[1]-start[1])/(end[0]-start[0]));
886 
887  w[dSqMaxIdx] = 0.;
888  droppedAny = true;
889  ++ndrop;
890  }
891  }// end while
892 
893  for(size_t n = prevWinSize; n < x.size(); ++n){
894  // Already part of the fit
895  if(w[n] != 0) continue;
896 
897  const double dSq = geo::DsqrToLine(x[n], y[n],
898  start[0], start[1],
899  end[0], end[1]);
900  if(dSq < util::sqr(fDHitGood)){
901  //LOG_DEBUG("WindowTracking")
902  //<< "adding back hit at ("
903  //<< x[n] << ", " << y[n]
904  //<< ")" << " dSq = " << dSq;
905 
906  w[n] = 1.;
907  --ndrop;
908  }
909  } // end for n
910 
911  // remove any hits from the current window that have a weight of 0
912  // erase all entries from w that are from the previous window
913  auto witr = w.begin();
914  if(prevWinSize > 0)
915  witr = w.erase(w.begin(), w.begin()+prevWinSize);
916  auto citr = currWindow.begin();
917  while( witr != w.end() && citr != currWindow.end()){
918  if( (*witr) < 1. ){
919  //LOG_VERBATIM("WindowTrackingAlg")
920  //<< "dropping hit at "
921  //<< (*citr)->Plane() << ","
922  //<< (*citr)->Cell();
923  // add the hit to the removed hits vector
924  removedHits.push_back(*citr);
925 
926  citr = currWindow.erase(citr);
927  witr = w.erase(witr);
928  }
929  else{
930  ++witr;
931  ++citr;
932  }
933  }
934 
935  //LOG_DEBUG("WindowTracking") << "end points from fit "
936  // << start[0] << " " << start[1] << " "
937  // << end[0] << " " << end[1];
938 
939  return ndrop;
940  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const PlaneGeo * Plane(unsigned int i) const
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
float abs(float number)
Definition: d0nt_math.hpp:39
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
static const double pb
Definition: Units.h:90
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
const hit & b
Definition: hits.cxx:21
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
void trk::WindowTrackingAlg::LookForBremsstrahlungHits ( rb::Track track,
std::vector< art::Ptr< rb::CellHit > > &  removedHits 
)
private

Definition at line 943 of file WindowTrackingAlg.cxx.

References std::abs(), rb::Cluster::Add(), plot_validation_datamc::c, fDHitGood, fGeo, geo::kX, geo::kY, geo::GeometryBase::Plane(), rb::Track::Trajectory(), registry_explorer::v, and rb::Cluster::View().

Referenced by MakeViewTrack().

945  {
946  // map the z position of each trajectory point to the position in the view
947  std::map<double, double> zToViewPos;
948 
949  for(auto tp : track.Trajectory()){
950  if (track.View() == geo::kX) zToViewPos[tp.Z()] = tp.X();
951  else if(track.View() == geo::kY) zToViewPos[tp.Z()] = tp.Y();
952  }// end loop to map z to trajectory in view
953 
954  // loop over the rejected cells to see if any of them are close enough
955  // to the other cells in the plane to be added back to the track
956  double xyz[3] = {0.};
957  double &v = (track.View() == geo::kX) ? xyz[0] : xyz[1];
958  auto rhitr = removedHits.begin();
959  bool useHit = false;
960  while( rhitr != removedHits.end() ){
961  useHit = false;
962  auto c = *rhitr;
963  fGeo->Plane(c->Plane())->Cell(c->Cell())->GetCenter(xyz);
964 
965  if(zToViewPos.count(xyz[2]) > 0){
966 
967  // if the cell is close enough to the trajectory at this location
968  // add it to the track and remove it from the
969  if(std::abs(zToViewPos[xyz[2]] - v) < 8.*fDHitGood){
970  track.Add(c);
971  rhitr = removedHits.erase(rhitr);
972  useHit = true;
973  }
974  }// end if the z position is in the map
975 
976  // the erase method returns an iterator to the position
977  // after the removed item, ie effectively incrementing the iterator
978  // if we didn't remove the hit, we need to increment the iterator
979  if(!useHit) ++rhitr;
980 
981  }// end loop over previously rejected hits
982 
983  return;
984  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const PlaneGeo * Plane(unsigned int i) const
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
rb::Track trk::WindowTrackingAlg::Make3DTrack ( std::vector< art::Ptr< rb::CellHit > > &  allHits,
std::vector< art::Ptr< rb::CellHit > > &  removedHits,
double &  chiSqr,
bool(*)(unsigned short int, unsigned short int comp 
)
private

Definition at line 298 of file WindowTrackingAlg.cxx.

References std::abs(), CheckTrackDirectionInY(), geo::ClampRayToDetectorVolume(), DetermineInitialDirection(), e, febshutoff_auto::end, fGeo, FitView(), fMaxPlaneSeparation, fMinViewPlanes, makeTrainCVSamples::int, geo::kX, geo::kY, LOG_DEBUG, LOG_WARNING, next(), geo::GeometryBase::Plane(), prev(), ShortTrack(), febshutoff_auto::start, rb::Prong::Start(), rb::Track::Stop(), registry_explorer::v, ValidTrack(), rb::Cluster::View(), and rb::Track::ZipWith().

Referenced by MakeTrack().

302  {
303  // loop over all the hits for this slice and get the
304  // x or y cells
305  std::map<unsigned short int, std::vector< art::Ptr<rb::CellHit> >,
306  bool(*)(unsigned short int, unsigned short int) > xViewHits(comp);
307  std::map<unsigned short int, std::vector< art::Ptr<rb::CellHit> >,
308  bool(*)(unsigned short int, unsigned short int) > yViewHits(comp);
309  std::map<unsigned short int, std::pair<float, float>,
310  bool(*)(unsigned short int, unsigned short int) > planeCoG(comp);
311 
312  double xyz[3] = {0.};
313  double v = 0.;
314 
315  // the allCells vector is already sorted properly.
316  std::set<int> planes;
317 
318  //LOG_VERBATIM("WindowTrackingAlg") << "starting with these hits:";
319  for(auto const chit : allCells){
320 
321  planes.insert(chit->Plane());
322 
323  fGeo->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
324 
325  //LOG_VERBATIM("WindowTrackingAlg") << "\t" << chit->Plane() << " " << chit->Cell();
326 
327  if(chit->View() == geo::kX){
328  xViewHits[chit->Plane()].push_back(chit);
329  v = xyz[0];
330  }
331  if(chit->View() == geo::kY){
332  yViewHits[chit->Plane()].push_back(chit);
333  v = xyz[1];
334  }
335 
336  if(planeCoG.find(chit->Plane()) != planeCoG.end()){
337  planeCoG[chit->Plane()].first += chit->ADC()*v;
338  planeCoG[chit->Plane()].second += chit->ADC();
339  }
340  else
341  planeCoG[chit->Plane()] = std::pair<float,float>(chit->ADC()*v,
342  chit->ADC());
343 
344  }// end loop over hits in this slice
345 
346  // loop over planeCoG and see what we have
347  //for(auto pcog : planeCoG)
348  // LOG_VERBATIM("WindowTrackingAlg")
349  // << pcog.first
350  // << " "
351  // << pcog.second.first/pcog.second.second;
352 
353  // remove planes that are too far separated from others in the slice
354  auto prev = planes.begin();
355  auto next = prev;
356  bool removeThisPlane = false;
357 
358  for(auto itr = planes.begin(); itr != planes.end(); ++itr){
359  next = itr; ++next;
360  removeThisPlane = false;
361 
362  // is the first plane separated from the second by too much?
363  if(itr == planes.begin() && next != planes.end()){
364  if(std::abs(*itr - *next) > fMaxPlaneSeparation) removeThisPlane = true;
365  }// end if on the first plane in the set
366  // is the current plane separated by both the previous and next one
367  // too much
368  else if(itr != planes.begin() && next != planes.end()){
369  prev = itr; --prev;
370  if(std::abs(*itr - *prev) > fMaxPlaneSeparation &&
371  std::abs(*itr - *next) > fMaxPlaneSeparation) removeThisPlane = true;
372  }
373  else if(next == planes.end()){
374  prev = itr; --prev;
375  if(std::abs(*itr - *prev) > fMaxPlaneSeparation) removeThisPlane = true;
376  }
377 
378  if(removeThisPlane){
379  xViewHits.erase(*itr);
380  yViewHits.erase(*itr);
381  planeCoG .erase(*itr);
382  // LOG_VERBATIM("WindowTrackingAlg") << "removing plane " << *itr << " from the fit; prev: "
383  // << *prev << " next: " << *next;
384  }
385 
386  }//end loop to look for big gaps between planes.
387 
388  rb::Track retTrack;
389 
390  if(this->ValidTrack(xViewHits) &&
391  this->ValidTrack(yViewHits) ){
392 
393  double chiSqrX = 0.;
394  double chiSqrY = 0.;
395 
396  rb::Track xTrack = FitView(xViewHits, planeCoG, removedHits, geo::kX, chiSqrX);
397  rb::Track yTrack = FitView(yViewHits, planeCoG, removedHits, geo::kY, chiSqrY);
398 
399  chiSqr = chiSqrX + chiSqrY;
400 
401  if(xTrack.View() != geo::kX || yTrack.View() != geo::kY){
402  LOG_DEBUG("WindowTrackingAlg") << "no track found in at least one view, return empty track";
403  return retTrack;
404  }
405 
406  //LOG_DEBUG("WindowTrackingAlg") << " x track view: " << xTrack.View()
407  // << " y track view: " << yTrack.View();
408 
409  // try making a 3D track from these two 2D tracks
410  try{
411  retTrack = xTrack.ZipWith(yTrack);
412  }
413  catch(cet::exception &e){
414  LOG_WARNING("WindowTrackingAlg") << "attempt to zip tracks failed with exception\n"
415  << e
416  << "\n return empty track";
417  return retTrack;
418  }
419 
420  TVector3 start(retTrack.Start());
421  TVector3 end (retTrack.Stop());
422 
423  this->CheckTrackDirectionInY(retTrack, start, end);
424 
425  // determine the direction of the track by fitting the first
426  // several trajectory points - just using the direction from the
427  // first point to the next can be misleading - need a bigger lever arm
428  this->DetermineInitialDirection(retTrack);
429 
430  // restrict track to exist inside detector; if this isn't possible,
431  // the track is no good, so don't write
433  //std::cout << "-------------------------------" << std::endl;
434  //for(auto const& traj : retTrack.Trajectory()) traj.Print();
435  return retTrack;
436  }
437 
438  } // end if we can make a 3D track
439  // see if we can make a short track - only worry about tracks that have at most fMinViewPlanes
440  // in each view - if one view has 12 planes and the other 2, it is probably garbage anyway
441  else if(xViewHits.size() < fMinViewPlanes &&
442  yViewHits.size() < fMinViewPlanes &&
443  xViewHits.size() > 0 &&
444  yViewHits.size() > 0 ){
445  //LOG_VERBATIM("WindowTrackingAlg")
446  //<< "attempting to fit a short track: "
447  //<< xViewHits.size() << " " << yViewHits.size();
448 
449  return this->ShortTrack(xViewHits, yViewHits, planeCoG, chiSqr);
450  }
451  else
452  // can't make a track, so set the chiSqr value to something stupidly big
453  chiSqr = 1.e9;
454 
455  return retTrack;
456  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Vertical planes which measure X.
Definition: PlaneGeo.h:28
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
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
virtual TVector3 Start() const
Definition: Prong.h:73
rb::Track ShortTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &xViewHits, std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &yViewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, double &chiSqr)
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
rb::Track FitView(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &hits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, std::vector< art::Ptr< rb::CellHit > > &removedHits, geo::View_t view, double &chiSqr)
void CheckTrackDirectionInY(rb::Track &track, TVector3 &start, TVector3 &end)
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
void prev()
Definition: show_event.C:91
rb::Track ZipWith(const rb::Track &trk) const
Combine with a 2D track from the other view to make a 3D track.
Definition: Track.cxx:442
unsigned int fMinViewPlanes
Minimum number of planes in a view to fit.
#define LOG_WARNING(category)
void DetermineInitialDirection(rb::Track &track)
bool ValidTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &planeToHit)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
int fMaxPlaneSeparation
Maximum number of planes allowed to separate planes with hits in slice.
void next()
Definition: show_event.C:84
std::vector< rb::Track > trk::WindowTrackingAlg::MakeTrack ( const rb::Cluster slice)
virtual

Implements trk::TrackAlg.

Definition at line 74 of file WindowTrackingAlg.cxx.

References std::abs(), rb::Cluster::Add(), rb::Cluster::AllCells(), plot_validation_datamc::c, getBrightness::cell, geo::PlaneGeo::Cell(), fCellHalfDepth, fDHitGood, fGeo, fMinHitsTryAgain, gt(), geo::CellGeo::HalfD(), lt(), Make3DTrack(), rb::Cluster::NCell(), geo::GeometryBase::Plane(), r(), APDHVSetting::removed, SetTrackEndPoints(), art::PtrVector< T >::sort(), track, and Z.

75  {
76  if(fCellHalfDepth < 0.) fCellHalfDepth = fGeo->Plane(0)->Cell(0)->HalfD();
77 
78 
79  //LOG_DEBUG("WindowTrackingAlg") << "fitting new slice";
80 
81  std::vector<rb::Track> tracks;
82 
83  // get the PtrVector of all cells in this slice
84  art::PtrVector<rb::CellHit> allCells = slice->AllCells();
85  allCells.sort();
86  std::vector<art::Ptr<rb::CellHit> > allHitsFor;
87  std::vector<art::Ptr<rb::CellHit> > allHitsRev;
88 
89  for(auto cell : allCells){
90  allHitsFor.push_back(cell);
91  allHitsRev.push_back(cell);
92  }
93 
94  std::vector<art::Ptr<rb::CellHit> > removedHits;
95  std::vector<art::Ptr<rb::CellHit> > removedHitsFor;
96  std::vector<art::Ptr<rb::CellHit> > removedHitsRev;
97  bool tryAnother = true;
98 
99  double chiSqrFor = 0.;
100  double chiSqrRev = 0.;
101 
103 
104  while(tryAnother){
105  tryAnother = false;
106 
107  rb::Track trackForward = this->Make3DTrack(allHitsFor, removedHitsFor, chiSqrFor, lt);
108  rb::Track trackReverse = this->Make3DTrack(allHitsRev, removedHitsRev, chiSqrRev, gt);
109 
110  //LOG_DEBUG("WindowTrackingAlg") << "forward chi^2: " << chiSqrFor << "/" << trackForward.NCell()
111  // << "\nreverse chi^2: " << chiSqrRev << "/" << trackReverse.NCell();
112 
113  if(trackForward.NCell() > 0) chiSqrFor /= trackForward.NCell();
114  if(trackReverse.NCell() > 0) chiSqrRev /= trackReverse.NCell();
115 
116  if(chiSqrFor < chiSqrRev){
117  track = trackForward;
118  removedHits.swap(removedHitsFor);
119  }
120  else{
121  track = trackReverse;
122  removedHits.swap(removedHitsRev);
123  }
124 
125  removedHitsRev.clear();
126  removedHitsFor.clear();
127 
128  // check to see if we made a legit track, it will have at least 1 cell
129  if(track.NCell() > 1){
130  this->SetTrackEndPoints(track);
131 
132  tracks.push_back(track);
133 
134  /* Some debugging code commented out below, leaving it in place for future convenience
135  if(track.Dir().Z() > 0){
136  // get the cells for this track and sort them
137  art::PtrVector<rb::CellHit> trCells = track.AllCells();
138  trCells.sort();
139 
140  double xyzStart[3] = {0.};
141  double xyzEnd[3] = {0.};
142  fGeo->Plane(trCells[0]->Plane())->Cell(trCells[0]->Cell())->GetCenter(xyzStart);
143  fGeo->Plane(trCells[track.NCell()-1]->Plane())->Cell(trCells[track.NCell()-1]->Cell())->GetCenter(xyzEnd);
144 
145  LOG_VERBATIM("WindowTrackingAlg")
146  << "first cell " << xyzStart[2] //<< "," << *(trCells[0])
147  << "\nlast cell " << xyzEnd[2] //<< "," << *(trCells[track.NCell()-1])
148  << "\nstart: " << track.Start().X() << " " << track.Start().Y() << " " << track.Start().Z()
149  << "\nend: " << track.Stop().X() << " " << track.Stop().Y() << " " << track.Stop().Z();
150 
151  for(size_t h = 0; h < track.NCell(); ++h)
152  LOG_VERBATIM("WindowTrackingAlg") << *(track.Cell(h));
153  }*/
154 
155  if(removedHits.size() > fMinHitsTryAgain){
156  tryAnother = true;
157  //LOG_DEBUG("WindowTrackingAlg")
158  //<< "trying to make another track from the "
159  //<< removedHits.size()
160  //<< " removed hits";
161  allHitsFor.clear();
162  allHitsRev.clear();
163  for(auto removed : removedHits){
164  //LOG_DEBUG("WindowTrackingAlg")
165  //<< "removed "
166  //<< *removed;
167  allHitsFor.push_back(removed);
168  allHitsRev.push_back(removed);
169  }
170  removedHits.clear();
171  }// end if there are enough left over hits to look for another track
172  }// end if the returned track was valid
173 
174  }// end loop to look for tracks in the slice
175 
176  // check to see if any of the tracks should be merged
177  if(tracks.size() > 1){
178  //LOG_VERBATIM("WindowTrackingAlg")
179  //<< "there are "
180  //<< tracks.size()
181  //<< " tracks in the slice";
182  TVector3 t1Start, t2Start, t1End, t2End, diff12, diff21;
183  float t1ZDir = 0., t2ZDir = 0.;
184  size_t trks = tracks.size();
185  std::vector<bool> removed(trks, false);
186  for(size_t tr1 = 0; tr1 < trks; ++tr1){
187  // int ctr = 0;
188  if(removed[tr1]) continue;
189  for(size_t tr2 = 0; tr2 < trks; ++tr2){
190  if(tr1 == tr2 || removed[tr2]) continue;
191 
192  // LOG_VERBATIM("WindowTrackingAlg") << "inner loop " << ctr;
193  // ++ctr;
194 
195  // get the start and end points of each track
196  t1Start = tracks[tr1].Start();
197  t2Start = tracks[tr2].Start();
198  t1End = tracks[tr1].Stop();
199  t2End = tracks[tr2].Stop();
200 
201  t1ZDir = tracks[tr1].Dir().Z()/std::abs(tracks[tr1].Dir().Z());
202  t2ZDir = tracks[tr2].Dir().Z()/std::abs(tracks[tr2].Dir().Z());
203 
204  // LOG_VERBATIM("WindowTrackingAlg") << "t1 start: "
205  // << t1Start.X() << " " << t1Start.Y() << " " << t1Start.Z() << " end: "
206  // << t1End.X() << " " << t1End.Y() << " " << t1End.Z() << "\n"
207  // << "t2 start: "
208  // << t2Start.X() << " " << t2Start.Y() << " " << t2Start.Z() << " end: "
209  // << t2End.X() << " " << t2End.Y() << " " << t2End.Z();
210 
211  // now see if any of the end points match up
212  // that is they are within a cell of each other
213  // The total difference in the end points has to be
214  // less than 2*fDHitGood and the z component of that
215  // should be less than fDHitGood
216  // also make sure they are going in the same z direction
217  if(t1ZDir*t2ZDir > 0.){
218  diff12 = t1Start - t2End;
219  diff21 = t2Start - t1End;
220 
221  if( (diff12.Mag() < 2.*fDHitGood && std::abs(diff12.Z()) < fDHitGood) ||
222  (diff21.Mag() < 2.*fDHitGood && std::abs(diff21.Z()) < fDHitGood) ){
223 
224  // get the iterators in the correct order for the z direction of
225  // the track
226  auto firstTr = tracks[tr1];
227  auto secondTr = tracks[tr2];
228  if(t1Start.Z()*t1ZDir > t2Start.Z()*t2ZDir){
229  firstTr = tracks[tr2];
230  secondTr = tracks[tr1];
231  }
232 
233  // merge the tracks
234  std::vector<TVector3> mergedTrajectory;
235  // get a merged trajectory point for the intersection of the two tracks
236  // and set that to the be the first trajectory point in track 1, don't
237  // prepend last point from the trajectory of track 2
238  // the first two points and last two points in a track represent the edge
239  // of and the midpoint of the final plane, so remove the last two points
240  // in firstTr and the first two points in secondTr from the merged track
241  for(size_t tp = 0; tp < firstTr.Trajectory().size()-2; ++tp){
242  mergedTrajectory.push_back(firstTr.TrajectoryPoint(tp));
243  //LOG_VERBATIM("WindowTrackingAlg")
244  //<< mergedTrajectory.back().X() << " "
245  //<< mergedTrajectory.back().Y() << " "
246  //<< mergedTrajectory.back().Z();
247  }
248  mergedTrajectory.push_back(0.5*(secondTr.Trajectory().front() + firstTr.Trajectory().back()));
249  //LOG_VERBATIM("WindowTrackingAlg") << "connection point "
250  //<< mergedTrajectory.back().X() << " "
251  //<< mergedTrajectory.back().Y() << " "
252  //<< mergedTrajectory.back().Z();
253  for(size_t tp = 2; tp < secondTr.Trajectory().size(); ++tp){
254  mergedTrajectory.push_back(secondTr.TrajectoryPoint(tp));
255  //LOG_VERBATIM("WindowTrackingAlg")
256  //<< mergedTrajectory.back().X() << " "
257  //<< mergedTrajectory.back().Y() << " "
258  //<< mergedTrajectory.back().Z();
259  }
260 
261  // now create the merged track with cells from both tracks
262  rb::Track mergedTrack(firstTr.AllCells(), mergedTrajectory.front(), firstTr.Dir());
263  for(auto c : secondTr.AllCells()) mergedTrack.Add(c);
264  for(size_t tp = 1; tp < mergedTrajectory.size(); ++tp) mergedTrack.AppendTrajectoryPoint(mergedTrajectory[tp]);
265  tracks.push_back(mergedTrack);
266 
267  //LOG_VERBATIM("WindowTrackingAlg") << firstTr.Trajectory().size() << " "
268  //<< secondTr.Trajectory().size() << " "
269  //<< mergedTrajectory.size() << " trajectory points now: \n";
270  //for(auto tp : tracks.back().Trajectory())
271  // LOG_VERBATIM("WindowTrackingAlg") << tp.X() << " " << tp.Y() << " " << tp.Z();
272 
273  // remove both track 1 and 2
274  removed[tr1] = true;
275  removed[tr2] = true;
276  //LOG_VERBATIM("WindowTrackingAlg") << "removing tracks " << tr1 << " " << tr2;
277  }// end if the tracks end points are close enough to merge
278  }// end if tracks are going in the same direction in Z
279  }// end inner loop over tracks
280  }// end outer loop over tracks
281 
282  // loop over the removed vector and remove those tracks from the vector
283  size_t rmctr = 0;
284  for(size_t r = 0; r < removed.size(); ++r){
285  if(removed[r]){
286  tracks.erase(tracks.begin() + (r-rmctr));
287  ++rmctr;
288  }
289  }
290  //LOG_VERBATIM("WindowTrackingAlg") << "there are " << tracks.size() << " tracks in the slice";
291  }// end if more than one track in the slice
292 
293  return tracks;
294 
295  } // End of MakeTrack
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
unsigned int fMinHitsTryAgain
Minimum number of hits left to look for another track.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
float abs(float number)
Definition: d0nt_math.hpp:39
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Float_t Z
Definition: plot.C:38
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
rb::Track Make3DTrack(std::vector< art::Ptr< rb::CellHit > > &allHits, std::vector< art::Ptr< rb::CellHit > > &removedHits, double &chiSqr, bool(*comp)(unsigned short int, unsigned short int))
double fCellHalfDepth
half of cell size in z
void SetTrackEndPoints(rb::Track &track)
TRandom3 r(0)
bool gt(unsigned short int a, unsigned short int b)
Float_t track
Definition: plot.C:35
bool lt(unsigned short int a, unsigned short int b)
rb::Track trk::WindowTrackingAlg::MakeViewTrack ( std::vector< art::Ptr< rb::CellHit > > const &  trkHits,
geo::View_t const &  view,
double &  chiSqr,
bool const &  dir,
std::vector< art::Ptr< rb::CellHit > > &  removedHits 
)
private

Definition at line 556 of file WindowTrackingAlg.cxx.

References rb::Track::AppendTrajectoryPoint(), geo::PlaneGeo::Cell(), dz, fGeo, gt(), makeTrainCVSamples::int, geo::kX, LookForBremsstrahlungHits(), lt(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, cet::pow(), moon_position_table_new3::second, confusionMatrixTree::t, ValidTrack(), w, and weights.

Referenced by FitView().

561  {
562  //LOG_DEBUG("WindowTrackingAlg")
563  //<< "view: " << view
564  //<< " has " << trkHits.size()
565  //<< " cell hits";
566  rb::Track retTrack;
567 
568  // make some trajectory points for this track - use the signal weighted
569  // average point in each plane
570 
571  // first make a map of the planes to hits and then find the weighted average for
572  // each plane, set the direction in plane number based on the input bool
573  bool(*comp)(unsigned short int, unsigned short int) = (dir) ? gt : lt;
574  std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
575  bool(*)(unsigned short int, unsigned short int) > planesToHitLoc(comp);
576 
577  for(auto const hititr : trkHits)
578  planesToHitLoc[hititr->Plane()].push_back(hititr);
579 
580  // check that you ended up with a reasonable track in terms
581  // of the number of hit planes compared to the number of planes crossed
582  // set the chiSqr to something stupidly huge if that is not the case and
583  // return an empty track
584  if( !this->ValidTrack(planesToHitLoc) ){
585  chiSqr = 1.e6;
586  return retTrack;
587  }
588 
589  std::vector<std::pair<double, double> > traj;
590  double xyz[3] = {0.};
591  double &w = (view == geo::kX) ? xyz[0] : xyz[1];
592  double avg = 0.;
593  chiSqr = 0.;
594  unsigned short int plane = 0;
595 
596  for(auto const mapitr : planesToHitLoc){
597  double weightedCells = 0.;
598  double weights = 0.;
599  plane = mapitr.first;
600  for(auto const cellitr : mapitr.second){
601  //LOG_VERBATIM("WindowTrackingAlg")
602  //<< "makeviewtrack: "
603  //<< cellitr->Plane()
604  //<< " " << cellitr->Cell()
605  //<< " " << cellitr->ADC();
606  fGeo->Plane(plane)->Cell(cellitr->Cell())->GetCenter(xyz);
607  weightedCells += w*cellitr->ADC();
608  weights += 1.*cellitr->ADC();
609 
610  } //end loop over cell hits in the current plane
611 
612  // the z position is the same for all cells in a plane, and we are looping
613  // over all cells in each plane
614  if(weights > 0.){
615  avg = weightedCells/weights;
616 
617  for(auto const cellitr : mapitr.second){
618  fGeo->Plane(plane)->Cell(cellitr->Cell())->GetCenter(xyz);
619  chiSqr += std::pow(avg - w, 2.)/3.8; // cell height is 3.8 cm
620  }
621 
622  // LOG_DEBUG("WindowTrackingAlg") << "adding trajectory point in "
623  // << view << " view: " << xyz[2] << " " << avg;
624 
625  traj.push_back(std::pair<double, double>(avg, xyz[2]));
626  }
627 
628  } // end loop over planesToHitLoc
629 
630 
631  double dz = traj[1].second - traj[0].second;
632  double dv = traj[1].first - traj[0].first;
633 
634  //LOG_DEBUG("WindowTrackingAlg") << "there are " << traj.size() << " trajectory points in this view";
635 
636  retTrack = rb::Track(trkHits, view, traj[0].first, traj[0].second, dv, dz);
637  for(size_t t = 1; t < traj.size(); ++t)
638  retTrack.AppendTrajectoryPoint(traj[t].first, traj[t].second);
639 
640  // look for bremstrahlung hits to add back in
641  this->LookForBremsstrahlungHits(retTrack, removedHits);
642 
643  return retTrack;
644  }
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
double dz[NP][NC]
void LookForBremsstrahlungHits(rb::Track &track, std::vector< art::Ptr< rb::CellHit > > &removedHits)
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
Var weights
TDirectory * dir
Definition: macro.C:5
bool ValidTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &planeToHit)
bool gt(unsigned short int a, unsigned short int b)
bool lt(unsigned short int a, unsigned short int b)
Float_t w
Definition: plot.C:20
void trk::WindowTrackingAlg::SetTrackEndPoints ( rb::Track track)
private

Definition at line 1382 of file WindowTrackingAlg.cxx.

References rb::Cluster::AllCells(), rb::Prong::Dir(), fGeo, FindEndPoint(), geo::kX, makeBrightnessMap::maxPlane, rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), rb::Track::NTrajectoryPoints(), geo::GeometryBase::Plane(), sortCellPtr(), rb::Prong::Start(), rb::Track::Stop(), rb::Track::StopDir(), std::swap(), track, rb::Track::Trajectory(), and geo::PlaneGeo::View().

Referenced by MakeTrack().

1383  {
1384  // look for any trajectory points that have a Nan
1385 // std::cout << " checking out the track before endpoints set" << std::endl;
1386 // for(auto tp : track.Trajectory())
1387 // if(std::isnan(tp.X()) ||
1388 // std::isnan(tp.Y()) ||
1389 // std::isnan(tp.Z()) ) tp.Print();
1390 
1391  // extrapolate the track along the direction cosines
1392  // first get the cells from the first and last planes
1393  std::vector<art::Ptr<rb::CellHit> > startPlaneCells;
1394  std::vector<art::Ptr<rb::CellHit> > endPlaneCells;
1395  unsigned int minPlane = track.MinPlane();
1396  unsigned int maxPlane = track.MaxPlane();
1397  for(auto ch : track.AllCells() ){
1398  if(ch->Plane() == maxPlane) endPlaneCells.push_back(ch);
1399  if(ch->Plane() == minPlane) startPlaneCells.push_back(ch);
1400  }
1401 
1402  std::sort(startPlaneCells.begin(), startPlaneCells.end(), sortCellPtr);
1403  std::sort(endPlaneCells.begin(), endPlaneCells.end(), sortCellPtr);
1404 
1405  // figure out the direction from the next to last point to the last point
1406  TVector3 endPt = track.Trajectory().back();
1407  TVector3 endDir = endPt - track.Trajectory()[track.NTrajectoryPoints()-2];
1408 
1409  // figure out the direction from the second point to the first point
1410  TVector3 startPt = track.Trajectory().front();
1411  TVector3 startDir = track.Trajectory()[1] - startPt;
1412 
1413  geo::View_t endPlaneView = fGeo->Plane(maxPlane)->View();
1414  geo::View_t startPlaneView = fGeo->Plane(minPlane)->View();
1415 
1416  // determining which is the startPlane of the track and the endPlane
1417  // depends on the z direction. Swap the min and max vectors if the
1418  // direction cosine in z is negative
1419  if(endDir.Z() < 0.){
1420  startPlaneCells.swap(endPlaneCells);
1421  std::swap(startPlaneView, endPlaneView);
1422  }
1423 
1424  // get the highest cell in the start plane and the lowest in the
1425  // end plane
1426  unsigned short int startCell = startPlaneCells.back()->Cell();
1427  unsigned short int endCell = endPlaneCells.front()->Cell();
1428  // if the start or end plane is in the X view, then we need to
1429  // determine the extrema for the cells based on the track
1430  // trajectory
1431  if(startPlaneView == geo::kX && startDir.X() > 0.) startCell = startPlaneCells.front()->Cell();
1432  if(endPlaneView == geo::kX && endDir.X() > 0.) endCell = endPlaneCells.back()->Cell();
1433 
1434  // point the track back from the starting trajectory point to
1435  // a better location accounting for the potential difference between
1436  // the center of gravity in a plane and its final cell in the plane
1437  double xyz[3] = {0.};
1438  fGeo->Plane(startPlaneCells.front()->Plane())->Cell(startCell)->GetCenter(xyz);
1439  this->FindEndPoint(startPlaneView, xyz, track.Start(), track.Dir(), true, track);
1440 
1441  // carry the track further in its direction of travel
1442  fGeo->Plane(endPlaneCells.front()->Plane())->Cell(endCell)->GetCenter(xyz);
1443  this->FindEndPoint(endPlaneView, xyz, track.Stop(), track.StopDir(), false, track);
1444 
1445  return;
1446  }
size_t NTrajectoryPoints() const
Definition: Track.h:83
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
bool sortCellPtr(art::Ptr< rb::CellHit > const a, art::Ptr< rb::CellHit > const b)
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void FindEndPoint(geo::View_t const &view, double *cellCenter, TVector3 const &endPoint, TVector3 const &endDir, bool const &start, rb::Track &track)
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
virtual TVector3 StopDir() const
Get direction at the stop point.
Definition: Track.cxx:193
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t track
Definition: plot.C:35
rb::Track trk::WindowTrackingAlg::ShortTrack ( std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &  xViewHits,
std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &  yViewHits,
std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &  planeCoG,
double &  chiSqr 
)
private

Definition at line 1030 of file WindowTrackingAlg.cxx.

References std::abs(), rb::Cluster::Add(), rb::Cluster::AllCells(), rb::Track::AppendTrajectoryPoint(), CheckTrackDirectionInY(), geo::ClampRayToDetectorVolume(), rb::Prong::Dir(), fGeo, art::PtrVector< T >::front(), gtvec(), geo::kX, geo::kY, ltvec(), geo::GeometryBase::Plane(), ShortTrackExtraPlane(), ShortTrackRemoveExtraPlane(), ShortViewTrack(), and POTSpillRate::view.

Referenced by Make3DTrack().

1037  {
1038 
1039  // the assumption for these tracks is that they are very steep and leave lots of
1040  // hits in only at most 3 planes in a view. Just use the two planes with the
1041  // most hits to set the track direction and then add in the hits from the third
1042  // plane if they fall along that line
1043  // make a copy of the input maps and then drop the planes with
1044  // the fewest hits in each view. At this point there should
1045  // be only two planes in each view
1046  auto cogCopy(planeCoG);
1047  auto xHitsCopy(xViewHits);
1048  auto yHitsCopy(yViewHits);
1049  unsigned short int removedXPlane = 0;
1050  unsigned short int removedYPlane = 0;
1051  if(xViewHits.size() > 2) this->ShortTrackRemoveExtraPlane(xHitsCopy, cogCopy, removedXPlane);
1052  if(yViewHits.size() > 2) this->ShortTrackRemoveExtraPlane(yHitsCopy, cogCopy, removedYPlane);
1053 
1054  // make sure the two plane hit view has the planes separated by only 1 plane in the other view, and that
1055  // the one hit plane view is between them
1056  int firstXPlane = xHitsCopy.cbegin()->first;
1057  int secondXPlane = 0;
1058  int firstYPlane = yHitsCopy.cbegin()->first;
1059  int secondYPlane = 0;
1060  float firstXCoG = cogCopy.find(firstXPlane)->second.first/cogCopy.find(firstXPlane)->second.second;
1061  float firstYCoG = cogCopy.find(firstYPlane)->second.first/cogCopy.find(firstYPlane)->second.second;
1062  float secondXCoG = 0.;
1063  float secondYCoG = 0.;
1064  double xyz[3] = {0.};
1065  fGeo->Plane(xHitsCopy.cbegin()->first)->Cell((xHitsCopy.cbegin()->second.front())->Cell())->GetCenter(xyz);
1066  float firstXZ = xyz[2];
1067  fGeo->Plane(yHitsCopy.cbegin()->first)->Cell((yHitsCopy.cbegin()->second.front())->Cell())->GetCenter(xyz);
1068  float firstYZ = xyz[2];
1069  float secondXZ = xyz[2];
1070  float secondYZ = xyz[2];
1071  bool twoXPlanes = (xHitsCopy.size() == 2) ? true : false;
1072  bool twoYPlanes = (yHitsCopy.size() == 2) ? true : false;
1073 
1074  rb::Track retTrack;
1075 
1076 
1077  if(twoXPlanes){
1078  secondXPlane = xHitsCopy.crbegin()->first;
1079  secondXCoG = cogCopy.find(secondXPlane)->second.first/cogCopy.find(secondXPlane)->second.second;
1080  fGeo->Plane(xHitsCopy.crbegin()->first)->Cell((xHitsCopy.crbegin()->second.front())->Cell())->GetCenter(xyz);
1081  secondXZ = xyz[2];
1082  }
1083  if(twoYPlanes){
1084  secondYPlane = yHitsCopy.crbegin()->first;
1085  secondYCoG = cogCopy.find(secondYPlane)->second.first/cogCopy.find(secondYPlane)->second.second;
1086  fGeo->Plane(yHitsCopy.crbegin()->first)->Cell((yHitsCopy.crbegin()->second.front())->Cell())->GetCenter(xyz);
1087  secondYZ = xyz[2];
1088  }
1089 // LOG_VERBATIM("WindowTrackingAlg") << firstXZ << " " << firstXCoG << " " << secondXZ << " " << secondXCoG;
1090 // LOG_VERBATIM("WindowTrackingAlg") << firstYZ << " " << firstYCoG << " " << secondYZ << " " << secondYCoG;
1091 
1092 
1093  bool makeTrack = false;
1094 
1095  // check that the planes in each view are contiguous
1096  if(twoXPlanes && !twoYPlanes){
1097  if(std::abs(firstXPlane - secondXPlane) < 3 &&
1098  std::abs(firstXPlane - firstYPlane) < 2 &&
1099  std::abs(secondXPlane - firstYPlane) < 2)
1100  makeTrack = true;
1101  }
1102  else if(!twoXPlanes && twoYPlanes){
1103  if(std::abs(firstYPlane - secondYPlane) < 3 &&
1104  std::abs(firstYPlane - firstXPlane) < 2 &&
1105  std::abs(secondYPlane - firstXPlane) < 2)
1106  makeTrack = true;
1107  }
1108  else if(twoXPlanes && twoYPlanes){
1109  if(std::abs(firstYPlane - secondYPlane) < 3 &&
1110  std::abs(firstXPlane - secondXPlane) < 3 &&
1111  std::abs(firstYPlane - firstXPlane ) < 2 &&
1112  std::abs(secondYPlane - secondXPlane) < 2)
1113  makeTrack = true;
1114  }
1115 
1116  if(makeTrack){
1117 
1118 // LOG_VERBATIM("WindowTrackingAlg") << "making the short tracks in each view";
1119 
1120  rb::Track xTrack = this->ShortViewTrack(xHitsCopy, cogCopy, twoXPlanes, firstXPlane, secondXPlane, geo::kX);
1121  rb::Track yTrack = this->ShortViewTrack(yHitsCopy, cogCopy, twoYPlanes, firstYPlane, secondYPlane, geo::kY);
1122 
1123  // for a short track all the points are used, make the chiSqr 1
1124  chiSqr = 1.;
1125 
1126  // don't use ZipWith to make the 3D track - it is simple enough that
1127  // we can do it here
1128  // get the direction cosines in dx/dz and dy/dz
1129  TVector3 dxdz = xTrack.Dir();
1130  TVector3 dydz = yTrack.Dir();
1131 
1132  // figure out the trajectory points for each plane
1133  std::vector<std::pair<geo::View_t, TVector3> > trajPoints(4);
1134  trajPoints[0] = std::pair<geo::View_t, TVector3>(geo::kX, TVector3(firstXCoG, 0., firstXZ ));
1135  trajPoints[1] = std::pair<geo::View_t, TVector3>(geo::kY, TVector3(0., firstYCoG, firstYZ ));
1136  trajPoints[2] = std::pair<geo::View_t, TVector3>(geo::kX, TVector3(secondXCoG, 0., secondXZ));
1137  trajPoints[3] = std::pair<geo::View_t, TVector3>(geo::kY, TVector3(0., secondYCoG, secondYZ));
1138 
1139  // remove the second point in a view if it doesn't exist
1140  if(!twoXPlanes) trajPoints.erase(trajPoints.begin()+2);
1141  if(!twoYPlanes) trajPoints.erase(trajPoints.begin()+trajPoints.size()-1);
1142 
1143  // check if we need to add in an extra plane in either view
1144  if(xViewHits.size() > 2) this->ShortTrackExtraPlane(xViewHits, planeCoG, trajPoints, xTrack, removedXPlane);
1145  if(yViewHits.size() > 2) this->ShortTrackExtraPlane(yViewHits, planeCoG, trajPoints, yTrack, removedYPlane);
1146 
1147  // now sort the list according to the z value of the planes
1148  if(cogCopy.cbegin()->first < cogCopy.crbegin()->first)
1149  std::sort(trajPoints.begin(), trajPoints.end(), ltvec);
1150  else
1151  std::sort(trajPoints.begin(), trajPoints.end(), gtvec);
1152 
1153  double deltaZ = 0.;
1154  TVector3 neighbor;
1156  // loop over the trajectory points and set the x & y positions
1157  // be sure to get the closest neighboring trajectory point in the
1158  // opposite view to be able to set the 3D positions correctly
1159  // and use the direction cosines from the 2D tracks
1160  for(size_t tp = 0; tp < trajPoints.size(); ++tp){
1161  view = trajPoints[tp].first;
1162  if(tp < 1){
1163  deltaZ = trajPoints[tp].second.Z() - trajPoints[tp+1].second.Z();
1164  neighbor = trajPoints[tp+1].second;
1165  }
1166  else{
1167  deltaZ = trajPoints[tp].second.Z() - trajPoints[tp-1].second.Z();
1168  neighbor = trajPoints[tp-1].second;
1169  }
1170  if (view == geo::kX) trajPoints[tp].second.SetY(neighbor.Y() + deltaZ*dydz.Y()/dydz.Z());
1171  else if(view == geo::kY) trajPoints[tp].second.SetX(neighbor.X() + deltaZ*dxdz.X()/dxdz.Z());
1172 
1173  //trajPoints[tp].second.Print();
1174  } // end loop over trajectory points
1175 
1176  // make the 3D track
1177  rb::Track retTrack(xTrack.AllCells(),
1178  trajPoints.front().second,
1179  TVector3(trajPoints.back().second.X()-trajPoints.front().second.X(),
1180  trajPoints.back().second.Y()-trajPoints.front().second.Y(),
1181  trajPoints.back().second.Z()-trajPoints.front().second.Z())
1182  );
1183  retTrack.Add(yTrack.AllCells());
1184  for(size_t tp = 1; tp < trajPoints.size(); ++tp) retTrack.AppendTrajectoryPoint(trajPoints[tp].second);
1185 
1186 
1187  this->CheckTrackDirectionInY(retTrack, trajPoints.front().second, trajPoints.back().second);
1188 
1189  // restrict track to exist inside detector; if this isn't possible,
1190  // the track is no good, so don't write
1191  if(geo::ClampRayToDetectorVolume(&trajPoints.front().second,
1192  &trajPoints.back().second,
1193  &(*fGeo))) return retTrack;
1194  //else LOG_DEBUG("WindowTrackAlg") << "Cannot clamp short track to detector";
1195 
1196  }// end if this works
1197 
1198  return retTrack;
1199 
1200  }
rb::Track ShortViewTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, bool const &twoPlanes, int const &firstPlane, int const &secondPlane, geo::View_t const &view)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
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
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
void CheckTrackDirectionInY(rb::Track &track, TVector3 &start, TVector3 &end)
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void ShortTrackRemoveExtraPlane(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > &planeCoG, unsigned short int &removedPlane)
reference front()
Definition: PtrVector.h:379
bool ltvec(std::pair< geo::View_t, TVector3 > a, std::pair< geo::View_t, TVector3 > b)
bool gtvec(std::pair< geo::View_t, TVector3 > a, std::pair< geo::View_t, TVector3 > b)
void ShortTrackExtraPlane(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, std::vector< std::pair< geo::View_t, TVector3 > > &trajPoints, rb::Track &viewTrack, unsigned short int const &removedPlane)
void trk::WindowTrackingAlg::ShortTrackExtraPlane ( std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &  viewHits,
std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &  planeCoG,
std::vector< std::pair< geo::View_t, TVector3 > > &  trajPoints,
rb::Track viewTrack,
unsigned short int const &  removedPlane 
)
private

Definition at line 1289 of file WindowTrackingAlg.cxx.

References std::abs(), rb::Cluster::Add(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), fDHitGood, fGeo, geo::CellGeo::GetCenter(), make_syst_table_plots::h, geo::kX, geo::GeometryBase::Plane(), and rb::Cluster::View().

Referenced by ShortTrack().

1296  {
1297 
1298  // get the expected trajectory point for that view
1299  double xyz[3] = {0.};
1300  fGeo->Plane(removedPlane)->Cell(0)->GetCenter(xyz);
1301  TVector3 trajPoint;
1302  float dirZ = viewTrack.Dir().Z();
1303  float dirV = 0.;
1304  float cog = planeCoG.find(removedPlane)->second.first/planeCoG.find(removedPlane)->second.second;
1305  float expectedV = 0.;
1306  float dV = 0.;
1307  if(viewTrack.View() == geo::kX){
1308  trajPoint = trajPoints[2].second;
1309  dirV = viewTrack.Dir().X();
1310  expectedV = trajPoint.X() + (dirV/dirZ)*std::abs(xyz[2] - trajPoint.Z());
1311  dV = cog - trajPoint.X();
1312 
1313  // set the trajectory point for this plane in case we want to keep it
1314  trajPoint.SetXYZ(cog, 0., xyz[2]);
1315  }
1316  else{
1317  trajPoint = trajPoints[3].second;
1318  dirV = viewTrack.Dir().Y();
1319  expectedV = trajPoint.Y() + (dirV/dirZ)*std::abs(xyz[2] - trajPoint.Z());
1320  dV = cog - trajPoint.Y();
1321 
1322  // set the trajectory point for this plane in case we want to keep it
1323  trajPoint.SetXYZ(0., cog, xyz[2]);
1324  }
1325 
1326  // compare the expected trajectory location to the center of gravity for that plane, and if they are
1327  // within 8*fDHitGood, ie a very steep track, add it back in
1328  if(std::abs(expectedV - cog) < 8*fDHitGood && dV*dirV > 0){
1329  trajPoints.push_back(std::pair<geo::View_t, TVector3>(viewTrack.View(), trajPoint));
1330  for(auto h : viewHits.at(removedPlane)) viewTrack.Add(h);
1331  }
1332 
1333  return;
1334  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const PlaneGeo * Plane(unsigned int i) const
float abs(float number)
Definition: d0nt_math.hpp:39
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void trk::WindowTrackingAlg::ShortTrackRemoveExtraPlane ( std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &  viewHits,
std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > &  planeCoG,
unsigned short int removedPlane 
)
private

Definition at line 1205 of file WindowTrackingAlg.cxx.

Referenced by ShortTrack().

1210  {
1211  unsigned short int minPlane = 0;
1212 
1213  size_t minHits = 1e6;
1214  for(auto vh : viewHits){
1215  if(vh.second.size() < minHits){
1216  minPlane = vh.first;
1217  minHits = vh.second.size();
1218  }
1219  }
1220 
1221  for(auto cogItr : planeCoG){
1222  if(cogItr.first == minPlane){
1223  removedPlane = cogItr.first;
1224  viewHits.erase(removedPlane);
1225  planeCoG.erase(removedPlane);
1226  }
1227  }
1228 
1229  return;
1230  }
rb::Track trk::WindowTrackingAlg::ShortViewTrack ( std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &  viewHits,
std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &  planeCoG,
bool const &  twoPlanes,
int const &  firstPlane,
int const &  secondPlane,
geo::View_t const &  view 
)
private

Definition at line 1234 of file WindowTrackingAlg.cxx.

References rb::Track::AppendTrajectoryPoint(), geo::PlaneGeo::Cell(), fCellHalfDepth, fGeo, hits(), geo::kX, geo::GeometryBase::Plane(), and track.

Referenced by ShortTrack().

1243  {
1244  // make 2d tracks using the center of gravity for each plane
1245  std::vector<art::Ptr<rb::CellHit> > hits;
1246  for(auto itr : viewHits){
1247  for(auto ch : itr.second) hits.push_back(ch);
1248  }
1249 
1250  double xyz1[3] = {0.};
1251  double xyz2[3] = {0.};
1252  auto chVecItr = viewHits.cbegin()->second.cbegin();
1253  fGeo->Plane(firstPlane) ->Cell((*chVecItr)->Cell())->GetCenter(xyz1);
1254  float cog1 = planeCoG.find(firstPlane) ->second.first/planeCoG.find(firstPlane) ->second.second;
1255  float deltaW = 0.;
1256  float deltaZ = 2.*fCellHalfDepth;
1257  if(twoPlanes){
1258  chVecItr = viewHits.crbegin()->second.cbegin();
1259  fGeo->Plane(secondPlane)->Cell((*chVecItr)->Cell())->GetCenter(xyz2);
1260  deltaW = (planeCoG.find(secondPlane)->second.first/planeCoG.find(secondPlane)->second.second);
1261  deltaZ = xyz2[2]-xyz1[2];
1262  }
1263  else{
1264  // hits only on one plane, so use the center of the first and last cells in the plane to find the
1265  // spread in the current view
1266  chVecItr = viewHits.cbegin()->second.cbegin();
1267  fGeo->Plane(firstPlane)->Cell((*chVecItr)->Cell())->GetCenter(xyz1);
1268  cog1 = (view == geo::kX) ? xyz1[0] : xyz1[1];
1269  fGeo->Plane(firstPlane)->Cell((viewHits.cbegin()->second.back())->Cell())->GetCenter(xyz2);
1270  deltaW = (view == geo::kX) ? xyz2[0] : xyz2[1];
1271  }
1272  deltaW -= cog1;
1273 
1274  // be sure we don't just have one cell in the view
1275  if(deltaW == 0.) deltaW = fCellHalfDepth;
1276 
1277 // LOG_VERBATIM("WindowTrackingAlg") << view << " " << cog1 << " " << xyz1[2] << " " << deltaW << " " << deltaZ;
1278 
1279  rb::Track track(hits, view, cog1, xyz1[2], deltaW, deltaZ);
1280  track.AppendTrajectoryPoint(deltaW+cog1, deltaZ);
1281 
1282  return track;
1283  }
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
void hits()
Definition: readHits.C:15
double fCellHalfDepth
half of cell size in z
Float_t track
Definition: plot.C:35
bool trk::WindowTrackingAlg::ValidTrack ( std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &  planeToHit)
private

Definition at line 1364 of file WindowTrackingAlg.cxx.

References std::abs(), fMinHitFraction, and fMinViewPlanes.

Referenced by Make3DTrack(), and MakeViewTrack().

1366  {
1367  if(planeToHit.size() < fMinViewPlanes) return false;
1368 
1369  unsigned int const& startPlane = planeToHit.cbegin()->first;
1370  unsigned int const& endPlane = planeToHit.crbegin()->first;
1371 
1372  // the 0.5 in the MinHitFraction accounts for the fact that
1373  // every other plane is X or Y as we are using the full extent
1374  // of the slice
1375  if(planeToHit.size() < fMinHitFraction*0.5*std::abs(1.*startPlane - 1.*endPlane) )
1376  return false;
1377 
1378  return true;
1379  }
float abs(float number)
Definition: d0nt_math.hpp:39
unsigned int fMinViewPlanes
Minimum number of planes in a view to fit.
double fMinHitFraction
Minimum value of hits/plane to look for a track.

Member Data Documentation

double trk::WindowTrackingAlg::fCellHalfDepth
private

half of cell size in z

Definition at line 118 of file WindowTrackingAlg.h.

Referenced by MakeTrack(), and ShortViewTrack().

double trk::WindowTrackingAlg::fDHitGood
private

Maximum distance from the line for a hit to be considered part of it

Definition at line 119 of file WindowTrackingAlg.h.

Referenced by FitWindow(), LookForBremsstrahlungHits(), MakeTrack(), ShortTrackExtraPlane(), and WindowTrackingAlg().

art::ServiceHandle<geo::Geometry> trk::WindowTrackingAlg::fGeo
private
int trk::WindowTrackingAlg::fMaxPlaneSeparation
private

Maximum number of planes allowed to separate planes with hits in slice.

Definition at line 126 of file WindowTrackingAlg.h.

Referenced by Make3DTrack(), and WindowTrackingAlg().

double trk::WindowTrackingAlg::fMinHitFraction
private

Minimum value of hits/plane to look for a track.

Definition at line 124 of file WindowTrackingAlg.h.

Referenced by ValidTrack(), and WindowTrackingAlg().

unsigned int trk::WindowTrackingAlg::fMinHitsTryAgain
private

Minimum number of hits left to look for another track.

Definition at line 123 of file WindowTrackingAlg.h.

Referenced by MakeTrack(), and WindowTrackingAlg().

unsigned int trk::WindowTrackingAlg::fMinViewPlanes
private

Minimum number of planes in a view to fit.

Definition at line 122 of file WindowTrackingAlg.h.

Referenced by Make3DTrack(), ValidTrack(), and WindowTrackingAlg().

unsigned int trk::WindowTrackingAlg::fNTrajectory
private

number of trajectory points to use in determining incoming direction

Definition at line 125 of file WindowTrackingAlg.h.

Referenced by DetermineInitialDirection(), and WindowTrackingAlg().

unsigned int trk::WindowTrackingAlg::fWindowSize
private

Number of planes in the sliding window.

Definition at line 121 of file WindowTrackingAlg.h.

Referenced by FitView(), and WindowTrackingAlg().


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