LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
pma::ProjectionMatchingAlg Class Reference

#include "ProjectionMatchingAlg.h"

Classes

struct  Config
 

Public Member Functions

 ProjectionMatchingAlg (const Config &config)
 
 ProjectionMatchingAlg (const fhicl::ParameterSet &pset)
 
double validate_on_adc (const detinfo::DetectorPropertiesData &detProp, const lariov::ChannelStatusProvider &channelStatus, const pma::Track3D &trk, const img::DataProviderAlg &adcImage, float thr) const
 
double validate_on_adc_test (const detinfo::DetectorPropertiesData &detProp, const lariov::ChannelStatusProvider &channelStatus, const pma::Track3D &trk, const img::DataProviderAlg &adcImage, const std::vector< art::Ptr< recob::Hit >> &hits, TH1F *histoPassing, TH1F *histoRejected) const
 
double validate (const detinfo::DetectorPropertiesData &detProp, const lariov::ChannelStatusProvider &channelStatus, const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit >> &hits) const
 
double validate (const detinfo::DetectorPropertiesData &detProp, const lariov::ChannelStatusProvider &channelStatus, const TVector3 &p0, const TVector3 &p1, const std::vector< art::Ptr< recob::Hit >> &hits, unsigned int testView, unsigned int tpc, unsigned int cryo) const
 
double twoViewFraction (pma::Track3D &trk) const
 
unsigned int testHits (detinfo::DetectorPropertiesData const &detProp, const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit >> &hits, double eps=1.0) const
 Count the number of hits that are closer than eps * fHitTestingDist2D to the track 2D projection. More...
 
bool isContained (const pma::Track3D &trk, float margin=0.0F) const
 
pma::Track3DbuildTrack (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
 
pma::Track3DbuildMultiTPCTrack (const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits) const
 
pma::Track3DbuildShowerSeg (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, const pma::Vector3D &vtx) const
 
pma::Track3DbuildSegment (const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
 
pma::Track3DbuildSegment (const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2, const geo::Point_t &point) const
 
pma::Track3DbuildSegment (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, const TVector3 &point) const
 
void FilterOutSmallParts (const detinfo::DetectorPropertiesData &detProp, double r2d, const std::vector< art::Ptr< recob::Hit >> &hits_in, std::vector< art::Ptr< recob::Hit >> &hits_out, const TVector2 &vtx2d) const
 
void RemoveNotEnabledHits (pma::Track3D &trk) const
 
pma::Track3DextendTrack (const detinfo::DetectorPropertiesData &clockData, const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit >> &hits, bool add_nodes) const
 Add more hits to an existing track, reoptimize, optionally add more nodes. More...
 
void guideEndpoints (const detinfo::DetectorPropertiesData &clockData, pma::Track3D &trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &hits) const
 
void guideEndpoints (const detinfo::DetectorPropertiesData &clockData, pma::Track3D &trk, pma::Track3D::ETrackEnd endpoint, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &hits) const
 Add 3D reference points to clean endpoint of a track. More...
 
std::vector< pma::Hit3D * > trimTrackToVolume (pma::Track3D &trk, TVector3 p0, TVector3 p1) const
 
bool alignTracks (pma::Track3D &first, pma::Track3D &second) const
 
void mergeTracks (const detinfo::DetectorPropertiesData &detProp, pma::Track3D &dst, pma::Track3D &src, bool reopt) const
 
void autoFlip (pma::Track3D &trk, pma::Track3D::EDirection dir=Track3D::kForward, double thr=0.0, unsigned int n=0) const
 
double selectInitialHits (pma::Track3D &trk, unsigned int view=geo::kZ, unsigned int *nused=0) const
 

Private Member Functions

bool chkEndpointHits_ (const detinfo::DetectorPropertiesData &detProp, int wire, int wdir, double drift_x, int view, unsigned int tpc, unsigned int cryo, const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit >> &hits) const
 
bool addEndpointRef_ (const detinfo::DetectorPropertiesData &detProp, pma::Track3D &trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &hits, std::pair< int, int > const *wires, double const *xPos, unsigned int tpc, unsigned int cryo) const
 
bool GetCloseHits_ (const detinfo::DetectorPropertiesData &detProp, double r2d, const std::vector< art::Ptr< recob::Hit >> &hits_in, std::vector< size_t > &used, std::vector< art::Ptr< recob::Hit >> &hits_out) const
 
bool Has_ (const std::vector< size_t > &v, size_t idx) const
 
void ShortenSeg_ (const detinfo::DetectorPropertiesData &detProp, pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
 
bool TestTrk_ (pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
 

Static Private Member Functions

static size_t getSegCount_ (size_t trk_size)
 

Private Attributes

double const fOptimizationEps
 
double const fFineTuningEps
 
double const fTrkValidationDist2D
 
double const fHitTestingDist2D
 
double const fMinTwoViewFraction
 
geo::GeometryCore const * fGeom
 

Detailed Description

Definition at line 62 of file ProjectionMatchingAlg.h.

Constructor & Destructor Documentation

pma::ProjectionMatchingAlg::ProjectionMatchingAlg ( const Config config)

Definition at line 31 of file ProjectionMatchingAlg.cxx.

References fFineTuningEps, fGeom, fHitTestingDist2D, fMinTwoViewFraction, fTrkValidationDist2D, geo::kU, geo::kV, geo::kZ, pma::ProjectionMatchingAlg::Config::OptimizationEps, pma::Node3D::SetMargin(), and pma::Element3D::SetOptFactor().

32  : fOptimizationEps{config.OptimizationEps()}
33  , fFineTuningEps{config.FineTuningEps()}
34  , fTrkValidationDist2D{config.TrkValidationDist2D()}
35  , fHitTestingDist2D{config.HitTestingDist2D()}
36  , fMinTwoViewFraction{config.MinTwoViewFraction()}
37  , fGeom{lar::providerFrom<geo::Geometry>()}
38 {
39  pma::Node3D::SetMargin(config.NodeMargin3D());
40 
41  pma::Element3D::SetOptFactor(geo::kU, config.HitWeightU());
42  pma::Element3D::SetOptFactor(geo::kV, config.HitWeightV());
43  pma::Element3D::SetOptFactor(geo::kZ, config.HitWeightZ());
44 }
static void SetOptFactor(unsigned int view, float value)
Definition: PmaElement3D.h:113
geo::GeometryCore const * fGeom
Planes which measure V.
Definition: geo_types.h:136
Planes which measure Z direction.
Definition: geo_types.h:138
static void SetMargin(double m)
Set allowed node position margin around TPC.
Definition: PmaNode3D.h:124
Planes which measure U.
Definition: geo_types.h:135
pma::ProjectionMatchingAlg::ProjectionMatchingAlg ( const fhicl::ParameterSet pset)
inline

Definition at line 102 of file ProjectionMatchingAlg.h.

References hits(), and lar::dump::vector().

104  {}
ProjectionMatchingAlg(const Config &config)

Member Function Documentation

bool pma::ProjectionMatchingAlg::addEndpointRef_ ( const detinfo::DetectorPropertiesData detProp,
pma::Track3D trk,
const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &  hits,
std::pair< int, int > const *  wires,
double const *  xPos,
unsigned int  tpc,
unsigned int  cryo 
) const
private

Definition at line 1031 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddRefPoint(), chkEndpointHits_(), fGeom, hits(), geo::GeometryCore::IntersectionPoint(), pma::Track3D::size(), x, y, and z.

Referenced by guideEndpoints().

1039 {
1040  double x = 0.0, y = 0.0, z = 0.0;
1041  std::vector<std::pair<int, unsigned int>> wire_view;
1042  for (unsigned int i = 0; i < 3; i++)
1043  if (wires[i].first >= 0) {
1044  const auto hiter = hits.find(i);
1045  if (hiter != hits.end()) {
1046  if (chkEndpointHits_(detProp,
1047  wires[i].first,
1048  wires[i].second,
1049  xPos[i],
1050  i,
1051  tpc,
1052  cryo,
1053  trk,
1054  hiter->second)) {
1055  x += xPos[i];
1056  wire_view.emplace_back(wires[i].first, i);
1057  }
1058  }
1059  }
1060  if (wire_view.size() > 1) {
1061  x /= wire_view.size();
1062  auto const [wire0, plane0] = wire_view[0];
1063  auto const [wire1, plane1] = wire_view[1];
1065  geo::WireID(cryo, tpc, plane0, wire0), geo::WireID(cryo, tpc, plane1, wire1), y, z);
1066  trk.AddRefPoint(x, y, z);
1067  mf::LogVerbatim("ProjectionMatchingAlg")
1068  << "trk tpc:" << tpc << " size:" << trk.size() << " add ref.point (" << x << "; " << y << "; "
1069  << z << ")";
1070  return true;
1071  }
1072 
1073  mf::LogVerbatim("ProjectionMatchingAlg")
1074  << "trk tpc:" << tpc << " size:" << trk.size()
1075  << " wire-plane-parallel track, but need two clean views of endpoint";
1076  return false;
1077 }
Float_t x
Definition: compare.C:6
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
geo::GeometryCore const * fGeom
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:276
void AddRefPoint(const TVector3 &p)
Definition: PmaTrack3D.h:218
bool chkEndpointHits_(const detinfo::DetectorPropertiesData &detProp, int wire, int wdir, double drift_x, int view, unsigned int tpc, unsigned int cryo, const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit >> &hits) const
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
size_t size() const
Definition: PmaTrack3D.h:89
bool pma::ProjectionMatchingAlg::alignTracks ( pma::Track3D first,
pma::Track3D second 
) const

Flip tracks to get second as a continuation of first; returns false if not possible (tracks in reversed order).

Definition at line 1275 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), larg4::dist(), pma::Dist2(), pma::Track3D::Flip(), pma::Track3D::front(), and pma::Hit3D::Point3D().

Referenced by mergeTracks().

1276 {
1277  unsigned int k = 0;
1278  double const distFF = pma::Dist2(first.front()->Point3D(), second.front()->Point3D());
1279  double dist = distFF;
1280 
1281  double const distFB = pma::Dist2(first.front()->Point3D(), second.back()->Point3D());
1282  if (distFB < dist) {
1283  k = 1;
1284  dist = distFB;
1285  }
1286 
1287  double const distBF = pma::Dist2(first.back()->Point3D(), second.front()->Point3D());
1288  if (distBF < dist) {
1289  k = 2;
1290  dist = distBF;
1291  }
1292 
1293  double distBB = pma::Dist2(first.back()->Point3D(), second.back()->Point3D());
1294  if (distBB < dist) {
1295  k = 3;
1296  dist = distBB;
1297  }
1298 
1299  switch (k) // flip to get dst end before src start, do not merge if track's order reversed
1300  {
1301  case 0:
1302  first.Flip(); // detProp);
1303  break;
1304  case 1: mf::LogError("PMAlgTrackMaker") << "Tracks in reversed order."; return false;
1305  case 2: break;
1306  case 3:
1307  second.Flip(); // detProp);
1308  break;
1309  default:
1310  throw cet::exception("pma::ProjectionMatchingAlg")
1311  << "alignTracks: should never happen." << std::endl;
1312  }
1313  return true;
1314 }
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:88
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:527
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void pma::ProjectionMatchingAlg::autoFlip ( pma::Track3D trk,
pma::Track3D::EDirection  dir = Track3D::kForward,
double  thr = 0.0,
unsigned int  n = 0 
) const
inline

Try to correct track direction of the stopping particle: dir: kForward - particle stop is at the end of the track; kBackward - particle stop is at the beginning of the track; dQ/dx difference has to be above thr to actually flip the track; compares dQ/dx of n hits at each end of the track (default is based on the track length).

Definition at line 264 of file ProjectionMatchingAlg.h.

References pma::Track3D::AutoFlip(), dir, hits(), geo::kZ, n, and lar::dump::vector().

268  {
269  trk.AutoFlip(dir, thr, n);
270  };
void AutoFlip(pma::Track3D::EDirection dir, double thr=0.0, unsigned int n=0)
Definition: PmaTrack3D.cxx:662
TDirectory * dir
Definition: macro.C:5
Char_t n[5]
pma::Track3D * pma::ProjectionMatchingAlg::buildMultiTPCTrack ( const detinfo::DetectorPropertiesData clockData,
const std::vector< art::Ptr< recob::Hit >> &  hits 
) const

Build a track from sets of hits, multiple TPCs are OK (like taken from PFParticles), as far as hits origin from at least two wire planes.

Definition at line 499 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), buildSegment(), d, pma::Dist2(), fFineTuningEps, fOptimizationEps, pma::Track3D::front(), getSegCount_(), hits(), mergeTracks(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), pma::Track3D::Segments(), pma::Track3D::SelectAllHits(), pma::Track3D::size(), and pma::Track3D::SortHits().

Referenced by pma::PMAlgFitter::buildTracks().

502 {
503  std::map<unsigned int, std::vector<art::Ptr<recob::Hit>>> hits_by_tpc;
504  for (auto const& h : hits) {
505  hits_by_tpc[h->WireID().TPC].push_back(h);
506  }
507 
508  std::vector<pma::Track3D*> tracks;
509  for (auto const& hsel : hits_by_tpc) {
510  pma::Track3D* trk = buildSegment(detProp, hsel.second);
511  if (trk) tracks.push_back(trk);
512  }
513 
514  bool need_reopt = false;
515  while (tracks.size() > 1) {
516  need_reopt = true;
517 
518  pma::Track3D* first = tracks.front();
519  pma::Track3D* best = 0;
520  double d, dmin = 1.0e12;
521  size_t t_best = 0, cfg = 0;
522  for (size_t t = 1; t < tracks.size(); ++t) {
523  pma::Track3D* second = tracks[t];
524 
525  d = pma::Dist2(first->front()->Point3D(), second->front()->Point3D());
526  if (d < dmin) {
527  dmin = d;
528  best = second;
529  t_best = t;
530  cfg = 0;
531  }
532 
533  d = pma::Dist2(first->front()->Point3D(), second->back()->Point3D());
534  if (d < dmin) {
535  dmin = d;
536  best = second;
537  t_best = t;
538  cfg = 1;
539  }
540 
541  d = pma::Dist2(first->back()->Point3D(), second->front()->Point3D());
542  if (d < dmin) {
543  dmin = d;
544  best = second;
545  t_best = t;
546  cfg = 2;
547  }
548 
549  d = pma::Dist2(first->back()->Point3D(), second->back()->Point3D());
550  if (d < dmin) {
551  dmin = d;
552  best = second;
553  t_best = t;
554  cfg = 3;
555  }
556  }
557  if (best) {
558  switch (cfg) {
559  default:
560  case 0:
561  case 1:
562  mergeTracks(detProp, *best, *first, false);
563  tracks[0] = best;
564  delete first;
565  break;
566 
567  case 2:
568  case 3:
569  mergeTracks(detProp, *first, *best, false);
570  delete best;
571  break;
572  }
573  tracks.erase(tracks.begin() + t_best);
574  }
575  else
576  break; // should not happen
577  }
578 
579  pma::Track3D* trk = 0;
580  if (!tracks.empty()) {
581  trk = tracks.front();
582  if (need_reopt) {
583  double g = trk->Optimize(detProp, 0, fOptimizationEps);
584  mf::LogVerbatim("ProjectionMatchingAlg")
585  << " reopt after merging initial tpc segments: done, g = " << g;
586  }
587 
588  int nSegments = getSegCount_(trk->size()) - trk->Segments().size() + 1;
589  if (nSegments > 0) // need to add segments
590  {
591  double g = 0.0;
592  size_t nNodes = (size_t)(nSegments - 1); // n nodes to add
593  if (nNodes) {
594  mf::LogVerbatim("ProjectionMatchingAlg") << " optimize trk (add " << nSegments << " seg)";
595 
596  g = trk->Optimize(detProp, nNodes, fOptimizationEps, false, true, 25,
597  10); // build nodes
598  mf::LogVerbatim("ProjectionMatchingAlg") << " nodes done, g = " << g;
599  trk->SelectAllHits();
600  }
601  g = trk->Optimize(detProp, 0, fFineTuningEps); // final tuning
602  mf::LogVerbatim("ProjectionMatchingAlg") << " tune done, g = " << g;
603  }
604  trk->SortHits();
605  }
606  return trk;
607 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool SelectAllHits()
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
static size_t getSegCount_(size_t trk_size)
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:269
Float_t d
Definition: plot.C:235
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
void mergeTracks(const detinfo::DetectorPropertiesData &detProp, pma::Track3D &dst, pma::Track3D &src, bool reopt) const
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:88
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
size_t size() const
Definition: PmaTrack3D.h:89
pma::Track3D * pma::ProjectionMatchingAlg::buildSegment ( const detinfo::DetectorPropertiesData clockData,
const std::vector< art::Ptr< recob::Hit >> &  hits_1,
const std::vector< art::Ptr< recob::Hit >> &  hits_2 = {} 
) const

Build a straight segment from two sets of hits (they should origin from two wire planes); method is intendet for short tracks or shower initial parts, where only a few hits per plane are available and there is no chance to see a curvature or any other features.

Definition at line 889 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddHits(), fFineTuningEps, pma::Track3D::HasTwoViews(), pma::Track3D::Initialize(), pma::Track3D::Optimize(), pma::Track3D::SetEndSegWeight(), pma::Track3D::SortHits(), and pma::Track3D::TPCs().

Referenced by buildMultiTPCTrack(), buildSegment(), buildShowerSeg(), shower::EMShowerAlg::ConstructTrack(), ems::EMShower3D::Make3DSeg(), shower::EMShowerAlg::MakeShower(), ems::EMShower3D::Reoptimize(), and ems::EMShower3D::Validate().

893 {
894  pma::Track3D* trk = new pma::Track3D();
895  trk->SetEndSegWeight(0.001F);
896  trk->AddHits(detProp, hits_1);
897  trk->AddHits(detProp, hits_2);
898 
899  if (trk->HasTwoViews() && (trk->TPCs().size() == 1)) // now only in single tpc
900  {
901  if (!trk->Initialize(detProp, 0.001F)) {
902  mf::LogWarning("ProjectionMatchingAlg") << "initialization failed, delete segment";
903  delete trk;
904  return 0;
905  }
906  trk->Optimize(detProp, 0, fFineTuningEps);
907 
908  trk->SortHits();
909  return trk;
910  }
911  else {
912  mf::LogWarning("ProjectionMatchingAlg") << "need at least two views in single tpc";
913  delete trk;
914  return 0;
915  }
916 }
void AddHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
Definition: PmaTrack3D.cxx:404
bool Initialize(detinfo::DetectorPropertiesData const &detProp, float initEndSegW=0.05F)
Definition: PmaTrack3D.cxx:77
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:448
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
void SetEndSegWeight(float value) noexcept
Definition: PmaTrack3D.h:320
bool HasTwoViews(size_t nmin=1) const
Definition: PmaTrack3D.cxx:439
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
pma::Track3D * pma::ProjectionMatchingAlg::buildSegment ( const detinfo::DetectorPropertiesData clockData,
const std::vector< art::Ptr< recob::Hit >> &  hits_1,
const std::vector< art::Ptr< recob::Hit >> &  hits_2,
const geo::Point_t point 
) const

Build a straight segment from two sets of hits (they should origin from two wire planes), starting from a given point (like vertex known from another algorithm); method is intendet for short tracks or shower initial parts, where only a few hits per plane are available and there is no chance to see a curvature or any other features.

Definition at line 919 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), buildSegment(), pma::Dist2(), fFineTuningEps, pma::Track3D::Flip(), pma::Track3D::front(), pma::Track3D::Nodes(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), and pma::Track3D::SortHits().

924 {
925  pma::Track3D* trk = buildSegment(detProp, hits_1, hits_2);
926 
927  if (trk) {
928  double dfront = pma::Dist2(trk->front()->Point3D(), point);
929  double dback = pma::Dist2(trk->back()->Point3D(), point);
930  if (dfront > dback) trk->Flip();
931 
932  trk->Nodes().front()->SetPoint3D({point.X(), point.Y(), point.Z()});
933  trk->Nodes().front()->SetFrozen(true);
934  trk->Optimize(detProp, 0, fFineTuningEps);
935 
936  trk->SortHits();
937  }
938  return trk;
939 }
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:88
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:527
pma::Track3D * pma::ProjectionMatchingAlg::buildSegment ( const detinfo::DetectorPropertiesData detProp,
const std::vector< art::Ptr< recob::Hit >> &  hits,
const TVector3 &  point 
) const

Build a straight segment from set of hits (they should origin from two wire planes at least), starting from a given point.

Definition at line 942 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), buildSegment(), pma::Dist2(), fFineTuningEps, pma::Track3D::Flip(), pma::Track3D::front(), hits(), pma::Track3D::Nodes(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), and pma::Track3D::SortHits().

946 {
947  pma::Track3D* trk = buildSegment(detProp, hits);
948 
949  if (trk) {
950  double dfront = pma::Dist2(trk->front()->Point3D(), point);
951  double dback = pma::Dist2(trk->back()->Point3D(), point);
952  if (dfront > dback) trk->Flip(); // detProp);
953 
954  trk->Nodes().front()->SetPoint3D(point);
955  trk->Nodes().front()->SetFrozen(true);
956  trk->Optimize(detProp, 0, fFineTuningEps);
957 
958  trk->SortHits();
959  }
960  return trk;
961 }
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:88
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:527
pma::Track3D * pma::ProjectionMatchingAlg::buildShowerSeg ( const detinfo::DetectorPropertiesData detProp,
const std::vector< art::Ptr< recob::Hit >> &  hits,
const pma::Vector3D vtx 
) const

Build a shower segment from sets of hits and attached to the provided vertex.

Definition at line 611 of file ProjectionMatchingAlg.cxx.

References buildSegment(), geo::CryostatID::Cryostat, fGeom, FilterOutSmallParts(), geo::GeometryCore::FindTPCAtPosition(), pma::GetProjectionToPlane(), geo::TPCGeo::HasPlane(), geo::GeometryCore::HasTPC(), hits(), geo::GeometryCore::PositionToCryostatID(), ShortenSeg_(), pma::Track3D::size(), geo::TPCID::TPC, and geo::GeometryCore::TPC().

Referenced by pma::PMAlgFitter::buildShowers().

615 {
616  geo::Point_t const vtxpoint{vtx.X(), vtx.Y(), vtx.Z()};
617 
618  if (!fGeom->HasTPC(fGeom->FindTPCAtPosition(vtxpoint))) return 0;
619 
620  TVector3 vtxv3(vtx.X(), vtx.Y(), vtx.Z());
621 
622  const size_t tpc = fGeom->FindTPCAtPosition(vtxpoint).TPC;
623  const size_t cryo = fGeom->PositionToCryostatID(vtxpoint).Cryostat;
624  const geo::TPCGeo& tpcgeom = fGeom->TPC(geo::TPCID(cryo, tpc));
625 
626  // use only hits from tpc where the vtx is
627  std::vector<art::Ptr<recob::Hit>> hitstpc;
628  for (size_t h = 0; h < hits.size(); ++h)
629  if (hits[h]->WireID().TPC == tpc) hitstpc.push_back(hits[h]);
630 
631  if (!hitstpc.size()) return 0;
632 
633  std::vector<art::Ptr<recob::Hit>> hitstrk;
634  size_t view = 0;
635  size_t countviews = 0;
636  while (view < 3) {
637  mf::LogVerbatim("ProjectionMatchinAlg") << "collecting hits from view: " << view;
638  if (!tpcgeom.HasPlane(view)) {
639  ++view;
640  continue;
641  }
642 
643  // select hits only for a single view
644  std::vector<art::Ptr<recob::Hit>> hitsview;
645  for (size_t h = 0; h < hitstpc.size(); ++h)
646  if (hitstpc[h]->WireID().Plane == view) hitsview.push_back(hitstpc[h]);
647  if (!hitsview.size()) {
648  ++view;
649  continue;
650  }
651 
652  // filter our small groups of hits, far from main shower
653  std::vector<art::Ptr<recob::Hit>> hitsfilter;
654  TVector2 proj_pr = pma::GetProjectionToPlane(vtxv3, view, tpc, cryo);
655 
656  mf::LogVerbatim("ProjectionMatchinAlg") << "start filter out: ";
657  FilterOutSmallParts(detProp, 2.0, hitsview, hitsfilter, proj_pr);
658  mf::LogVerbatim("ProjectionMatchingAlg") << "after filter out";
659 
660  for (size_t h = 0; h < hitsfilter.size(); ++h)
661  hitstrk.push_back(hitsfilter[h]);
662 
663  if (hitsfilter.size() >= 5) countviews++;
664 
665  ++view;
666  }
667 
668  if (!hitstrk.size() || (countviews < 2)) {
669  mf::LogWarning("ProjectionMatchinAlg") << "too few hits, segment not built";
670  return 0;
671  }
672 
673  // hits are prepared, finally segment is built
674 
675  pma::Track3D* trk = new pma::Track3D();
676  trk = buildSegment(detProp, hitstrk, vtxv3);
677 
678  // make shorter segment to estimate direction more precise
679  ShortenSeg_(detProp, *trk, tpcgeom);
680 
681  if (trk->size() < 10) {
682  mf::LogWarning("ProjectionMatchingAlg") << "too short segment, delete segment";
683  delete trk;
684  return 0;
685  }
686 
687  return trk;
688 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:147
geo::GeometryCore const * fGeom
Geometry information for a single TPC.
Definition: TPCGeo.h:36
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void FilterOutSmallParts(const detinfo::DetectorPropertiesData &detProp, double r2d, const std::vector< art::Ptr< recob::Hit >> &hits_in, std::vector< art::Ptr< recob::Hit >> &hits_out, const TVector2 &vtx2d) const
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
CryostatID PositionToCryostatID(Point_t const &point) const
Returns the ID of the cryostat at specified location.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
void ShortenSeg_(const detinfo::DetectorPropertiesData &detProp, pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:263
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
size_t size() const
Definition: PmaTrack3D.h:89
bool HasTPC(TPCID const &tpcid) const
Returns whether we have the specified TPC.
Definition: GeometryCore.h:699
pma::Track3D * pma::ProjectionMatchingAlg::buildTrack ( const detinfo::DetectorPropertiesData detProp,
const std::vector< art::Ptr< recob::Hit >> &  hits_1,
const std::vector< art::Ptr< recob::Hit >> &  hits_2 = {} 
) const

Build a track from two sets of hits from single TPC, hits should origin from at least two wire planes; number of segments used to create the track depends on the number of hits.

Definition at line 446 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddHits(), f, fFineTuningEps, fMinTwoViewFraction, fOptimizationEps, getSegCount_(), pma::Track3D::Initialize(), geo::kU, geo::kV, geo::kZ, pma::Track3D::NHits(), pma::Track3D::Optimize(), pma::Track3D::SelectAllHits(), pma::Track3D::size(), pma::Track3D::SortHits(), pma::Track3D::TPCs(), and twoViewFraction().

Referenced by pma::PMAlgTracker::matchCluster().

450 {
451  pma::Track3D* trk = new pma::Track3D(); // track candidate
452  trk->AddHits(detProp, hits_1);
453  trk->AddHits(detProp, hits_2);
454 
455  mf::LogVerbatim("ProjectionMatchingAlg") << "track size: " << trk->size();
456  std::vector<unsigned int> tpcs = trk->TPCs();
457  for (size_t t = 0; t < tpcs.size(); ++t) {
458  mf::LogVerbatim("ProjectionMatchingAlg") << " tpc:" << tpcs[t];
459  }
460  mf::LogVerbatim("ProjectionMatchingAlg")
461  << " #coll:" << trk->NHits(geo::kZ) << " #ind2:" << trk->NHits(geo::kV)
462  << " #ind1:" << trk->NHits(geo::kU);
463 
464  size_t nSegments = getSegCount_(trk->size());
465  size_t nNodes = (size_t)(nSegments - 1); // n nodes to add
466 
467  mf::LogVerbatim("ProjectionMatchingAlg") << " initialize trk";
468  if (!trk->Initialize(detProp)) {
469  mf::LogWarning("ProjectionMatchingAlg") << " initialization failed, delete trk";
470  delete trk;
471  return 0;
472  }
473 
474  double f = twoViewFraction(*trk);
475  if (f > fMinTwoViewFraction) {
476  double g = 0.0;
477  mf::LogVerbatim("ProjectionMatchingAlg") << " optimize trk (" << nSegments << " seg)";
478  if (nNodes) {
479  g = trk->Optimize(detProp, nNodes, fOptimizationEps, false, true, 25,
480  10); // build nodes
481  mf::LogVerbatim("ProjectionMatchingAlg") << " nodes done, g = " << g;
482  trk->SelectAllHits();
483  }
484  g = trk->Optimize(detProp, 0, fFineTuningEps); // final tuning
485  mf::LogVerbatim("ProjectionMatchingAlg") << " tune done, g = " << g;
486 
487  trk->SortHits();
488  // trk->ShiftEndsToHits(); // not sure if useful already here
489  return trk;
490  }
491  else {
492  mf::LogVerbatim("ProjectionMatchingAlg") << " clusters do not match, f = " << f;
493  delete trk;
494  return 0;
495  }
496 }
void AddHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
Definition: PmaTrack3D.cxx:404
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool SelectAllHits()
Planes which measure V.
Definition: geo_types.h:136
bool Initialize(detinfo::DetectorPropertiesData const &detProp, float initEndSegW=0.05F)
Definition: PmaTrack3D.cxx:77
Planes which measure Z direction.
Definition: geo_types.h:138
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:448
static size_t getSegCount_(size_t trk_size)
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
TFile f
Definition: plotHisto.C:6
Planes which measure U.
Definition: geo_types.h:135
unsigned int NHits(unsigned int view) const
Definition: PmaTrack3D.cxx:421
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
size_t size() const
Definition: PmaTrack3D.h:89
double twoViewFraction(pma::Track3D &trk) const
bool pma::ProjectionMatchingAlg::chkEndpointHits_ ( const detinfo::DetectorPropertiesData detProp,
int  wire,
int  wdir,
double  drift_x,
int  view,
unsigned int  tpc,
unsigned int  cryo,
const pma::Track3D trk,
const std::vector< art::Ptr< recob::Hit >> &  hits 
) const
private

Definition at line 994 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), hits(), pma::Track3D::size(), and x.

Referenced by addEndpointRef_().

1004 {
1005  size_t nCloseHits = 0;
1006  int forwWires = 3, backWires = -1;
1007  double xMargin = 0.4;
1008  for (auto h : hits)
1009  if ((view == (int)h->WireID().Plane) && (tpc == h->WireID().TPC) &&
1010  (cryo == h->WireID().Cryostat)) {
1011  bool found = false;
1012  for (size_t ht = 0; ht < trk.size(); ht++)
1013  if (trk[ht]->Hit2DPtr().key() == h.key()) {
1014  found = true;
1015  break;
1016  }
1017  if (found) continue;
1018 
1019  int dw = wdir * (h->WireID().Wire - wire);
1020  if ((dw <= forwWires) && (dw >= backWires)) {
1021  double x = detProp.ConvertTicksToX(h->PeakTime(), view, tpc, cryo);
1022  if (fabs(x - drift_x) < xMargin) nCloseHits++;
1023  }
1024  }
1025  if (nCloseHits > 1)
1026  return false;
1027  else
1028  return true;
1029 }
Float_t x
Definition: compare.C:6
double ConvertTicksToX(double ticks, int p, int t, int c) const
size_t size() const
Definition: PmaTrack3D.h:89
pma::Track3D * pma::ProjectionMatchingAlg::extendTrack ( const detinfo::DetectorPropertiesData clockData,
const pma::Track3D trk,
const std::vector< art::Ptr< recob::Hit >> &  hits,
bool  add_nodes 
) const

Add more hits to an existing track, reoptimize, optionally add more nodes.

Definition at line 964 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddHits(), fFineTuningEps, fOptimizationEps, getSegCount_(), hits(), geo::kU, geo::kV, geo::kZ, pma::Track3D::NHits(), pma::Track3D::Nodes(), pma::Track3D::Optimize(), and pma::Track3D::size().

Referenced by pma::PMAlgTracker::extendTrack(), and pma::PMAlgTracker::reassignHits_1().

969 {
970  pma::Track3D* copy = new pma::Track3D(trk);
971  copy->AddHits(detProp, hits);
972 
973  mf::LogVerbatim("ProjectionMatchingAlg")
974  << "ext. track size: " << copy->size() << " #coll:" << copy->NHits(geo::kZ)
975  << " #ind2:" << copy->NHits(geo::kV) << " #ind1:" << copy->NHits(geo::kU);
976 
977  if (add_nodes) {
978  size_t nSegments = getSegCount_(copy->size());
979  int nNodes = nSegments - copy->Nodes().size() + 1; // n nodes to add
980  if (nNodes < 0) nNodes = 0;
981 
982  if (nNodes) {
983  mf::LogVerbatim("ProjectionMatchingAlg") << " add " << nNodes << " nodes";
984  copy->Optimize(detProp, nNodes, fOptimizationEps);
985  }
986  }
987  double g = copy->Optimize(detProp, 0, fFineTuningEps);
988  mf::LogVerbatim("ProjectionMatchingAlg") << " reopt done, g = " << g;
989 
990  return copy;
991 }
void AddHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
Definition: PmaTrack3D.cxx:404
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Planes which measure V.
Definition: geo_types.h:136
Planes which measure Z direction.
Definition: geo_types.h:138
static size_t getSegCount_(size_t trk_size)
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
Planes which measure U.
Definition: geo_types.h:135
unsigned int NHits(unsigned int view) const
Definition: PmaTrack3D.cxx:421
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
size_t size() const
Definition: PmaTrack3D.h:89
void pma::ProjectionMatchingAlg::FilterOutSmallParts ( const detinfo::DetectorPropertiesData detProp,
double  r2d,
const std::vector< art::Ptr< recob::Hit >> &  hits_in,
std::vector< art::Ptr< recob::Hit >> &  hits_out,
const TVector2 &  vtx2d 
) const

Get rid of small groups of hits around cascades; used to calculate cascade starting direction using the compact core cluster.

Definition at line 692 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), GetCloseHits_(), util::size(), and pma::WireDriftToCm().

Referenced by buildShowerSeg().

698 {
699  size_t min_size = hits_in.size() / 5;
700  if (min_size < 3) min_size = 3;
701 
702  std::vector<size_t> used;
703  std::vector<std::vector<art::Ptr<recob::Hit>>> sub_groups;
704  std::vector<art::Ptr<recob::Hit>> close_hits;
705 
706  float mindist2 = 99999.99;
707  size_t id_sub_groups_save = 0;
708  size_t id = 0;
709  while (GetCloseHits_(detProp, r2d, hits_in, used, close_hits)) {
710 
711  sub_groups.emplace_back(close_hits);
712 
713  for (size_t h = 0; h < close_hits.size(); ++h) {
714  TVector2 hi_cm = pma::WireDriftToCm(detProp,
715  close_hits[h]->WireID().Wire,
716  close_hits[h]->PeakTime(),
717  close_hits[h]->WireID().Plane,
718  close_hits[h]->WireID().TPC,
719  close_hits[h]->WireID().Cryostat);
720 
721  float dist2 = pma::Dist2(hi_cm, vtx2d);
722  if (dist2 < mindist2) {
723  id_sub_groups_save = id;
724  mindist2 = dist2;
725  }
726  }
727 
728  id++;
729  }
730 
731  for (size_t i = 0; i < sub_groups.size(); ++i) {
732  if (i == id_sub_groups_save) {
733  for (auto h : sub_groups[i])
734  hits_out.push_back(h);
735  }
736  }
737 
738  for (size_t i = 0; i < sub_groups.size(); ++i) {
739  if ((i != id_sub_groups_save) && (hits_out.size() < 10) && (sub_groups[i].size() > min_size))
740  for (auto h : sub_groups[i])
741  hits_out.push_back(h);
742  }
743 }
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:286
bool GetCloseHits_(const detinfo::DetectorPropertiesData &detProp, double r2d, const std::vector< art::Ptr< recob::Hit >> &hits_in, std::vector< size_t > &used, std::vector< art::Ptr< recob::Hit >> &hits_out) const
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
recob::tracking::Plane Plane
Definition: TrackState.h:17
bool pma::ProjectionMatchingAlg::GetCloseHits_ ( const detinfo::DetectorPropertiesData detProp,
double  r2d,
const std::vector< art::Ptr< recob::Hit >> &  hits_in,
std::vector< size_t > &  used,
std::vector< art::Ptr< recob::Hit >> &  hits_out 
) const
private

Definition at line 747 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), fGeom, detinfo::DetectorPropertiesData::GetXTicksCoefficient(), Has_(), recob::Hit::PeakTime(), geo::GeometryCore::Plane(), geo::WireID::Wire, recob::Hit::WireID(), and geo::PlaneGeo::WirePitch().

Referenced by FilterOutSmallParts().

752 {
753  hits_out.clear();
754 
755  const double gapMargin = 5.0; // can be changed to f(id_tpc1, id_tpc2)
756  size_t idx = 0;
757 
758  while ((idx < hits_in.size()) && Has_(used, idx))
759  idx++;
760 
761  if (idx < hits_in.size()) {
762  hits_out.push_back(hits_in[idx]);
763  used.push_back(idx);
764 
765  unsigned int tpc = hits_in[idx]->WireID().TPC;
766  unsigned int cryo = hits_in[idx]->WireID().Cryostat;
767  unsigned int view = hits_in[idx]->WireID().Plane;
768  double wirePitch = fGeom->Plane(geo::PlaneID(cryo, tpc, view)).WirePitch();
769  double driftPitch = detProp.GetXTicksCoefficient(tpc, cryo);
770 
771  double r2d2 = r2d * r2d;
772  double gapMargin2 = sqrt(2 * gapMargin * gapMargin);
773  gapMargin2 = (gapMargin2 + r2d) * (gapMargin2 + r2d);
774 
775  bool collect = true;
776  while (collect) {
777  collect = false;
778  for (size_t i = 0; i < hits_in.size(); i++)
779  if (!Has_(used, i)) {
780  art::Ptr<recob::Hit> const& hi = hits_in[i];
781  TVector2 hi_cm(wirePitch * hi->WireID().Wire, driftPitch * hi->PeakTime());
782 
783  bool accept = false;
784 
785  for (auto const& ho : hits_out) {
786 
787  TVector2 ho_cm(wirePitch * ho->WireID().Wire, driftPitch * ho->PeakTime());
788  double d2 = pma::Dist2(hi_cm, ho_cm);
789 
790  if (d2 < r2d2) {
791  accept = true;
792  break;
793  }
794  }
795  if (accept) {
796  collect = true;
797  hits_out.push_back(hi);
798  used.push_back(i);
799  }
800  }
801  }
802  return true;
803  }
804  else
805  return false;
806 }
bool Has_(const std::vector< size_t > &v, size_t idx) const
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
double GetXTicksCoefficient(int t, int c) const
geo::GeometryCore const * fGeom
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:378
size_t pma::ProjectionMatchingAlg::getSegCount_ ( size_t  trk_size)
staticprivate

Definition at line 439 of file ProjectionMatchingAlg.cxx.

Referenced by buildMultiTPCTrack(), buildTrack(), and extendTrack().

440 {
441  int const nSegments = 0.8 * trk_size / sqrt(trk_size);
442  return std::max(size_t{1}, static_cast<size_t>(nSegments));
443 }
void pma::ProjectionMatchingAlg::guideEndpoints ( const detinfo::DetectorPropertiesData clockData,
pma::Track3D trk,
const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &  hits 
) const

Add 3D reference points to clean endpoints of a track (both need to be in the same TPC).

Definition at line 1079 of file ProjectionMatchingAlg.cxx.

References addEndpointRef_(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), detinfo::DetectorPropertiesData::ConvertTicksToX(), fFineTuningEps, fGeom, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), pma::Segment3D::GetDirection3D(), geo::TPCGeo::HasPlane(), hits(), pma::Element3D::Length(), pma::Track3D::NextHit(), pma::Track3D::Optimize(), pma::Track3D::PrevHit(), pma::Track3D::Segments(), pma::Track3D::size(), and geo::GeometryCore::TPC().

Referenced by pma::PMAlgTrackingBase::guideEndpoints().

1083 {
1084  unsigned int tpc = trk.FrontTPC(), cryo = trk.FrontCryo();
1085  if ((tpc != trk.BackTPC()) || (cryo != trk.BackCryo())) {
1086  mf::LogWarning("ProjectionMatchingAlg") << "Please, apply before TPC stitching.";
1087  return;
1088  }
1089 
1090  const double maxCosXZ = 0.992546; // 7 deg
1091 
1092  pma::Segment3D* segFront = trk.Segments().front();
1093  if (trk.Segments().size() > 2) {
1094  pma::Segment3D* segFront1 = trk.Segments()[1];
1095  if ((segFront->Length() < 0.8) && (segFront1->Length() > 5.0)) segFront = segFront1;
1096  }
1097  pma::Vector3D dirFront = segFront->GetDirection3D();
1098  pma::Vector3D dirFrontXZ(dirFront.X(), 0., dirFront.Z());
1099  dirFrontXZ *= 1.0 / dirFrontXZ.R();
1100 
1101  pma::Segment3D* segBack = trk.Segments().back();
1102  if (trk.Segments().size() > 2) {
1103  pma::Segment3D* segBack1 = trk.Segments()[trk.Segments().size() - 2];
1104  if ((segBack->Length() < 0.8) && (segBack1->Length() > 5.0)) segBack = segBack1;
1105  }
1106  pma::Vector3D dirBack = segBack->GetDirection3D();
1107  pma::Vector3D dirBackXZ(dirBack.X(), 0., dirBack.Z());
1108  dirBackXZ *= 1.0 / dirBackXZ.R();
1109 
1110  if ((fabs(dirFrontXZ.Z()) < maxCosXZ) && (fabs(dirBackXZ.Z()) < maxCosXZ)) {
1111  return; // front & back are not parallel to wire planes => exit
1112  }
1113 
1114  unsigned int nPlanesFront = 0, nPlanesBack = 0;
1115  std::pair<int, int> wiresFront[3], wiresBack[3]; // wire index; index direction
1116  double xFront[3], xBack[3];
1117 
1118  for (unsigned int i = 0; i < 3; i++) {
1119  bool frontPresent = false, backPresent = false;
1120  if (fGeom->TPC(geo::TPCID(cryo, tpc)).HasPlane(i)) {
1121  int idxFront0 = trk.NextHit(-1, i);
1122  int idxBack0 = trk.PrevHit(trk.size(), i);
1123  if ((idxFront0 >= 0) && (idxFront0 < (int)trk.size()) && (idxBack0 >= 0) &&
1124  (idxBack0 < (int)trk.size())) {
1125  int idxFront1 = trk.NextHit(idxFront0, i);
1126  int idxBack1 = trk.PrevHit(idxBack0, i);
1127  if ((idxFront1 >= 0) && (idxFront1 < (int)trk.size()) && (idxBack1 >= 0) &&
1128  (idxBack1 < (int)trk.size())) {
1129  int wFront0 = trk[idxFront0]->Wire();
1130  int wBack0 = trk[idxBack0]->Wire();
1131 
1132  int wFront1 = trk[idxFront1]->Wire();
1133  int wBack1 = trk[idxBack1]->Wire();
1134 
1135  wiresFront[i].first = wFront0;
1136  wiresFront[i].second = wFront0 - wFront1;
1137  xFront[i] = detProp.ConvertTicksToX(trk[idxFront0]->PeakTime(), i, tpc, cryo);
1138 
1139  wiresBack[i].first = wBack0;
1140  wiresBack[i].second = wBack0 - wBack1;
1141  xBack[i] = detProp.ConvertTicksToX(trk[idxBack0]->PeakTime(), i, tpc, cryo);
1142 
1143  if (wiresFront[i].second) {
1144  if (wiresFront[i].second > 0)
1145  wiresFront[i].second = 1;
1146  else
1147  wiresFront[i].second = -1;
1148 
1149  frontPresent = true;
1150  nPlanesFront++;
1151  }
1152 
1153  if (wiresBack[i].second) {
1154  if (wiresBack[i].second > 0)
1155  wiresBack[i].second = 1;
1156  else
1157  wiresBack[i].second = -1;
1158 
1159  backPresent = true;
1160  nPlanesBack++;
1161  }
1162  }
1163  }
1164  }
1165  if (!frontPresent) { wiresFront[i].first = -1; }
1166  if (!backPresent) { wiresBack[i].first = -1; }
1167  }
1168 
1169  bool refAdded = false;
1170  if ((nPlanesFront > 1) && (fabs(dirFrontXZ.Z()) >= maxCosXZ)) {
1171  refAdded |= addEndpointRef_(detProp, trk, hits, wiresFront, xFront, tpc, cryo);
1172  }
1173 
1174  if ((nPlanesBack > 1) && (fabs(dirBackXZ.Z()) >= maxCosXZ)) {
1175  refAdded |= addEndpointRef_(detProp, trk, hits, wiresBack, xBack, tpc, cryo);
1176  }
1177  if (refAdded) {
1178  mf::LogVerbatim("ProjectionMatchingAlg") << "guide wire-plane-parallel track endpoints";
1179  double g = trk.Optimize(detProp, 0, 0.1 * fFineTuningEps);
1180  mf::LogVerbatim("ProjectionMatchingAlg") << " done, g = " << g;
1181  }
1182 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:147
bool addEndpointRef_(const detinfo::DetectorPropertiesData &detProp, pma::Track3D &trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &hits, std::pair< int, int > const *wires, double const *xPos, unsigned int tpc, unsigned int cryo) const
int PrevHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:911
geo::GeometryCore const * fGeom
unsigned int BackTPC() const
Definition: PmaTrack3D.h:121
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:34
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
unsigned int BackCryo() const
Definition: PmaTrack3D.h:122
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:269
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:118
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:119
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
pma::Vector3D GetDirection3D(void) const override
Get 3D direction cosines of this segment.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:894
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
size_t size() const
Definition: PmaTrack3D.h:89
double Length(void) const
Definition: PmaElement3D.h:53
void pma::ProjectionMatchingAlg::guideEndpoints ( const detinfo::DetectorPropertiesData clockData,
pma::Track3D trk,
pma::Track3D::ETrackEnd  endpoint,
const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &  hits 
) const

Add 3D reference points to clean endpoint of a track.

Definition at line 1185 of file ProjectionMatchingAlg.cxx.

References addEndpointRef_(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), detinfo::DetectorPropertiesData::ConvertTicksToX(), fFineTuningEps, fGeom, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), pma::Segment3D::GetDirection3D(), geo::TPCGeo::HasPlane(), hits(), pma::Track3D::kBegin, pma::Element3D::Length(), pma::Track3D::NextHit(), pma::Track3D::Optimize(), art::productstatus::present(), pma::Track3D::PrevHit(), pma::Track3D::Segments(), pma::Track3D::size(), and geo::GeometryCore::TPC().

1190 {
1191  const double maxCosXZ = 0.992546; // 7 deg
1192 
1193  unsigned int tpc, cryo;
1194  pma::Segment3D* seg0 = 0;
1195  pma::Segment3D* seg1 = 0;
1196 
1197  if (endpoint == pma::Track3D::kBegin) {
1198  tpc = trk.FrontTPC(), cryo = trk.FrontCryo();
1199  seg0 = trk.Segments().front();
1200  if (trk.Segments().size() > 2) { seg1 = trk.Segments()[1]; }
1201  }
1202  else {
1203  tpc = trk.BackTPC(), cryo = trk.BackCryo();
1204  seg0 = trk.Segments().back();
1205  if (trk.Segments().size() > 2) { seg1 = trk.Segments()[trk.Segments().size() - 2]; }
1206  }
1207  if (seg1 && (seg0->Length() < 0.8) && (seg1->Length() > 5.0)) { seg0 = seg1; }
1208  pma::Vector3D dir0 = seg0->GetDirection3D();
1209  pma::Vector3D dir0XZ(dir0.X(), 0., dir0.Z());
1210  dir0XZ *= 1.0 / dir0XZ.R();
1211 
1212  if (fabs(dir0XZ.Z()) < maxCosXZ) { return; } // not parallel to wire planes => exit
1213 
1214  unsigned int nPlanes = 0;
1215  std::pair<int, int> wires[3]; // wire index; index direction
1216  double x0[3];
1217 
1218  for (unsigned int i = 0; i < 3; i++) {
1219  bool present = false;
1220  if (fGeom->TPC(geo::TPCID(cryo, tpc)).HasPlane(i)) {
1221  int idx0 = -1, idx1 = -1;
1222  if (endpoint == pma::Track3D::kBegin) { idx0 = trk.NextHit(-1, i); }
1223  else {
1224  idx0 = trk.PrevHit(trk.size(), i);
1225  }
1226 
1227  if ((idx0 >= 0) && (idx0 < (int)trk.size()) && (trk[idx0]->TPC() == tpc) &&
1228  (trk[idx0]->Cryo() == cryo)) {
1229  if (endpoint == pma::Track3D::kBegin) { idx1 = trk.NextHit(idx0, i); }
1230  else {
1231  idx1 = trk.PrevHit(idx0, i);
1232  }
1233 
1234  if ((idx1 >= 0) && (idx1 < (int)trk.size()) && (trk[idx1]->TPC() == tpc) &&
1235  (trk[idx1]->Cryo() == cryo)) {
1236  int w0 = trk[idx0]->Wire();
1237  int w1 = trk[idx1]->Wire();
1238 
1239  wires[i].first = w0;
1240  wires[i].second = w0 - w1;
1241  x0[i] = detProp.ConvertTicksToX(trk[idx0]->PeakTime(), i, tpc, cryo);
1242 
1243  if (wires[i].second) {
1244  if (wires[i].second > 0)
1245  wires[i].second = 1;
1246  else
1247  wires[i].second = -1;
1248 
1249  present = true;
1250  nPlanes++;
1251  }
1252  }
1253  }
1254  }
1255  if (!present) { wires[i].first = -1; }
1256  }
1257 
1258  if ((nPlanes > 1) && (fabs(dir0XZ.Z()) >= maxCosXZ) &&
1259  addEndpointRef_(detProp, trk, hits, wires, x0, tpc, cryo)) {
1260  mf::LogVerbatim("ProjectionMatchingAlg") << "guide wire-plane-parallel track endpoint";
1261  double g = trk.Optimize(detProp, 0, 0.1 * fFineTuningEps);
1262  mf::LogVerbatim("ProjectionMatchingAlg") << " done, g = " << g;
1263  }
1264 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:147
bool addEndpointRef_(const detinfo::DetectorPropertiesData &detProp, pma::Track3D &trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit >>> &hits, std::pair< int, int > const *wires, double const *xPos, unsigned int tpc, unsigned int cryo) const
int PrevHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:911
geo::GeometryCore const * fGeom
unsigned int BackTPC() const
Definition: PmaTrack3D.h:121
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:34
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
unsigned int BackCryo() const
Definition: PmaTrack3D.h:122
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:269
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:118
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:119
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
pma::Vector3D GetDirection3D(void) const override
Get 3D direction cosines of this segment.
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:894
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
size_t size() const
Definition: PmaTrack3D.h:89
constexpr ProductStatus present() noexcept
Definition: ProductStatus.h:10
double Length(void) const
Definition: PmaElement3D.h:53
bool pma::ProjectionMatchingAlg::Has_ ( const std::vector< size_t > &  v,
size_t  idx 
) const
private

Definition at line 867 of file ProjectionMatchingAlg.cxx.

Referenced by GetCloseHits_().

868 {
869  for (auto c : v)
870  if (c == idx) return true;
871  return false;
872 }
bool pma::ProjectionMatchingAlg::isContained ( const pma::Track3D trk,
float  margin = 0.0F 
) const
inline

Test if hits at the track endpoinds do not stick out of TPC which they belong to. Here one can implement some configurable margin if needed for real data imeprfections.

Definition at line 167 of file ProjectionMatchingAlg.h.

References pma::Track3D::back(), pma::Track3D::FirstElement(), pma::Track3D::front(), hits(), pma::Track3D::LastElement(), pma::Hit3D::Point3D(), pma::Node3D::SameTPC(), and lar::dump::vector().

Referenced by pma::PMAlgTracker::matchCluster(), and pma::PMAlgTracker::reassignHits_1().

168  {
169  return (trk.FirstElement()->SameTPC(trk.front()->Point3D(), margin) &&
170  trk.LastElement()->SameTPC(trk.back()->Point3D(), margin));
171  }
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
pma::Node3D * FirstElement() const
Definition: PmaTrack3D.h:275
pma::Node3D * LastElement() const
Definition: PmaTrack3D.h:276
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:100
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:88
void pma::ProjectionMatchingAlg::mergeTracks ( const detinfo::DetectorPropertiesData detProp,
pma::Track3D dst,
pma::Track3D src,
bool  reopt 
) const

Add src to dst as it was its continuation; nodes of src are added to dst after its own nodes, hits of src are added to hits of dst, then dst is reoptimized.

Definition at line 1316 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddNode(), alignTracks(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), pma::Element3D::Cryo(), pma::Dist2(), fFineTuningEps, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), pma::Element3D::IsFrozen(), pma::Track3D::Length(), pma::Track3D::MakeProjection(), n, pma::Track3D::Nodes(), pma::Track3D::Optimize(), pma::Track3D::push_back(), pma::Track3D::ShiftEndsToHits(), pma::Track3D::size(), pma::Track3D::SortHits(), and pma::Element3D::TPC().

Referenced by buildMultiTPCTrack(), and pma::PMAlgTracker::mergeCoLinear().

1320 {
1321  if (!alignTracks(dst, src)) return;
1322 
1323  unsigned int tpc = src.FrontTPC();
1324  unsigned int cryo = src.FrontCryo();
1325  double lmean = dst.Length() / (dst.Nodes().size() - 1);
1326  if ((pma::Dist2(dst.Nodes().back()->Point3D(), src.Nodes().front()->Point3D()) > 0.5 * lmean) ||
1327  (tpc != dst.BackTPC()) || (cryo != dst.BackCryo())) {
1328  dst.AddNode(detProp, src.Nodes().front()->Point3D(), tpc, cryo);
1329  if (src.Nodes().front()->IsFrozen()) dst.Nodes().back()->SetFrozen(true);
1330  }
1331  for (size_t n = 1; n < src.Nodes().size(); n++) {
1332  pma::Node3D* node = src.Nodes()[n];
1333 
1334  dst.AddNode(detProp, src.Nodes()[n]->Point3D(), node->TPC(), node->Cryo());
1335 
1336  if (node->IsFrozen()) dst.Nodes().back()->SetFrozen(true);
1337  }
1338  for (size_t h = 0; h < src.size(); h++) {
1339  dst.push_back(detProp, src[h]->Hit2DPtr());
1340  }
1341  if (reopt) {
1342  double g = dst.Optimize(detProp, 0, fFineTuningEps);
1343  mf::LogVerbatim("ProjectionMatchingAlg") << " reopt after merging done, g = " << g;
1344  }
1345  else {
1346  dst.MakeProjection();
1347  }
1348 
1349  dst.SortHits();
1350  dst.ShiftEndsToHits();
1351 
1352  dst.MakeProjection();
1353  dst.SortHits();
1354 }
void MakeProjection()
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
unsigned int BackTPC() const
Definition: PmaTrack3D.h:121
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
unsigned int BackCryo() const
Definition: PmaTrack3D.h:122
void AddNode(pma::Node3D *node)
bool alignTracks(pma::Track3D &first, pma::Track3D &second) const
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:118
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:119
void push_back(pma::Hit3D *hit)
Definition: PmaTrack3D.h:77
int Cryo(void) const
Cryostat index or -1 if out of any cryostat.
Definition: PmaElement3D.h:37
double Length(size_t step=1) const
Definition: PmaTrack3D.h:94
bool IsFrozen(void) const
Check if the vertex 3D position is fixed.
Definition: PmaElement3D.h:105
Char_t n[5]
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
size_t size() const
Definition: PmaTrack3D.h:89
bool ShiftEndsToHits()
void pma::ProjectionMatchingAlg::RemoveNotEnabledHits ( pma::Track3D trk) const

Definition at line 876 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::release_at(), and pma::Track3D::size().

Referenced by ShortenSeg_().

877 {
878  size_t h = 0;
879  while (h < trk.size()) {
880  if (trk[h]->IsEnabled())
881  ++h;
882  else
883  (trk.release_at(h));
884  }
885 }
pma::Hit3D * release_at(size_t index)
Definition: PmaTrack3D.cxx:348
size_t size() const
Definition: PmaTrack3D.h:89
double pma::ProjectionMatchingAlg::selectInitialHits ( pma::Track3D trk,
unsigned int  view = geo::kZ,
unsigned int *  nused = 0 
) const

Intendet to calculate dQ/dx in the initial part of EM cascade; collection view is used by default, but it works also with other projections.

Definition at line 1357 of file ProjectionMatchingAlg.cxx.

References lar::util::absDiff(), pma::Hit3D::GetDistToProj(), pma::Hit3D::GetSegFraction(), pma::Track3D::HitDxByView(), pma::Track3D::NextHit(), pma::Track3D::size(), pma::Hit3D::SummedADC(), pma::Hit3D::TagOutlier(), pma::Hit3D::View2D(), and pma::Hit3D::Wire().

Referenced by ems::EMShower3D::ConvertFrom(), and ems::EMShower3D::ConvertFrom2().

1360 {
1361  for (size_t i = 0; i < trk.size(); i++) {
1362  pma::Hit3D* hit = trk[i];
1363  if (hit->View2D() == view) {
1364  if ((hit->GetDistToProj() > 0.5) || // more than 0.5cm away away from the segment
1365  (hit->GetSegFraction() < -1.0)) // projects before segment start (to check!!!)
1366  hit->TagOutlier(true);
1367  else
1368  hit->TagOutlier(false);
1369  }
1370  }
1371 
1372  unsigned int nhits = 0;
1373  double last_x, dx = 0.0, last_q, dq = 0.0, dqdx = 0.0;
1374  int ih = trk.NextHit(-1, view);
1375 
1376  pma::Hit3D* hit = trk[ih];
1377  pma::Hit3D* lastHit = hit;
1378 
1379  if ((ih >= 0) && (ih < (int)trk.size())) {
1380  hit->TagOutlier(true);
1381 
1382  ih = trk.NextHit(ih, view);
1383  while ((dx < 2.5) && (ih >= 0) && (ih < (int)trk.size())) {
1384  hit = trk[ih];
1385 
1386  if (lar::util::absDiff(hit->Wire(), lastHit->Wire()) > 2)
1387  break; // break on gap in wire direction
1388 
1389  last_x = trk.HitDxByView(ih, view);
1390  last_q = hit->SummedADC();
1391  if (dx + last_x < 3.0) {
1392  dq += last_q;
1393  dx += last_x;
1394  nhits++;
1395  }
1396  else
1397  break;
1398 
1399  lastHit = hit;
1400  ih = trk.NextHit(ih, view);
1401  }
1402  while ((ih >= 0) && (ih < (int)trk.size())) {
1403  hit = trk[ih];
1404  hit->TagOutlier(true);
1405  ih = trk.NextHit(ih, view);
1406  }
1407  }
1408  else {
1409  mf::LogError("ProjectionMatchingAlg") << "Initial part selection failed.";
1410  }
1411 
1412  if (!nhits) {
1413  mf::LogError("ProjectionMatchingAlg") << "Initial part too short to select useful hits.";
1414  }
1415 
1416  if (dx > 0.0) dqdx = dq / dx;
1417 
1418  if (nused) (*nused) = nhits;
1419 
1420  return dqdx;
1421 }
void TagOutlier(bool state) noexcept
Definition: PmaHit3D.h:90
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
unsigned int Wire() const noexcept
Definition: PmaHit3D.h:61
float SummedADC() const noexcept
Definition: PmaHit3D.h:64
float GetSegFraction() const noexcept
Definition: PmaHit3D.h:74
Detector simulation of raw signals on wires.
unsigned int View2D() const noexcept
Definition: PmaHit3D.h:60
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:69
double GetDistToProj() const
Definition: PmaHit3D.h:71
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:894
double HitDxByView(size_t index, unsigned int view) const
Definition: PmaTrack3D.cxx:993
size_t size() const
Definition: PmaTrack3D.h:89
void pma::ProjectionMatchingAlg::ShortenSeg_ ( const detinfo::DetectorPropertiesData detProp,
pma::Track3D trk,
const geo::TPCGeo tpcgeom 
) const
private

Definition at line 810 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), pma::Track3D::front(), pma::Track3D::GetMse(), pma::Track3D::Length(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), RemoveNotEnabledHits(), pma::Track3D::size(), pma::Track3D::SortHits(), and TestTrk_().

Referenced by buildShowerSeg().

813 {
814  double mse = trk.GetMse();
815  mf::LogWarning("ProjectionMatchingAlg") << "initial value of mse: " << mse;
816  while ((mse > 0.5) && TestTrk_(trk, tpcgeom)) {
817  mse = trk.GetMse();
818  for (size_t h = 0; h < trk.size(); ++h)
819  if (std::sqrt(pma::Dist2(trk[h]->Point3D(), trk.front()->Point3D())) > 0.8 * trk.Length())
820  trk[h]->SetEnabled(false);
821 
823 
824  // trk.Optimize(0.0001, false); // BUG: first argument missing; tentatively:
825  trk.Optimize(detProp, 0, 0.0001, false);
826  trk.SortHits();
827 
828  mf::LogWarning("ProjectionMatchingAlg") << " mse: " << mse;
829  if (mse == trk.GetMse()) break;
830 
831  mse = trk.GetMse();
832  }
833 
835 }
double GetMse(unsigned int view=geo::kUnknown) const
MSE of hits weighted with hit amplidudes and wire plane coefficients.
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
double Length(size_t step=1) const
Definition: PmaTrack3D.h:94
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
size_t size() const
Definition: PmaTrack3D.h:89
void RemoveNotEnabledHits(pma::Track3D &trk) const
bool TestTrk_(pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
unsigned int pma::ProjectionMatchingAlg::testHits ( detinfo::DetectorPropertiesData const &  detProp,
const pma::Track3D trk,
const std::vector< art::Ptr< recob::Hit >> &  hits,
double  eps = 1.0 
) const
inline

Count the number of hits that are closer than eps * fHitTestingDist2D to the track 2D projection.

Definition at line 157 of file ProjectionMatchingAlg.h.

References hits(), and pma::Track3D::TestHits().

Referenced by pma::PMAlgTracker::matchCluster(), pma::PMAlgTracker::matchTrack(), and pma::PMAlgTracker::reassignHits_1().

161  {
162  return trk.TestHits(detProp, hits, eps * fHitTestingDist2D);
163  }
unsigned int TestHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, double dist=0.4) const
Count close 2D hits.
Definition: PmaTrack3D.cxx:845
bool pma::ProjectionMatchingAlg::TestTrk_ ( pma::Track3D trk,
const geo::TPCGeo tpcgeom 
) const
private

Definition at line 839 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), pma::Track3D::front(), geo::TPCGeo::HasPlane(), pma::Track3D::NEnabledHits(), pma::Hit3D::Point3D(), and pma::Track3D::size().

Referenced by ShortenSeg_().

840 {
841  bool test = false;
842 
843  if (tpcgeom.HasPlane(0) && tpcgeom.HasPlane(1)) {
844  if ((trk.NEnabledHits(0) > 5) && (trk.NEnabledHits(1) > 5)) test = true;
845  }
846  else if (tpcgeom.HasPlane(0) && tpcgeom.HasPlane(2)) {
847  if ((trk.NEnabledHits(0) > 5) && (trk.NEnabledHits(2) > 5)) test = true;
848  }
849  else if (tpcgeom.HasPlane(1) && tpcgeom.HasPlane(2)) {
850  if ((trk.NEnabledHits(1) > 5) && (trk.NEnabledHits(2) > 5)) test = true;
851  }
852 
853  double length = 0.0;
854  for (size_t h = 0; h < trk.size(); ++h) {
855  if (!trk[h]->IsEnabled()) break;
856  length = std::sqrt(pma::Dist2(trk.front()->Point3D(), trk[h]->Point3D()));
857  }
858 
859  mf::LogWarning("ProjectionMatchingAlg") << "length of segment: " << length;
860  if (length < 3.0) test = false; // cm
861 
862  return test;
863 }
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:147
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
unsigned int NEnabledHits(unsigned int view=geo::kUnknown) const
Definition: PmaTrack3D.cxx:430
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
size_t size() const
Definition: PmaTrack3D.h:89
std::vector< pma::Hit3D * > pma::ProjectionMatchingAlg::trimTrackToVolume ( pma::Track3D trk,
TVector3  p0,
TVector3  p1 
) const

Definition at line 1267 of file ProjectionMatchingAlg.cxx.

1270 {
1271  return {};
1272 }
double pma::ProjectionMatchingAlg::twoViewFraction ( pma::Track3D trk) const

Calculate the fraction of trajectory seen by two 2D projections at least; even a prfect track starts/stops with the hit from one 2D view, then hits from other views come, which results with the fraction value high, but always < 1.0; wrong cluster matchings or incomplete tracks give significantly lower values.

Definition at line 418 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::DisableSingleViewEnds(), pma::Track3D::Length(), pma::Track3D::SelectHits(), and pma::Track3D::size().

Referenced by buildTrack(), and pma::PMAlgTracker::matchCluster().

419 {
420  trk.SelectHits();
421  trk.DisableSingleViewEnds();
422 
423  size_t idx = 0;
424  while ((idx < trk.size() - 1) && !trk[idx]->IsEnabled())
425  idx++;
426  double l0 = trk.Length(0, idx + 1);
427 
428  idx = trk.size() - 1;
429  while ((idx > 1) && !trk[idx]->IsEnabled())
430  idx--;
431  double l1 = trk.Length(idx - 1, trk.size() - 1);
432 
433  trk.SelectHits();
434 
435  return 1.0 - (l0 + l1) / trk.Length();
436 }
bool SelectHits(float fraction=1.0F)
unsigned int DisableSingleViewEnds()
double Length(size_t step=1) const
Definition: PmaTrack3D.h:94
size_t size() const
Definition: PmaTrack3D.h:89
double pma::ProjectionMatchingAlg::validate ( const detinfo::DetectorPropertiesData detProp,
const lariov::ChannelStatusProvider &  channelStatus,
const pma::Track3D trk,
const std::vector< art::Ptr< recob::Hit >> &  hits 
) const

Calculate the fraction of the track that is closer than fTrkValidationDist2D to any hit from hits in their plane (a plane that was not used to build the track). Hits should be preselected, so all belong to the same plane.

Definition at line 252 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), pma::Track3D::Cryos(), pma::Dist2(), pma::Track3D::Dist2(), f, fGeom, pma::Track3D::front(), fTrkValidationDist2D, pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), hits(), geo::TPCGeo::Plane(), geo::GeometryCore::Plane(), geo::GeometryCore::PlaneWireToChannel(), pma::Hit3D::Point3D(), pma::Node3D::SameTPC(), pma::Track3D::Segments(), lar::to_element, geo::vect::toPoint(), pma::Element3D::TPC(), geo::GeometryCore::TPC(), pma::Track3D::TPCs(), lar::dump::vector(), geo::GeometryCore::WireCoordinate(), pma::Track3D::WireDriftRange(), pma::WireDriftToCm(), and geo::PlaneGeo::WirePitch().

Referenced by pma::PMAlgTracker::validate().

256 {
257  if (hits.empty()) { return 0; }
258 
259  double max_d = fTrkValidationDist2D;
260  double const max_d2 = max_d * max_d;
261  unsigned int nAll = 0, nPassed = 0;
262  unsigned int testPlane = hits.front()->WireID().Plane;
263 
264  std::vector<unsigned int> trkTPCs = trk.TPCs();
265  std::vector<unsigned int> trkCryos = trk.Cryos();
266  std::map<std::pair<unsigned int, unsigned int>, std::pair<TVector2, TVector2>> ranges;
267  std::map<std::pair<unsigned int, unsigned int>, double> wirePitch;
268  for (auto const& [t, c] : views::cartesian_product(trkTPCs, trkCryos)) {
269  ranges[{t, c}] = trk.WireDriftRange(detProp, testPlane, t, c);
270  wirePitch[{t, c}] = fGeom->TPC(geo::TPCID(c, t)).Plane(testPlane).WirePitch();
271  }
272 
273  unsigned int tpc, cryo;
274  std::map<std::pair<unsigned int, unsigned int>, std::vector<pma::Vector2D>> all_close_points;
275 
276  for (auto const& h :
277  hits | views::transform(to_element) | views::filter(hits_on_plane{testPlane})) {
278  tpc = h.WireID().TPC;
279  cryo = h.WireID().Cryostat;
280  std::pair<unsigned int, unsigned int> tpc_cryo(tpc, cryo);
281  std::pair<TVector2, TVector2> rect = ranges[tpc_cryo];
282 
283  if ((h.WireID().Wire > rect.first.X() - 10) && // chceck only hits in the rectangle around
284  (h.WireID().Wire < rect.second.X() + 10) && // the track projection, it is faster than
285  (h.PeakTime() > rect.first.Y() - 100) && // calculation of trk.Dist2(p2d, testPlane)
286  (h.PeakTime() < rect.second.Y() + 100)) {
287  TVector2 p2d(wirePitch[tpc_cryo] * h.WireID().Wire,
288  detProp.ConvertTicksToX(h.PeakTime(), testPlane, tpc, cryo));
289 
290  double const d2 = trk.Dist2(p2d, testPlane, tpc, cryo);
291 
292  if (d2 < max_d2) all_close_points[tpc_cryo].emplace_back(p2d.X(), p2d.Y());
293  }
294  }
295 
296  // then check how points-close-to-the-track-projection are distributed along
297  // the track, namely: are there track sections crossing empty spaces, except
298  // dead wires?
299  pma::Vector3D p(
300  trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
301  for (auto const* seg : trk.Segments()) {
302  if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
303  {
304  p = seg->End();
305  continue;
306  }
307  pma::Vector3D p0 = seg->Start();
308  pma::Vector3D p1 = seg->End();
309 
310  pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());
311 
312  tpc = seg->TPC();
313  cryo = seg->Cryo();
314 
315  pma::Vector3D dc = step * seg->GetDirection3D();
316 
317  auto const& points = all_close_points[{tpc, cryo}];
318 
319  double f = pma::GetSegmentProjVector(p, p0, p1);
320 
321  geo::PlaneID const planeID{cryo, tpc, testPlane};
322  double const wirepitch = fGeom->Plane(planeID).WirePitch();
323  while ((f < 1.0) && node->SameTPC(p)) {
324  pma::Vector2D p2d(fGeom->WireCoordinate(toPoint(p), planeID), p.X());
325  geo::WireID const wireID{planeID, static_cast<unsigned int>(p2d.X())};
326  if (fGeom->HasWire(wireID)) {
328  if (channelStatus.IsGood(ch)) {
329  if (points.size()) {
330  p2d.SetX(wirepitch * p2d.X());
331  for (const auto& h : points) {
332  if (pma::Dist2(p2d, h) < max_d2) {
333  nPassed++;
334  break;
335  }
336  }
337  }
338  nAll++;
339  }
340  //else mf::LogVerbatim("ProjectionMatchingAlg")
341  // << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
342  }
343 
344  p += dc;
345  f = pma::GetSegmentProjVector(p, p0, p1);
346  }
347  p = seg->End();
348  }
349 
350  if (nAll > 0) {
351  double v = nPassed / (double)nAll;
352  mf::LogVerbatim("ProjectionMatchingAlg") << " trk fraction ok: " << v;
353  return v;
354  }
355 
356  return 1.0;
357 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
constexpr to_element_t to_element
Definition: ToElement.h:25
double Dist2(const TVector2 &p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
bool HasWire(WireID const &wireid) const
Returns whether we have the specified wire.
geo::GeometryCore const * fGeom
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:448
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:34
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
std::vector< unsigned int > Cryos() const
Definition: PmaTrack3D.cxx:466
TFile f
Definition: plotHisto.C:6
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:269
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:100
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
double ConvertTicksToX(double ticks, int p, int t, int c) const
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:112
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
Definition: Utilities.h:33
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:252
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:378
std::pair< TVector2, TVector2 > WireDriftRange(detinfo::DetectorPropertiesData const &detProp, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:484
double pma::ProjectionMatchingAlg::validate ( const detinfo::DetectorPropertiesData detProp,
const lariov::ChannelStatusProvider &  channelStatus,
const TVector3 &  p0,
const TVector3 &  p1,
const std::vector< art::Ptr< recob::Hit >> &  hits,
unsigned int  testView,
unsigned int  tpc,
unsigned int  cryo 
) const

Calculate the fraction of the 3D segment that is closer than fTrkValidationDist2D to any hit from hits in the testPlane of TPC/Cryo. Hits from the testPlane are preselected by this function among all provided (so a bit slower than fn above).

double pma::ProjectionMatchingAlg::validate_on_adc ( const detinfo::DetectorPropertiesData detProp,
const lariov::ChannelStatusProvider &  channelStatus,
const pma::Track3D trk,
const img::DataProviderAlg adcImage,
float  thr 
) const

Calculate the fraction of the track that is close to non-empty pixel (above thr value) in the ADC image of the testView (a view that was not used to build the track).

Definition at line 47 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertXToTicks(), pma::Track3D::Cryos(), f, fGeom, pma::Track3D::front(), pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), img::DataProviderAlg::Plane(), geo::PlaneID::Plane, geo::GeometryCore::PlaneWireToChannel(), pma::Hit3D::Point3D(), img::DataProviderAlg::poolMax(), pma::Node3D::SameTPC(), pma::Track3D::Segments(), geo::vect::toPoint(), pma::Element3D::TPC(), pma::Track3D::TPCs(), geo::GeometryCore::WireCoordinate(), and recob::Hit::WireID().

Referenced by pma::PMAlgTracker::validate().

53 {
54  unsigned int nAll = 0, nPassed = 0;
55  unsigned int testPlane = adcImage.Plane();
56 
57  std::vector<unsigned int> trkTPCs = trk.TPCs();
58  std::vector<unsigned int> trkCryos = trk.Cryos();
59 
60  // check how pixels with a high signal are distributed along the track
61  // namely: are there track sections crossing empty spaces, except dead wires?
62  pma::Vector3D p(
63  trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
64  for (auto const* seg : trk.Segments()) {
65  if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
66  {
67  p = seg->End();
68  continue;
69  }
70  pma::Vector3D p0 = seg->Start();
71  pma::Vector3D p1 = seg->End();
72 
73  pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());
74 
75  unsigned tpc = seg->TPC();
76  unsigned cryo = seg->Cryo();
77 
78  pma::Vector3D dc = step * seg->GetDirection3D();
79 
80  double f = pma::GetSegmentProjVector(p, p0, p1);
81  while ((f < 1.0) && node->SameTPC(p)) {
82  pma::Vector2D p2d(fGeom->WireCoordinate(toPoint(p), geo::PlaneID{cryo, tpc, testPlane}),
83  p.X());
84  geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
85 
86  int widx = (int)p2d.X();
87  int didx = (int)detProp.ConvertXToTicks(p2d.Y(), testPlane, tpc, cryo);
88 
89  if (fGeom->HasWire(wireID)) {
91  if (channelStatus.IsGood(ch)) {
92  float max_adc = adcImage.poolMax(widx, didx, 2); // +/- 2 wires, can be parameterized
93  if (max_adc > thr) nPassed++;
94 
95  nAll++;
96  }
97  }
98 
99  p += dc;
100  f = pma::GetSegmentProjVector(p, p0, p1);
101  }
102 
103  p = seg->End(); // need to have it at the end due to the p in the first iter
104  // set to the hit position, not segment start
105  }
106 
107  if (nAll > 0) {
108  double v = nPassed / (double)nAll;
109  mf::LogVerbatim("ProjectionMatchingAlg") << " trk fraction ok: " << v;
110  return v;
111  }
112 
113  return 1.0;
114 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
bool HasWire(WireID const &wireid) const
Returns whether we have the specified wire.
geo::GeometryCore const * fGeom
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:448
unsigned int Plane() const
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:34
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
std::vector< unsigned int > Cryos() const
Definition: PmaTrack3D.cxx:466
TFile f
Definition: plotHisto.C:6
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:269
double ConvertXToTicks(double X, int p, int t, int c) const
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:100
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:112
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
Definition: Utilities.h:33
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
float poolMax(int wire, int drift, size_t r=0) const
Pool max value in a patch around the wire/drift pixel.
double pma::ProjectionMatchingAlg::validate_on_adc_test ( const detinfo::DetectorPropertiesData detProp,
const lariov::ChannelStatusProvider &  channelStatus,
const pma::Track3D trk,
const img::DataProviderAlg adcImage,
const std::vector< art::Ptr< recob::Hit >> &  hits,
TH1F *  histoPassing,
TH1F *  histoRejected 
) const

Calculate the fraction of the track that is closer than fTrkValidationDist2D to any hit from hits in the testView (a view that was not used to build the track). Creates also histograms of values in pixels for the passing and rejected points on the track, so the threshold value for the ADC-based calibration can be estimated.

Definition at line 125 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), detinfo::DetectorPropertiesData::ConvertXToTicks(), pma::Track3D::Cryos(), pma::Dist2(), pma::Track3D::Dist2(), f, fGeom, pma::Track3D::front(), fTrkValidationDist2D, pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), hits(), img::DataProviderAlg::Plane(), geo::TPCGeo::Plane(), geo::GeometryCore::Plane(), geo::GeometryCore::PlaneWireToChannel(), pma::Hit3D::Point3D(), img::DataProviderAlg::poolMax(), pma::Node3D::SameTPC(), pma::Track3D::Segments(), lar::to_element, geo::vect::toPoint(), pma::Element3D::TPC(), geo::GeometryCore::TPC(), pma::Track3D::TPCs(), geo::GeometryCore::WireCoordinate(), pma::Track3D::WireDriftRange(), and geo::PlaneGeo::WirePitch().

Referenced by pma::PMAlgTracker::validate().

133 {
134  double max_d = fTrkValidationDist2D;
135  double const max_d2 = max_d * max_d;
136  unsigned int nAll = 0, nPassed = 0;
137  unsigned int testPlane = adcImage.Plane();
138 
139  std::vector<unsigned int> trkTPCs = trk.TPCs();
140  std::vector<unsigned int> trkCryos = trk.Cryos();
141  std::map<std::pair<unsigned int, unsigned int>, std::pair<TVector2, TVector2>> ranges;
142  std::map<std::pair<unsigned int, unsigned int>, double> wirePitch;
143  for (auto const& [t, c] : views::cartesian_product(trkTPCs, trkCryos)) {
144  ranges[{t, c}] = trk.WireDriftRange(detProp, testPlane, t, c);
145  wirePitch[{t, c}] = fGeom->TPC(geo::TPCID(c, t)).Plane(testPlane).WirePitch();
146  }
147 
148  unsigned int tpc, cryo;
149  std::map<std::pair<unsigned int, unsigned int>, std::vector<pma::Vector2D>> all_close_points;
150 
151  for (auto const& h :
152  hits | views::transform(to_element) | views::filter(hits_on_plane{testPlane})) {
153  tpc = h.WireID().TPC;
154  cryo = h.WireID().Cryostat;
155  std::pair<unsigned int, unsigned int> tpc_cryo(tpc, cryo);
156  std::pair<TVector2, TVector2> rect = ranges[tpc_cryo];
157 
158  if ((h.WireID().Wire > rect.first.X() - 10) && // check only hits in the rectangle around
159  (h.WireID().Wire < rect.second.X() + 10) && // the track projection, it is faster than
160  (h.PeakTime() > rect.first.Y() - 100) && // calculation of trk.Dist2(p2d, testPlane)
161  (h.PeakTime() < rect.second.Y() + 100)) {
162  TVector2 p2d(wirePitch[tpc_cryo] * h.WireID().Wire,
163  detProp.ConvertTicksToX(h.PeakTime(), testPlane, tpc, cryo));
164 
165  double const d2 = trk.Dist2(p2d, testPlane, tpc, cryo);
166 
167  if (d2 < max_d2) { all_close_points[tpc_cryo].emplace_back(p2d.X(), p2d.Y()); }
168  }
169  }
170 
171  // then check how points-close-to-the-track-projection are distributed along
172  // the track, namely: are there track sections crossing empty spaces, except
173  // dead wires?
174  pma::Vector3D p(
175  trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
176  for (auto const* seg : trk.Segments()) {
177  if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
178  {
179  p = seg->End();
180  continue;
181  }
182  pma::Vector3D p0 = seg->Start();
183  pma::Vector3D p1 = seg->End();
184 
185  pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());
186 
187  tpc = seg->TPC();
188  cryo = seg->Cryo();
189 
190  pma::Vector3D dc = step * seg->GetDirection3D();
191 
192  auto const& points = all_close_points[{tpc, cryo}];
193 
194  double f = pma::GetSegmentProjVector(p, p0, p1);
195 
196  double wirepitch = fGeom->Plane(geo::PlaneID(cryo, tpc, testPlane)).WirePitch();
197  while ((f < 1.0) && node->SameTPC(p)) {
198  geo::PlaneID const planeID{cryo, tpc, testPlane};
199  pma::Vector2D p2d(fGeom->WireCoordinate(toPoint(p), planeID), p.X());
200  geo::WireID const wireID{planeID, static_cast<unsigned int>(p2d.X())};
201 
202  int widx = (int)p2d.X();
203  int didx = (int)detProp.ConvertXToTicks(p2d.Y(), planeID);
204 
205  if (fGeom->HasWire(wireID)) {
207  if (channelStatus.IsGood(ch)) {
208  bool is_close = false;
209  float max_adc = adcImage.poolMax(widx, didx, 2);
210 
211  if (points.size()) {
212  p2d.SetX(wirepitch * p2d.X());
213  for (const auto& h : points) {
214  if (pma::Dist2(p2d, h) < max_d2) {
215  is_close = true;
216  nPassed++;
217  break;
218  }
219  }
220  }
221  nAll++;
222 
223  // now fill the calibration histograms
224  if (is_close) {
225  if (histoPassing) histoPassing->Fill(max_adc);
226  }
227  else {
228  if (histoRejected) histoRejected->Fill(max_adc);
229  }
230  }
231  //else mf::LogVerbatim("ProjectionMatchingAlg")
232  // << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
233  }
234 
235  p += dc;
236  f = pma::GetSegmentProjVector(p, p0, p1);
237  }
238  p = seg->End();
239  }
240 
241  if (nAll > 0) {
242  double v = nPassed / (double)nAll;
243  mf::LogVerbatim("ProjectionMatchingAlg") << " trk fraction ok: " << v;
244  return v;
245  }
246 
247  return 1.0;
248 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
constexpr to_element_t to_element
Definition: ToElement.h:25
double Dist2(const TVector2 &p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
bool HasWire(WireID const &wireid) const
Returns whether we have the specified wire.
geo::GeometryCore const * fGeom
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:448
unsigned int Plane() const
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:34
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
std::vector< unsigned int > Cryos() const
Definition: PmaTrack3D.cxx:466
TFile f
Definition: plotHisto.C:6
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:269
double ConvertXToTicks(double X, int p, int t, int c) const
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:100
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
double ConvertTicksToX(double ticks, int p, int t, int c) const
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:112
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
Definition: Utilities.h:33
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:252
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
float poolMax(int wire, int drift, size_t r=0) const
Pool max value in a patch around the wire/drift pixel.
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:378
std::pair< TVector2, TVector2 > WireDriftRange(detinfo::DetectorPropertiesData const &detProp, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:484

Member Data Documentation

double const pma::ProjectionMatchingAlg::fFineTuningEps
private
double const pma::ProjectionMatchingAlg::fHitTestingDist2D
private

Definition at line 328 of file ProjectionMatchingAlg.h.

Referenced by ProjectionMatchingAlg().

double const pma::ProjectionMatchingAlg::fMinTwoViewFraction
private

Definition at line 331 of file ProjectionMatchingAlg.h.

Referenced by buildTrack(), and ProjectionMatchingAlg().

double const pma::ProjectionMatchingAlg::fOptimizationEps
private

Definition at line 319 of file ProjectionMatchingAlg.h.

Referenced by buildMultiTPCTrack(), buildTrack(), and extendTrack().

double const pma::ProjectionMatchingAlg::fTrkValidationDist2D
private

Definition at line 326 of file ProjectionMatchingAlg.h.

Referenced by ProjectionMatchingAlg(), validate(), and validate_on_adc_test().


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