All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
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_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::Track3D * buildTrack (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
 
pma::Track3D * buildMultiTPCTrack (const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits) const
 
pma::Track3D * buildShowerSeg (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, const pma::Vector3D &vtx) 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
 
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 TVector3 &point) const
 
pma::Track3D * buildSegment (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::Track3D * 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. 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 30 of file ProjectionMatchingAlg.cxx.

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

Definition at line 102 of file ProjectionMatchingAlg.h.

103  : ProjectionMatchingAlg(fhicl::Table<Config>(pset, {})())
pma::ProjectionMatchingAlg::ProjectionMatchingAlg
ProjectionMatchingAlg(const Config &config)
Definition: ProjectionMatchingAlg.cxx:30

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 1042 of file ProjectionMatchingAlg.cxx.

1050 {
1051  double x = 0.0, y = 0.0, z = 0.0;
1052  std::vector<std::pair<int, unsigned int>> wire_view;
1053  for (unsigned int i = 0; i < 3; i++)
1054  if (wires[i].first >= 0) {
1055  const auto hiter = hits.find(i);
1056  if (hiter != hits.end()) {
1057  if (chkEndpointHits_(detProp,
1058  wires[i].first,
1059  wires[i].second,
1060  xPos[i],
1061  i,
1062  tpc,
1063  cryo,
1064  trk,
1065  hiter->second)) {
1066  x += xPos[i];
1067  wire_view.emplace_back(wires[i].first, i);
1068  }
1069  }
1070  }
1071  if (wire_view.size() > 1) {
1072  x /= wire_view.size();
1073  fGeom->IntersectionPoint(wire_view[0].first,
1074  wire_view[1].first,
1075  wire_view[0].second,
1076  wire_view[1].second,
1077  cryo,
1078  tpc,
1079  y,
1080  z);
1081  trk.AddRefPoint(x, y, z);
1082  mf::LogVerbatim("ProjectionMatchingAlg")
1083  << "trk tpc:" << tpc << " size:" << trk.size() << " add ref.point (" << x << "; " << y << "; "
1084  << z << ")";
1085  return true;
1086  }
1087  else {
1088  mf::LogVerbatim("ProjectionMatchingAlg")
1089  << "trk tpc:" << tpc << " size:" << trk.size()
1090  << " wire-plane-parallel track, but need two clean views of endpoint";
1091  return false;
1092  }
1093 }
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
second
second
Definition: detectorproperties_icarus.fcl:8
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
pma::Track3D::AddRefPoint
void AddRefPoint(const TVector3 &p)
Definition: PmaTrack3D.h:265
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
geo::GeometryCore::IntersectionPoint
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Definition: GeometryCore.cxx:1746
pma::ProjectionMatchingAlg::chkEndpointHits_
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
Definition: ProjectionMatchingAlg.cxx:1005
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
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 1293 of file ProjectionMatchingAlg.cxx.

1294 {
1295  unsigned int k = 0;
1296  double const distFF = pma::Dist2(first.front()->Point3D(), second.front()->Point3D());
1297  double dist = distFF;
1298 
1299  double const distFB = pma::Dist2(first.front()->Point3D(), second.back()->Point3D());
1300  if (distFB < dist) {
1301  k = 1;
1302  dist = distFB;
1303  }
1304 
1305  double const distBF = pma::Dist2(first.back()->Point3D(), second.front()->Point3D());
1306  if (distBF < dist) {
1307  k = 2;
1308  dist = distBF;
1309  }
1310 
1311  double distBB = pma::Dist2(first.back()->Point3D(), second.back()->Point3D());
1312  if (distBB < dist) {
1313  k = 3;
1314  dist = distBB;
1315  }
1316 
1317  switch (k) // flip to get dst end before src start, do not merge if track's order reversed
1318  {
1319  case 0:
1320  first.Flip(); // detProp);
1321  break;
1322  case 1: mf::LogError("PMAlgTrackMaker") << "Tracks in reversed order."; return false;
1323  case 2: break;
1324  case 3:
1325  second.Flip(); // detProp);
1326  break;
1327  default:
1328  throw cet::exception("pma::ProjectionMatchingAlg")
1329  << "alignTracks: should never happen." << std::endl;
1330  }
1331  return true;
1332 }
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
pma::Track3D::back
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:106
k
pdgs k
Definition: selectors.fcl:22
pma::Track3D::Flip
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:534
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 267 of file ProjectionMatchingAlg.h.

271  {
272  trk.AutoFlip(dir, thr, n);
273  };
pma::Track3D::AutoFlip
void AutoFlip(pma::Track3D::EDirection dir, double thr=0.0, unsigned int n=0)
Definition: PmaTrack3D.cxx:673
dropbox.dir
tuple dir
Definition: dropbox.py:28
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 505 of file ProjectionMatchingAlg.cxx.

507 {
508  std::map<unsigned int, std::vector<art::Ptr<recob::Hit>>> hits_by_tpc;
509  for (auto const& h : hits) {
510  hits_by_tpc[h->WireID().TPC].push_back(h);
511  }
512 
513  std::vector<pma::Track3D*> tracks;
514  for (auto const& hsel : hits_by_tpc) {
515  pma::Track3D* trk = buildSegment(detProp, hsel.second);
516  if (trk) tracks.push_back(trk);
517  }
518 
519  bool need_reopt = false;
520  while (tracks.size() > 1) {
521  need_reopt = true;
522 
523  pma::Track3D* first = tracks.front();
524  pma::Track3D* best = 0;
525  double d, dmin = 1.0e12;
526  size_t t_best = 0, cfg = 0;
527  for (size_t t = 1; t < tracks.size(); ++t) {
528  pma::Track3D* second = tracks[t];
529 
530  d = pma::Dist2(first->front()->Point3D(), second->front()->Point3D());
531  if (d < dmin) {
532  dmin = d;
533  best = second;
534  t_best = t;
535  cfg = 0;
536  }
537 
538  d = pma::Dist2(first->front()->Point3D(), second->back()->Point3D());
539  if (d < dmin) {
540  dmin = d;
541  best = second;
542  t_best = t;
543  cfg = 1;
544  }
545 
546  d = pma::Dist2(first->back()->Point3D(), second->front()->Point3D());
547  if (d < dmin) {
548  dmin = d;
549  best = second;
550  t_best = t;
551  cfg = 2;
552  }
553 
554  d = pma::Dist2(first->back()->Point3D(), second->back()->Point3D());
555  if (d < dmin) {
556  dmin = d;
557  best = second;
558  t_best = t;
559  cfg = 3;
560  }
561  }
562  if (best) {
563  switch (cfg) {
564  default:
565  case 0:
566  case 1:
567  mergeTracks(detProp, *best, *first, false);
568  tracks[0] = best;
569  delete first;
570  break;
571 
572  case 2:
573  case 3:
574  mergeTracks(detProp, *first, *best, false);
575  delete best;
576  break;
577  }
578  tracks.erase(tracks.begin() + t_best);
579  }
580  else
581  break; // should not happen
582  }
583 
584  pma::Track3D* trk = 0;
585  if (!tracks.empty()) {
586  trk = tracks.front();
587  if (need_reopt) {
588  double g = trk->Optimize(detProp, 0, fOptimizationEps);
589  mf::LogVerbatim("ProjectionMatchingAlg")
590  << " reopt after merging initial tpc segments: done, g = " << g;
591  }
592 
593  int nSegments = getSegCount_(trk->size()) - trk->Segments().size() + 1;
594  if (nSegments > 0) // need to add segments
595  {
596  double g = 0.0;
597  size_t nNodes = (size_t)(nSegments - 1); // n nodes to add
598  if (nNodes) {
599  mf::LogVerbatim("ProjectionMatchingAlg") << " optimize trk (add " << nSegments << " seg)";
600 
601  g = trk->Optimize(detProp, nNodes, fOptimizationEps, false, true, 25,
602  10); // build nodes
603  mf::LogVerbatim("ProjectionMatchingAlg") << " nodes done, g = " << g;
604  trk->SelectAllHits();
605  }
606  g = trk->Optimize(detProp, 0, fFineTuningEps); // final tuning
607  mf::LogVerbatim("ProjectionMatchingAlg") << " tune done, g = " << g;
608  }
609  trk->SortHits();
610  }
611  return trk;
612 }
second
second
Definition: detectorproperties_icarus.fcl:8
tracks
ClusterModuleLabel join with tracks
Definition: trackfinderalgorithms.fcl:56
pma::Track3D::SelectAllHits
bool SelectAllHits()
Definition: PmaTrack3D.cxx:2942
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::ProjectionMatchingAlg::getSegCount_
static size_t getSegCount_(size_t trk_size)
Definition: ProjectionMatchingAlg.cxx:444
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
h
while getopts h
Definition: getStatistics.sh:13
pma::Track3D::Segments
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:329
pma::ProjectionMatchingAlg::buildSegment
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
Definition: ProjectionMatchingAlg.cxx:900
pma::ProjectionMatchingAlg::mergeTracks
void mergeTracks(const detinfo::DetectorPropertiesData &detProp, pma::Track3D &dst, pma::Track3D &src, bool reopt) const
Definition: ProjectionMatchingAlg.cxx:1335
pma::Track3D::back
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:106
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
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 900 of file ProjectionMatchingAlg.cxx.

903 {
904  pma::Track3D* trk = new pma::Track3D();
905  trk->SetEndSegWeight(0.001F);
906  trk->AddHits(detProp, hits_1);
907  trk->AddHits(detProp, hits_2);
908 
909  if (trk->HasTwoViews() && (trk->TPCs().size() == 1)) // now only in single tpc
910  {
911  if (!trk->Initialize(detProp, 0.001F)) {
912  mf::LogWarning("ProjectionMatchingAlg") << "initialization failed, delete segment";
913  delete trk;
914  return 0;
915  }
916  trk->Optimize(detProp, 0, fFineTuningEps);
917 
918  trk->SortHits();
919  return trk;
920  }
921  else {
922  mf::LogWarning("ProjectionMatchingAlg") << "need at least two views in single tpc";
923  delete trk;
924  return 0;
925  }
926 }
pma::Track3D::AddHits
void AddHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
Definition: PmaTrack3D.cxx:404
pma::Track3D::Initialize
bool Initialize(detinfo::DetectorPropertiesData const &detProp, float initEndSegW=0.05F)
Definition: PmaTrack3D.cxx:78
pma::Track3D::TPCs
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:453
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
pma::Track3D::SetEndSegWeight
void SetEndSegWeight(float value) noexcept
Definition: PmaTrack3D.h:401
pma::Track3D::HasTwoViews
bool HasTwoViews(size_t nmin=1) const
Definition: PmaTrack3D.cxx:443
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
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 TVector3 &  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 930 of file ProjectionMatchingAlg.cxx.

934 {
935  pma::Track3D* trk = buildSegment(detProp, hits_1, hits_2);
936 
937  if (trk) {
938  double dfront = pma::Dist2(trk->front()->Point3D(), point);
939  double dback = pma::Dist2(trk->back()->Point3D(), point);
940  if (dfront > dback) trk->Flip();
941 
942  trk->Nodes().front()->SetPoint3D(point);
943  trk->Nodes().front()->SetFrozen(true);
944  trk->Optimize(detProp, 0, fFineTuningEps);
945 
946  trk->SortHits();
947  }
948  return trk;
949 }
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:338
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
pma::ProjectionMatchingAlg::buildSegment
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
Definition: ProjectionMatchingAlg.cxx:900
pma::Track3D::back
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:106
pma::Track3D::Flip
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:534
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
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 953 of file ProjectionMatchingAlg.cxx.

956 {
957  pma::Track3D* trk = buildSegment(detProp, hits);
958 
959  if (trk) {
960  double dfront = pma::Dist2(trk->front()->Point3D(), point);
961  double dback = pma::Dist2(trk->back()->Point3D(), point);
962  if (dfront > dback) trk->Flip(); // detProp);
963 
964  trk->Nodes().front()->SetPoint3D(point);
965  trk->Nodes().front()->SetFrozen(true);
966  trk->Optimize(detProp, 0, fFineTuningEps);
967 
968  trk->SortHits();
969  }
970  return trk;
971 }
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:338
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
pma::ProjectionMatchingAlg::buildSegment
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
Definition: ProjectionMatchingAlg.cxx:900
pma::Track3D::back
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:106
pma::Track3D::Flip
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:534
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 617 of file ProjectionMatchingAlg.cxx.

620 {
621  double vtxarray[3]{vtx.X(), vtx.Y(), vtx.Z()};
622 
623  if (!fGeom->HasTPC(fGeom->FindTPCAtPosition(vtxarray))) return 0;
624 
625  TVector3 vtxv3(vtx.X(), vtx.Y(), vtx.Z());
626 
627  const size_t tpc = fGeom->FindTPCAtPosition(vtxarray).TPC;
628  const size_t cryo = fGeom->FindCryostatAtPosition(vtxarray);
629  const geo::TPCGeo& tpcgeom = fGeom->Cryostat(cryo).TPC(tpc);
630 
631  // use only hits from tpc where the vtx is
632  std::vector<art::Ptr<recob::Hit>> hitstpc;
633  for (size_t h = 0; h < hits.size(); ++h)
634  if (hits[h]->WireID().TPC == tpc) hitstpc.push_back(hits[h]);
635 
636  if (!hitstpc.size()) return 0;
637 
638  std::vector<art::Ptr<recob::Hit>> hitstrk;
639  size_t view = 0;
640  size_t countviews = 0;
641  while (view < 3) {
642  mf::LogVerbatim("ProjectionMatchinAlg") << "collecting hits from view: " << view;
643  if (!tpcgeom.HasPlane(view)) {
644  ++view;
645  continue;
646  }
647 
648  // select hits only for a single view
649  std::vector<art::Ptr<recob::Hit>> hitsview;
650  for (size_t h = 0; h < hitstpc.size(); ++h)
651  if (hitstpc[h]->WireID().Plane == view) hitsview.push_back(hitstpc[h]);
652  if (!hitsview.size()) {
653  ++view;
654  continue;
655  }
656 
657  // filter our small groups of hits, far from main shower
658  std::vector<art::Ptr<recob::Hit>> hitsfilter;
659  TVector2 proj_pr = pma::GetProjectionToPlane(vtxv3, view, tpc, cryo);
660 
661  mf::LogVerbatim("ProjectionMatchinAlg") << "start filter out: ";
662  FilterOutSmallParts(detProp, 2.0, hitsview, hitsfilter, proj_pr);
663  mf::LogVerbatim("ProjectionMatchingAlg") << "after filter out";
664 
665  for (size_t h = 0; h < hitsfilter.size(); ++h)
666  hitstrk.push_back(hitsfilter[h]);
667 
668  if (hitsfilter.size() >= 5) countviews++;
669 
670  ++view;
671  }
672 
673  if (!hitstrk.size() || (countviews < 2)) {
674  mf::LogWarning("ProjectionMatchinAlg") << "too few hits, segment not built";
675  return 0;
676  }
677 
678  // hits are prepared, finally segment is built
679 
680  pma::Track3D* trk = new pma::Track3D();
681  trk = buildSegment(detProp, hitstrk, vtxv3);
682 
683  // make shorter segment to estimate direction more precise
684  ShortenSeg_(detProp, *trk, tpcgeom);
685 
686  if (trk->size() < 10) {
687  mf::LogWarning("ProjectionMatchingAlg") << "too short segment, delete segment";
688  delete trk;
689  return 0;
690  }
691 
692  return trk;
693 }
geo::TPCGeo::HasPlane
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:175
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
h
while getopts h
Definition: getStatistics.sh:13
geo::GeometryCore::FindTPCAtPosition
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.h:2462
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
pma::ProjectionMatchingAlg::FilterOutSmallParts
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
Definition: ProjectionMatchingAlg.cxx:698
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:321
pma::ProjectionMatchingAlg::buildSegment
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
Definition: ProjectionMatchingAlg.cxx:900
geo::GeometryCore::HasTPC
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
Definition: GeometryCore.h:2405
pma::ProjectionMatchingAlg::ShortenSeg_
void ShortenSeg_(const detinfo::DetectorPropertiesData &detProp, pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
Definition: ProjectionMatchingAlg.cxx:817
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
pma::GetProjectionToPlane
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:273
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
geo::GeometryCore::FindCryostatAtPosition
geo::CryostatID::CryostatID_t FindCryostatAtPosition(geo::Point_t const &worldLoc) const
Returns the index of the cryostat at specified location.
Definition: GeometryCore.cxx:393
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 452 of file ProjectionMatchingAlg.cxx.

455 {
456  pma::Track3D* trk = new pma::Track3D(); // track candidate
457  trk->AddHits(detProp, hits_1);
458  trk->AddHits(detProp, hits_2);
459 
460  mf::LogVerbatim("ProjectionMatchingAlg") << "track size: " << trk->size();
461  std::vector<unsigned int> tpcs = trk->TPCs();
462  for (size_t t = 0; t < tpcs.size(); ++t) {
463  mf::LogVerbatim("ProjectionMatchingAlg") << " tpc:" << tpcs[t];
464  }
465  mf::LogVerbatim("ProjectionMatchingAlg")
466  << " #coll:" << trk->NHits(geo::kZ) << " #ind2:" << trk->NHits(geo::kV)
467  << " #ind1:" << trk->NHits(geo::kU);
468 
469  size_t nSegments = getSegCount_(trk->size());
470  size_t nNodes = (size_t)(nSegments - 1); // n nodes to add
471 
472  mf::LogVerbatim("ProjectionMatchingAlg") << " initialize trk";
473  if (!trk->Initialize(detProp)) {
474  mf::LogWarning("ProjectionMatchingAlg") << " initialization failed, delete trk";
475  delete trk;
476  return 0;
477  }
478 
479  double f = twoViewFraction(*trk);
480  if (f > fMinTwoViewFraction) {
481  double g = 0.0;
482  mf::LogVerbatim("ProjectionMatchingAlg") << " optimize trk (" << nSegments << " seg)";
483  if (nNodes) {
484  g = trk->Optimize(detProp, nNodes, fOptimizationEps, false, true, 25,
485  10); // build nodes
486  mf::LogVerbatim("ProjectionMatchingAlg") << " nodes done, g = " << g;
487  trk->SelectAllHits();
488  }
489  g = trk->Optimize(detProp, 0, fFineTuningEps); // final tuning
490  mf::LogVerbatim("ProjectionMatchingAlg") << " tune done, g = " << g;
491 
492  trk->SortHits();
493  // trk->ShiftEndsToHits(); // not sure if useful already here
494  return trk;
495  }
496  else {
497  mf::LogVerbatim("ProjectionMatchingAlg") << " clusters do not match, f = " << f;
498  delete trk;
499  return 0;
500  }
501 }
pma::Track3D::AddHits
void AddHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
Definition: PmaTrack3D.cxx:404
pma::Track3D::SelectAllHits
bool SelectAllHits()
Definition: PmaTrack3D.cxx:2942
geo::kV
Planes which measure V.
Definition: geo_types.h:130
pma::Track3D::Initialize
bool Initialize(detinfo::DetectorPropertiesData const &detProp, float initEndSegW=0.05F)
Definition: PmaTrack3D.cxx:78
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:132
pma::Track3D::TPCs
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:453
pma::ProjectionMatchingAlg::getSegCount_
static size_t getSegCount_(size_t trk_size)
Definition: ProjectionMatchingAlg.cxx:444
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
geo::kU
Planes which measure U.
Definition: geo_types.h:129
pma::Track3D::NHits
unsigned int NHits(unsigned int view) const
Definition: PmaTrack3D.cxx:423
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
pma::ProjectionMatchingAlg::twoViewFraction
double twoViewFraction(pma::Track3D &trk) const
Definition: ProjectionMatchingAlg.cxx:422
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 1005 of file ProjectionMatchingAlg.cxx.

1014 {
1015  size_t nCloseHits = 0;
1016  int forwWires = 3, backWires = -1;
1017  double xMargin = 0.4;
1018  for (auto h : hits)
1019  if ((view == (int)h->WireID().Plane) && (tpc == h->WireID().TPC) &&
1020  (cryo == h->WireID().Cryostat)) {
1021  bool found = false;
1022  for (size_t ht = 0; ht < trk.size(); ht++)
1023  if (trk[ht]->Hit2DPtr().key() == h.key()) {
1024  found = true;
1025  break;
1026  }
1027  if (found) continue;
1028 
1029  int dw = wdir * (h->WireID().Wire - wire);
1030  if ((dw <= forwWires) && (dw >= backWires)) {
1031  double x = detProp.ConvertTicksToX(h->PeakTime(), view, tpc, cryo);
1032  if (fabs(x - drift_x) < xMargin) nCloseHits++;
1033  }
1034  }
1035  if (nCloseHits > 1)
1036  return false;
1037  else
1038  return true;
1039 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
h
while getopts h
Definition: getStatistics.sh:13
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
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 975 of file ProjectionMatchingAlg.cxx.

979 {
980  pma::Track3D* copy = new pma::Track3D(trk);
981  copy->AddHits(detProp, hits);
982 
983  mf::LogVerbatim("ProjectionMatchingAlg")
984  << "ext. track size: " << copy->size() << " #coll:" << copy->NHits(geo::kZ)
985  << " #ind2:" << copy->NHits(geo::kV) << " #ind1:" << copy->NHits(geo::kU);
986 
987  if (add_nodes) {
988  size_t nSegments = getSegCount_(copy->size());
989  int nNodes = nSegments - copy->Nodes().size() + 1; // n nodes to add
990  if (nNodes < 0) nNodes = 0;
991 
992  if (nNodes) {
993  mf::LogVerbatim("ProjectionMatchingAlg") << " add " << nNodes << " nodes";
994  copy->Optimize(detProp, nNodes, fOptimizationEps);
995  }
996  }
997  double g = copy->Optimize(detProp, 0, fFineTuningEps);
998  mf::LogVerbatim("ProjectionMatchingAlg") << " reopt done, g = " << g;
999 
1000  return copy;
1001 }
pma::Track3D::AddHits
void AddHits(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
Definition: PmaTrack3D.cxx:404
geo::kV
Planes which measure V.
Definition: geo_types.h:130
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:338
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:132
pma::ProjectionMatchingAlg::getSegCount_
static size_t getSegCount_(size_t trk_size)
Definition: ProjectionMatchingAlg.cxx:444
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
geo::kU
Planes which measure U.
Definition: geo_types.h:129
pma::Track3D::NHits
unsigned int NHits(unsigned int view) const
Definition: PmaTrack3D.cxx:423
copy
T copy(T const &v)
Definition: RangeForWrapper_test.cc:98
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
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 698 of file ProjectionMatchingAlg.cxx.

703 {
704  size_t min_size = hits_in.size() / 5;
705  if (min_size < 3) min_size = 3;
706 
707  std::vector<size_t> used;
708  std::vector<std::vector<art::Ptr<recob::Hit>>> sub_groups;
709  std::vector<art::Ptr<recob::Hit>> close_hits;
710 
711  float mindist2 = 99999.99;
712  size_t id_sub_groups_save = 0;
713  size_t id = 0;
714  while (GetCloseHits_(detProp, r2d, hits_in, used, close_hits)) {
715 
716  sub_groups.emplace_back(close_hits);
717 
718  for (size_t h = 0; h < close_hits.size(); ++h) {
719  TVector2 hi_cm = pma::WireDriftToCm(detProp,
720  close_hits[h]->WireID().Wire,
721  close_hits[h]->PeakTime(),
722  close_hits[h]->WireID().Plane,
723  close_hits[h]->WireID().TPC,
724  close_hits[h]->WireID().Cryostat);
725 
726  float dist2 = pma::Dist2(hi_cm, vtx2d);
727  if (dist2 < mindist2) {
728  id_sub_groups_save = id;
729  mindist2 = dist2;
730  }
731  }
732 
733  id++;
734  }
735 
736  for (size_t i = 0; i < sub_groups.size(); ++i) {
737  if (i == id_sub_groups_save) {
738  for (auto h : sub_groups[i])
739  hits_out.push_back(h);
740  }
741  }
742 
743  for (size_t i = 0; i < sub_groups.size(); ++i) {
744  if ((i != id_sub_groups_save) && (hits_out.size() < 10) && (sub_groups[i].size() > min_size))
745  for (auto h : sub_groups[i])
746  hits_out.push_back(h);
747  }
748 }
Plane
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
Definition: rawdigitfilter_icarus.fcl:8
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
pma::WireDriftToCm
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:297
TPC
BEGIN_PROLOG TPC
Definition: sbnd_flashalgo.fcl:14
pma::ProjectionMatchingAlg::GetCloseHits_
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
Definition: ProjectionMatchingAlg.cxx:753
h
while getopts h
Definition: getStatistics.sh:13
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
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 753 of file ProjectionMatchingAlg.cxx.

758 {
759  hits_out.clear();
760 
761  const double gapMargin = 5.0; // can be changed to f(id_tpc1, id_tpc2)
762  size_t idx = 0;
763 
764  while ((idx < hits_in.size()) && Has_(used, idx))
765  idx++;
766 
767  if (idx < hits_in.size()) {
768  hits_out.push_back(hits_in[idx]);
769  used.push_back(idx);
770 
771  unsigned int tpc = hits_in[idx]->WireID().TPC;
772  unsigned int cryo = hits_in[idx]->WireID().Cryostat;
773  unsigned int view = hits_in[idx]->WireID().Plane;
774  double wirePitch = fGeom->TPC(tpc, cryo).Plane(view).WirePitch();
775  double driftPitch = detProp.GetXTicksCoefficient(tpc, cryo);
776 
777  double r2d2 = r2d * r2d;
778  double gapMargin2 = sqrt(2 * gapMargin * gapMargin);
779  gapMargin2 = (gapMargin2 + r2d) * (gapMargin2 + r2d);
780 
781  bool collect = true;
782  while (collect) {
783  collect = false;
784  for (size_t i = 0; i < hits_in.size(); i++)
785  if (!Has_(used, i)) {
786  art::Ptr<recob::Hit> const& hi = hits_in[i];
787  TVector2 hi_cm(wirePitch * hi->WireID().Wire, driftPitch * hi->PeakTime());
788 
789  bool accept = false;
790 
791  for (auto const& ho : hits_out) {
792 
793  TVector2 ho_cm(wirePitch * ho->WireID().Wire, driftPitch * ho->PeakTime());
794  double d2 = pma::Dist2(hi_cm, ho_cm);
795 
796  if (d2 < r2d2) {
797  accept = true;
798  break;
799  }
800  }
801  if (accept) {
802  collect = true;
803  hits_out.push_back(hi);
804  used.push_back(i);
805  }
806  }
807  }
808  return true;
809  }
810  else
811  return false;
812 }
pma::ProjectionMatchingAlg::Has_
bool Has_(const std::vector< size_t > &v, size_t idx) const
Definition: ProjectionMatchingAlg.cxx:876
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
detinfo::DetectorPropertiesData::GetXTicksCoefficient
double GetXTicksCoefficient(int t, int c) const
Definition: DetectorPropertiesData.cc:160
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
geo::PlaneGeo::WirePitch
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
size_t pma::ProjectionMatchingAlg::getSegCount_ ( size_t  trk_size)
staticprivate

Definition at line 444 of file ProjectionMatchingAlg.cxx.

445 {
446  int const nSegments = 0.8 * trk_size / sqrt(trk_size);
447  return std::max(size_t{1}, static_cast<size_t>(nSegments));
448 }
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 1096 of file ProjectionMatchingAlg.cxx.

1100 {
1101  unsigned int tpc = trk.FrontTPC(), cryo = trk.FrontCryo();
1102  if ((tpc != trk.BackTPC()) || (cryo != trk.BackCryo())) {
1103  mf::LogWarning("ProjectionMatchingAlg") << "Please, apply before TPC stitching.";
1104  return;
1105  }
1106 
1107  const double maxCosXZ = 0.992546; // 7 deg
1108 
1109  pma::Segment3D* segFront = trk.Segments().front();
1110  if (trk.Segments().size() > 2) {
1111  pma::Segment3D* segFront1 = trk.Segments()[1];
1112  if ((segFront->Length() < 0.8) && (segFront1->Length() > 5.0)) segFront = segFront1;
1113  }
1114  pma::Vector3D dirFront = segFront->GetDirection3D();
1115  pma::Vector3D dirFrontXZ(dirFront.X(), 0., dirFront.Z());
1116  dirFrontXZ *= 1.0 / dirFrontXZ.R();
1117 
1118  pma::Segment3D* segBack = trk.Segments().back();
1119  if (trk.Segments().size() > 2) {
1120  pma::Segment3D* segBack1 = trk.Segments()[trk.Segments().size() - 2];
1121  if ((segBack->Length() < 0.8) && (segBack1->Length() > 5.0)) segBack = segBack1;
1122  }
1123  pma::Vector3D dirBack = segBack->GetDirection3D();
1124  pma::Vector3D dirBackXZ(dirBack.X(), 0., dirBack.Z());
1125  dirBackXZ *= 1.0 / dirBackXZ.R();
1126 
1127  if ((fabs(dirFrontXZ.Z()) < maxCosXZ) && (fabs(dirBackXZ.Z()) < maxCosXZ)) {
1128  return; // front & back are not parallel to wire planes => exit
1129  }
1130 
1131  unsigned int nPlanesFront = 0, nPlanesBack = 0;
1132  std::pair<int, int> wiresFront[3], wiresBack[3]; // wire index; index direction
1133  double xFront[3], xBack[3];
1134 
1135  for (unsigned int i = 0; i < 3; i++) {
1136  bool frontPresent = false, backPresent = false;
1137  if (fGeom->TPC(tpc, cryo).HasPlane(i)) {
1138  int idxFront0 = trk.NextHit(-1, i);
1139  int idxBack0 = trk.PrevHit(trk.size(), i);
1140  if ((idxFront0 >= 0) && (idxFront0 < (int)trk.size()) && (idxBack0 >= 0) &&
1141  (idxBack0 < (int)trk.size())) {
1142  int idxFront1 = trk.NextHit(idxFront0, i);
1143  int idxBack1 = trk.PrevHit(idxBack0, i);
1144  if ((idxFront1 >= 0) && (idxFront1 < (int)trk.size()) && (idxBack1 >= 0) &&
1145  (idxBack1 < (int)trk.size())) {
1146  int wFront0 = trk[idxFront0]->Wire();
1147  int wBack0 = trk[idxBack0]->Wire();
1148 
1149  int wFront1 = trk[idxFront1]->Wire();
1150  int wBack1 = trk[idxBack1]->Wire();
1151 
1152  wiresFront[i].first = wFront0;
1153  wiresFront[i].second = wFront0 - wFront1;
1154  xFront[i] = detProp.ConvertTicksToX(trk[idxFront0]->PeakTime(), i, tpc, cryo);
1155 
1156  wiresBack[i].first = wBack0;
1157  wiresBack[i].second = wBack0 - wBack1;
1158  xBack[i] = detProp.ConvertTicksToX(trk[idxBack0]->PeakTime(), i, tpc, cryo);
1159 
1160  if (wiresFront[i].second) {
1161  if (wiresFront[i].second > 0)
1162  wiresFront[i].second = 1;
1163  else
1164  wiresFront[i].second = -1;
1165 
1166  frontPresent = true;
1167  nPlanesFront++;
1168  }
1169 
1170  if (wiresBack[i].second) {
1171  if (wiresBack[i].second > 0)
1172  wiresBack[i].second = 1;
1173  else
1174  wiresBack[i].second = -1;
1175 
1176  backPresent = true;
1177  nPlanesBack++;
1178  }
1179  }
1180  }
1181  }
1182  if (!frontPresent) { wiresFront[i].first = -1; }
1183  if (!backPresent) { wiresBack[i].first = -1; }
1184  }
1185 
1186  bool refAdded = false;
1187  if ((nPlanesFront > 1) && (fabs(dirFrontXZ.Z()) >= maxCosXZ)) {
1188  refAdded |= addEndpointRef_(detProp, trk, hits, wiresFront, xFront, tpc, cryo);
1189  }
1190 
1191  if ((nPlanesBack > 1) && (fabs(dirBackXZ.Z()) >= maxCosXZ)) {
1192  refAdded |= addEndpointRef_(detProp, trk, hits, wiresBack, xBack, tpc, cryo);
1193  }
1194  if (refAdded) {
1195  mf::LogVerbatim("ProjectionMatchingAlg") << "guide wire-plane-parallel track endpoints";
1196  double g = trk.Optimize(detProp, 0, 0.1 * fFineTuningEps);
1197  mf::LogVerbatim("ProjectionMatchingAlg") << " done, g = " << g;
1198  }
1199 }
geo::TPCGeo::HasPlane
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:175
pma::ProjectionMatchingAlg::addEndpointRef_
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
Definition: ProjectionMatchingAlg.cxx:1042
second
second
Definition: detectorproperties_icarus.fcl:8
pma::Track3D::PrevHit
int PrevHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:928
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
pma::Track3D::BackTPC
unsigned int BackTPC() const
Definition: PmaTrack3D.h:159
pma::Vector3D
recob::tracking::Vector_t Vector3D
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.h:34
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
pma::Track3D::BackCryo
unsigned int BackCryo() const
Definition: PmaTrack3D.h:164
pma::Track3D::Segments
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:329
pma::Track3D::FrontTPC
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:148
pma::Track3D::FrontCryo
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:153
pma::Track3D::NextHit
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:910
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
pma::Element3D::Length
double Length(void) const
Definition: PmaElement3D.h:53
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
pma::Segment3D::GetDirection3D
pma::Vector3D GetDirection3D(void) const override
Get 3D direction cosines of this segment.
Definition: PmaSegment3D.cxx:61
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 1203 of file ProjectionMatchingAlg.cxx.

1208 {
1209  const double maxCosXZ = 0.992546; // 7 deg
1210 
1211  unsigned int tpc, cryo;
1212  pma::Segment3D* seg0 = 0;
1213  pma::Segment3D* seg1 = 0;
1214 
1215  if (endpoint == pma::Track3D::kBegin) {
1216  tpc = trk.FrontTPC(), cryo = trk.FrontCryo();
1217  seg0 = trk.Segments().front();
1218  if (trk.Segments().size() > 2) { seg1 = trk.Segments()[1]; }
1219  }
1220  else {
1221  tpc = trk.BackTPC(), cryo = trk.BackCryo();
1222  seg0 = trk.Segments().back();
1223  if (trk.Segments().size() > 2) { seg1 = trk.Segments()[trk.Segments().size() - 2]; }
1224  }
1225  if (seg1 && (seg0->Length() < 0.8) && (seg1->Length() > 5.0)) { seg0 = seg1; }
1226  pma::Vector3D dir0 = seg0->GetDirection3D();
1227  pma::Vector3D dir0XZ(dir0.X(), 0., dir0.Z());
1228  dir0XZ *= 1.0 / dir0XZ.R();
1229 
1230  if (fabs(dir0XZ.Z()) < maxCosXZ) { return; } // not parallel to wire planes => exit
1231 
1232  unsigned int nPlanes = 0;
1233  std::pair<int, int> wires[3]; // wire index; index direction
1234  double x0[3];
1235 
1236  for (unsigned int i = 0; i < 3; i++) {
1237  bool present = false;
1238  if (fGeom->TPC(tpc, cryo).HasPlane(i)) {
1239  int idx0 = -1, idx1 = -1;
1240  if (endpoint == pma::Track3D::kBegin) { idx0 = trk.NextHit(-1, i); }
1241  else {
1242  idx0 = trk.PrevHit(trk.size(), i);
1243  }
1244 
1245  if ((idx0 >= 0) && (idx0 < (int)trk.size()) && (trk[idx0]->TPC() == tpc) &&
1246  (trk[idx0]->Cryo() == cryo)) {
1247  if (endpoint == pma::Track3D::kBegin) { idx1 = trk.NextHit(idx0, i); }
1248  else {
1249  idx1 = trk.PrevHit(idx0, i);
1250  }
1251 
1252  if ((idx1 >= 0) && (idx1 < (int)trk.size()) && (trk[idx1]->TPC() == tpc) &&
1253  (trk[idx1]->Cryo() == cryo)) {
1254  int w0 = trk[idx0]->Wire();
1255  int w1 = trk[idx1]->Wire();
1256 
1257  wires[i].first = w0;
1258  wires[i].second = w0 - w1;
1259  x0[i] = detProp.ConvertTicksToX(trk[idx0]->PeakTime(), i, tpc, cryo);
1260 
1261  if (wires[i].second) {
1262  if (wires[i].second > 0)
1263  wires[i].second = 1;
1264  else
1265  wires[i].second = -1;
1266 
1267  present = true;
1268  nPlanes++;
1269  }
1270  }
1271  }
1272  }
1273  if (!present) { wires[i].first = -1; }
1274  }
1275 
1276  if ((nPlanes > 1) && (fabs(dir0XZ.Z()) >= maxCosXZ) &&
1277  addEndpointRef_(detProp, trk, hits, wires, x0, tpc, cryo)) {
1278  mf::LogVerbatim("ProjectionMatchingAlg") << "guide wire-plane-parallel track endpoint";
1279  double g = trk.Optimize(detProp, 0, 0.1 * fFineTuningEps);
1280  mf::LogVerbatim("ProjectionMatchingAlg") << " done, g = " << g;
1281  }
1282 }
geo::TPCGeo::HasPlane
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:175
pma::ProjectionMatchingAlg::addEndpointRef_
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
Definition: ProjectionMatchingAlg.cxx:1042
second
second
Definition: detectorproperties_icarus.fcl:8
pma::Track3D::PrevHit
int PrevHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:928
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
pma::Track3D::BackTPC
unsigned int BackTPC() const
Definition: PmaTrack3D.h:159
pma::Vector3D
recob::tracking::Vector_t Vector3D
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.h:34
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
TPC
BEGIN_PROLOG TPC
Definition: sbnd_flashalgo.fcl:14
pma::Track3D::BackCryo
unsigned int BackCryo() const
Definition: PmaTrack3D.h:164
pma::Track3D::Segments
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:329
pma::Track3D::FrontTPC
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:148
pma::Track3D::FrontCryo
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:153
pma::Track3D::NextHit
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:910
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
pma::Element3D::Length
double Length(void) const
Definition: PmaElement3D.h:53
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
pma::Segment3D::GetDirection3D
pma::Vector3D GetDirection3D(void) const override
Get 3D direction cosines of this segment.
Definition: PmaSegment3D.cxx:61
bool pma::ProjectionMatchingAlg::Has_ ( const std::vector< size_t > &  v,
size_t  idx 
) const
private

Definition at line 876 of file ProjectionMatchingAlg.cxx.

877 {
878  for (auto c : v)
879  if (c == idx) return true;
880  return false;
881 }
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 169 of file ProjectionMatchingAlg.h.

170  {
171  return (trk.FirstElement()->SameTPC(trk.front()->Point3D(), margin) &&
172  trk.LastElement()->SameTPC(trk.back()->Point3D(), margin));
173  }
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::FirstElement
pma::Node3D * FirstElement() const
Definition: PmaTrack3D.h:343
pma::Track3D::LastElement
pma::Node3D * LastElement() const
Definition: PmaTrack3D.h:348
pma::Node3D::SameTPC
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:102
pma::Track3D::back
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:106
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 1335 of file ProjectionMatchingAlg.cxx.

1339 {
1340  if (!alignTracks(dst, src)) return;
1341 
1342  unsigned int tpc = src.FrontTPC();
1343  unsigned int cryo = src.FrontCryo();
1344  double lmean = dst.Length() / (dst.Nodes().size() - 1);
1345  if ((pma::Dist2(dst.Nodes().back()->Point3D(), src.Nodes().front()->Point3D()) > 0.5 * lmean) ||
1346  (tpc != dst.BackTPC()) || (cryo != dst.BackCryo())) {
1347  dst.AddNode(detProp, src.Nodes().front()->Point3D(), tpc, cryo);
1348  if (src.Nodes().front()->IsFrozen()) dst.Nodes().back()->SetFrozen(true);
1349  }
1350  for (size_t n = 1; n < src.Nodes().size(); n++) {
1351  pma::Node3D* node = src.Nodes()[n];
1352 
1353  dst.AddNode(detProp, src.Nodes()[n]->Point3D(), node->TPC(), node->Cryo());
1354 
1355  if (node->IsFrozen()) dst.Nodes().back()->SetFrozen(true);
1356  }
1357  for (size_t h = 0; h < src.size(); h++) {
1358  dst.push_back(detProp, src[h]->Hit2DPtr());
1359  }
1360  if (reopt) {
1361  double g = dst.Optimize(detProp, 0, fFineTuningEps);
1362  mf::LogVerbatim("ProjectionMatchingAlg") << " reopt after merging done, g = " << g;
1363  }
1364  else {
1365  dst.MakeProjection();
1366  }
1367 
1368  dst.SortHits();
1369  dst.ShiftEndsToHits();
1370 
1371  dst.MakeProjection();
1372  dst.SortHits();
1373 }
pma::Track3D::MakeProjection
void MakeProjection()
Definition: PmaTrack3D.cxx:2953
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:338
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
pma::Track3D::BackTPC
unsigned int BackTPC() const
Definition: PmaTrack3D.h:159
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
pma::Element3D::TPC
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
pma::Track3D::BackCryo
unsigned int BackCryo() const
Definition: PmaTrack3D.h:164
h
while getopts h
Definition: getStatistics.sh:13
pma::Track3D::AddNode
void AddNode(pma::Node3D *node)
Definition: PmaTrack3D.cxx:1265
pma::ProjectionMatchingAlg::alignTracks
bool alignTracks(pma::Track3D &first, pma::Track3D &second) const
Definition: ProjectionMatchingAlg.cxx:1293
pma::Track3D::FrontTPC
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:148
pma::Track3D::FrontCryo
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:153
pma::Track3D::push_back
void push_back(pma::Hit3D *hit)
Definition: PmaTrack3D.h:90
pma::Element3D::Cryo
int Cryo(void) const
Cryostat index or -1 if out of any cryostat.
Definition: PmaElement3D.h:37
pma::Track3D::Length
double Length(size_t step=1) const
Definition: PmaTrack3D.h:120
pma::Element3D::IsFrozen
bool IsFrozen(void) const
Check if the vertex 3D position is fixed.
Definition: PmaElement3D.h:100
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
pma::Track3D::ShiftEndsToHits
bool ShiftEndsToHits()
Definition: PmaTrack3D.cxx:2480
void pma::ProjectionMatchingAlg::RemoveNotEnabledHits ( pma::Track3D &  trk) const

Definition at line 886 of file ProjectionMatchingAlg.cxx.

887 {
888  size_t h = 0;
889  while (h < trk.size()) {
890  if (trk[h]->IsEnabled())
891  ++h;
892  else
893  (trk.release_at(h));
894  }
895 }
h
while getopts h
Definition: getStatistics.sh:13
pma::Track3D::release_at
pma::Hit3D * release_at(size_t index)
Definition: PmaTrack3D.cxx:344
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
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 1377 of file ProjectionMatchingAlg.cxx.

1380 {
1381  for (size_t i = 0; i < trk.size(); i++) {
1382  pma::Hit3D* hit = trk[i];
1383  if (hit->View2D() == view) {
1384  if ((hit->GetDistToProj() > 0.5) || // more than 0.5cm away away from the segment
1385  (hit->GetSegFraction() < -1.0)) // projects before segment start (to check!!!)
1386  hit->TagOutlier(true);
1387  else
1388  hit->TagOutlier(false);
1389  }
1390  }
1391 
1392  unsigned int nhits = 0;
1393  double last_x, dx = 0.0, last_q, dq = 0.0, dqdx = 0.0;
1394  int ih = trk.NextHit(-1, view);
1395 
1396  pma::Hit3D* hit = trk[ih];
1397  pma::Hit3D* lastHit = hit;
1398 
1399  if ((ih >= 0) && (ih < (int)trk.size())) {
1400  hit->TagOutlier(true);
1401 
1402  ih = trk.NextHit(ih, view);
1403  while ((dx < 2.5) && (ih >= 0) && (ih < (int)trk.size())) {
1404  hit = trk[ih];
1405 
1406  if (util::absDiff(hit->Wire(), lastHit->Wire()) > 2) break; // break on gap in wire direction
1407 
1408  last_x = trk.HitDxByView(ih, view);
1409  last_q = hit->SummedADC();
1410  if (dx + last_x < 3.0) {
1411  dq += last_q;
1412  dx += last_x;
1413  nhits++;
1414  }
1415  else
1416  break;
1417 
1418  lastHit = hit;
1419  ih = trk.NextHit(ih, view);
1420  }
1421  while ((ih >= 0) && (ih < (int)trk.size())) {
1422  hit = trk[ih];
1423  hit->TagOutlier(true);
1424  ih = trk.NextHit(ih, view);
1425  }
1426  }
1427  else {
1428  mf::LogError("ProjectionMatchingAlg") << "Initial part selection failed.";
1429  }
1430 
1431  if (!nhits) {
1432  mf::LogError("ProjectionMatchingAlg") << "Initial part too short to select useful hits.";
1433  }
1434 
1435  if (dx > 0.0) dqdx = dq / dx;
1436 
1437  if (nused) (*nused) = nhits;
1438 
1439  return dqdx;
1440 }
pma::Hit3D::GetSegFraction
float GetSegFraction() const noexcept
Definition: PmaHit3D.h:143
pma::Hit3D::TagOutlier
void TagOutlier(bool state) noexcept
Definition: PmaHit3D.h:177
hit
process_name hit
Definition: cheaterreco.fcl:51
pma::Hit3D::View2D
unsigned int View2D() const noexcept
Definition: PmaHit3D.h:93
util::absDiff
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:43
pma::Hit3D::GetDistToProj
double GetDistToProj() const
Definition: PmaHit3D.h:136
pma::Track3D::NextHit
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:910
pma::Track3D::HitDxByView
double HitDxByView(size_t index, unsigned int view) const
Definition: PmaTrack3D.cxx:1012
pma::Hit3D::Wire
unsigned int Wire() const noexcept
Definition: PmaHit3D.h:98
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
pma::Hit3D::SummedADC
float SummedADC() const noexcept
Definition: PmaHit3D.h:109
void pma::ProjectionMatchingAlg::ShortenSeg_ ( const detinfo::DetectorPropertiesData &  detProp,
pma::Track3D &  trk,
const geo::TPCGeo &  tpcgeom 
) const
private

Definition at line 817 of file ProjectionMatchingAlg.cxx.

820 {
821  double mse = trk.GetMse();
822  mf::LogWarning("ProjectionMatchingAlg") << "initial value of mse: " << mse;
823  while ((mse > 0.5) && TestTrk_(trk, tpcgeom)) {
824  mse = trk.GetMse();
825  for (size_t h = 0; h < trk.size(); ++h)
826  if (std::sqrt(pma::Dist2(trk[h]->Point3D(), trk.front()->Point3D())) > 0.8 * trk.Length())
827  trk[h]->SetEnabled(false);
828 
829  RemoveNotEnabledHits(trk);
830 
831  // trk.Optimize(0.0001, false); // BUG: first argument missing; tentatively:
832  trk.Optimize(detProp, 0, 0.0001, false);
833  trk.SortHits();
834 
835  mf::LogWarning("ProjectionMatchingAlg") << " mse: " << mse;
836  if (mse == trk.GetMse()) break;
837 
838  mse = trk.GetMse();
839  }
840 
841  RemoveNotEnabledHits(trk);
842 }
pma::Track3D::GetMse
double GetMse(unsigned int view=geo::kUnknown) const
MSE of hits weighted with hit amplidudes and wire plane coefficients.
Definition: PmaTrack3D.cxx:1838
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
lar::fhicl::geo::Point3D
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
Definition: geo_vectors_fhicl.h:38
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::Optimize
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.
Definition: PmaTrack3D.cxx:1869
h
while getopts h
Definition: getStatistics.sh:13
pma::Track3D::Length
double Length(size_t step=1) const
Definition: PmaTrack3D.h:120
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
pma::ProjectionMatchingAlg::RemoveNotEnabledHits
void RemoveNotEnabledHits(pma::Track3D &trk) const
Definition: ProjectionMatchingAlg.cxx:886
pma::ProjectionMatchingAlg::TestTrk_
bool TestTrk_(pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
Definition: ProjectionMatchingAlg.cxx:847
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 158 of file ProjectionMatchingAlg.h.

162  {
163  return trk.TestHits(detProp, hits, eps * fHitTestingDist2D);
164  }
pma::Track3D::TestHits
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:859
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
bool pma::ProjectionMatchingAlg::TestTrk_ ( pma::Track3D &  trk,
const geo::TPCGeo &  tpcgeom 
) const
private

Definition at line 847 of file ProjectionMatchingAlg.cxx.

848 {
849  bool test = false;
850 
851  if (tpcgeom.HasPlane(0) && tpcgeom.HasPlane(1)) {
852  if ((trk.NEnabledHits(0) > 5) && (trk.NEnabledHits(1) > 5)) test = true;
853  }
854  else if (tpcgeom.HasPlane(0) && tpcgeom.HasPlane(2)) {
855  if ((trk.NEnabledHits(0) > 5) && (trk.NEnabledHits(2) > 5)) test = true;
856  }
857  else if (tpcgeom.HasPlane(1) && tpcgeom.HasPlane(2)) {
858  if ((trk.NEnabledHits(1) > 5) && (trk.NEnabledHits(2) > 5)) test = true;
859  }
860 
861  double length = 0.0;
862  for (size_t h = 0; h < trk.size(); ++h) {
863  if (!trk[h]->IsEnabled()) break;
864  length = std::sqrt(pma::Dist2(trk.front()->Point3D(), trk[h]->Point3D()));
865  }
866 
867  mf::LogWarning("ProjectionMatchingAlg") << "length of segment: " << length;
868  if (length < 3.0) test = false; // cm
869 
870  return test;
871 }
geo::TPCGeo::HasPlane
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:175
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::NEnabledHits
unsigned int NEnabledHits(unsigned int view=geo::kUnknown) const
Definition: PmaTrack3D.cxx:433
h
while getopts h
Definition: getStatistics.sh:13
test
BEGIN_PROLOG hitmakerfive clustermakerfour pfparticlemakerthree showermakertwo END_PROLOG hitmakerfive clustermakerfour pfparticlemakerthree sequence::inline_paths sequence::inline_paths sequence::inline_paths showermakers test
Definition: assnschainutils_test.fcl:19
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
std::vector< pma::Hit3D * > pma::ProjectionMatchingAlg::trimTrackToVolume ( pma::Track3D &  trk,
TVector3  p0,
TVector3  p1 
) const

Definition at line 1286 of file ProjectionMatchingAlg.cxx.

1287 {
1288  return {};
1289 }
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 422 of file ProjectionMatchingAlg.cxx.

423 {
424  trk.SelectHits();
425  trk.DisableSingleViewEnds();
426 
427  size_t idx = 0;
428  while ((idx < trk.size() - 1) && !trk[idx]->IsEnabled())
429  idx++;
430  double l0 = trk.Length(0, idx + 1);
431 
432  idx = trk.size() - 1;
433  while ((idx > 1) && !trk[idx]->IsEnabled())
434  idx--;
435  double l1 = trk.Length(idx - 1, trk.size() - 1);
436 
437  trk.SelectHits();
438 
439  return 1.0 - (l0 + l1) / trk.Length();
440 }
pma::Track3D::SelectHits
bool SelectHits(float fraction=1.0F)
Definition: PmaTrack3D.cxx:2889
pma::Track3D::DisableSingleViewEnds
unsigned int DisableSingleViewEnds()
Definition: PmaTrack3D.cxx:2750
pma::Track3D::Length
double Length(size_t step=1) const
Definition: PmaTrack3D.h:120
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:111
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 255 of file ProjectionMatchingAlg.cxx.

259 {
260  if (hits.empty()) { return 0; }
261 
262  double max_d = fTrkValidationDist2D;
263  double d2, max_d2 = max_d * max_d;
264  unsigned int nAll = 0, nPassed = 0;
265  unsigned int testPlane = hits.front()->WireID().Plane;
266 
267  std::vector<unsigned int> trkTPCs = trk.TPCs();
268  std::vector<unsigned int> trkCryos = trk.Cryos();
269  std::map<std::pair<unsigned int, unsigned int>, std::pair<TVector2, TVector2>> ranges;
270  std::map<std::pair<unsigned int, unsigned int>, double> wirePitch;
271  for (auto const& [t, c] : views::cartesian_product(trkTPCs, trkCryos)) {
272  ranges[{t, c}] = trk.WireDriftRange(detProp, testPlane, t, c);
273  wirePitch[{t, c}] = fGeom->TPC(t, c).Plane(testPlane).WirePitch();
274  }
275 
276  unsigned int tpc, cryo;
277  std::map<std::pair<unsigned int, unsigned int>, std::vector<pma::Vector2D>> all_close_points;
278 
279  for (auto const& h :
280  hits | views::transform(to_element) | views::filter(hits_on_plane{testPlane})) {
281  tpc = h.WireID().TPC;
282  cryo = h.WireID().Cryostat;
283  std::pair<unsigned int, unsigned int> tpc_cryo(tpc, cryo);
284  std::pair<TVector2, TVector2> rect = ranges[tpc_cryo];
285 
286  if ((h.WireID().Wire > rect.first.X() - 10) && // chceck only hits in the rectangle around
287  (h.WireID().Wire < rect.second.X() + 10) && // the track projection, it is faster than
288  (h.PeakTime() > rect.first.Y() - 100) && // calculation of trk.Dist2(p2d, testPlane)
289  (h.PeakTime() < rect.second.Y() + 100)) {
290  TVector2 p2d(wirePitch[tpc_cryo] * h.WireID().Wire,
291  detProp.ConvertTicksToX(h.PeakTime(), testPlane, tpc, cryo));
292 
293  d2 = trk.Dist2(p2d, testPlane, tpc, cryo);
294 
295  if (d2 < max_d2) all_close_points[tpc_cryo].emplace_back(p2d.X(), p2d.Y());
296  }
297  }
298 
299  // then check how points-close-to-the-track-projection are distributed along
300  // the track, namely: are there track sections crossing empty spaces, except
301  // dead wires?
302  pma::Vector3D p(
303  trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
304  for (auto const* seg : trk.Segments()) {
305  if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
306  {
307  p = seg->End();
308  continue;
309  }
310  pma::Vector3D p0 = seg->Start();
311  pma::Vector3D p1 = seg->End();
312 
313  pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());
314 
315  tpc = seg->TPC();
316  cryo = seg->Cryo();
317 
318  pma::Vector3D dc = step * seg->GetDirection3D();
319 
320  auto const& points = all_close_points[{tpc, cryo}];
321 
322  double f = pma::GetSegmentProjVector(p, p0, p1);
323 
324  double wirepitch = fGeom->TPC(tpc, cryo).Plane(testPlane).WirePitch();
325  while ((f < 1.0) && node->SameTPC(p)) {
326  pma::Vector2D p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
327  geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
328  if (fGeom->HasWire(wireID)) {
329  raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
330  if (channelStatus.IsGood(ch)) {
331  if (points.size()) {
332  p2d.SetX(wirepitch * p2d.X());
333  for (const auto& h : points) {
334  d2 = pma::Dist2(p2d, h);
335  if (d2 < max_d2) {
336  nPassed++;
337  break;
338  }
339  }
340  }
341  nAll++;
342  }
343  //else mf::LogVerbatim("ProjectionMatchingAlg")
344  // << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
345  }
346 
347  p += dc;
348  f = pma::GetSegmentProjVector(p, p0, p1);
349  }
350  p = seg->End();
351  }
352 
353  if (nAll > 0) {
354  double v = nPassed / (double)nAll;
355  mf::LogVerbatim("ProjectionMatchingAlg") << " trk fraction ok: " << v;
356  return v;
357  }
358 
359  return 1.0;
360 }
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1217
lar::to_element
constexpr to_element_t to_element
Definition: ToElement.h:24
pma::Track3D::Dist2
double Dist2(const TVector2 &p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:2576
icarus::waveform_operations::details::transform
static constexpr Sample_t transform(Sample_t sample)
Definition: WaveformOperations.h:145
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
p
pdgs p
Definition: selectors.fcl:22
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::TPCs
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:453
pma::Vector3D
recob::tracking::Vector_t Vector3D
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.h:34
pma::Element3D::TPC
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
pma::Track3D::Cryos
std::vector< unsigned int > Cryos() const
Definition: PmaTrack3D.cxx:472
h
while getopts h
Definition: getStatistics.sh:13
pma::Track3D::Segments
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:329
lariov::ChannelStatusProvider::IsGood
virtual bool IsGood(raw::ChannelID_t channel) const
Returns whether the specified channel is physical and good.
Definition: ChannelStatusProvider.h:84
pma::Node3D::SameTPC
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:102
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
filter
physics filters filter
Definition: g4_3drift_faketrigger_filter.fcl:8
geo::GeometryCore::PlaneWireToChannel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Definition: GeometryCore.cxx:1315
pma::GetSegmentProjVector
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:120
pma::Vector2D
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.h:33
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
p1
physics associatedGroupsWithLeft p1
Definition: test_associatedgroupswithleft_full.fcl:24
geo::PlaneGeo::WirePitch
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
pma::Track3D::WireDriftRange
std::pair< TVector2, TVector2 > WireDriftRange(detinfo::DetectorPropertiesData const &detProp, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:491
geo::GeometryCore::HasWire
bool HasWire(geo::WireID const &wireid) const
Returns whether we have the specified wire.
Definition: GeometryCore.h:3386
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).

Definition at line 364 of file ProjectionMatchingAlg.cxx.

372 {
373  double max_d = fTrkValidationDist2D;
374  double d2, max_d2 = max_d * max_d;
375  unsigned int nAll = 0, nPassed = 0;
376 
377  TVector3 p(p0);
378  TVector3 dc(p1);
379  dc -= p;
380  dc *= step / dc.Mag();
381 
382  double f = pma::GetSegmentProjVector(p, p0, p1);
383  double wirepitch = fGeom->TPC(tpc, cryo).Plane(testPlane).WirePitch();
384  while (f < 1.0) {
385  TVector2 p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
386  geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
387  if (fGeom->HasWire(wireID)) {
388  raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
389  if (channelStatus.IsGood(ch)) {
390  p2d.Set(wirepitch * p2d.X(), p2d.Y());
391  for (const auto& h : hits)
392  if ((h->WireID().Plane == testPlane) && (h->WireID().TPC == tpc) &&
393  (h->WireID().Cryostat == cryo)) {
394  d2 = pma::Dist2(
395  p2d,
396  pma::WireDriftToCm(detProp, h->WireID().Wire, h->PeakTime(), testPlane, tpc, cryo));
397  if (d2 < max_d2) {
398  nPassed++;
399  break;
400  }
401  }
402  nAll++;
403  }
404  }
405 
406  p += dc;
407  f = pma::GetSegmentProjVector(p, p0, p1);
408  }
409 
410  if (nAll > 3) // validate actually only if 2D projection in testPlane has some minimum length
411  {
412  double v = nPassed / (double)nAll;
413  mf::LogVerbatim("ProjectionMatchingAlg") << " segment fraction ok: " << v;
414  return v;
415  }
416  else
417  return 1.0;
418 }
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1217
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
p
pdgs p
Definition: selectors.fcl:22
pma::WireDriftToCm
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:297
h
while getopts h
Definition: getStatistics.sh:13
lariov::ChannelStatusProvider::IsGood
virtual bool IsGood(raw::ChannelID_t channel) const
Returns whether the specified channel is physical and good.
Definition: ChannelStatusProvider.h:84
geo::GeometryCore::PlaneWireToChannel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Definition: GeometryCore.cxx:1315
pma::GetSegmentProjVector
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:120
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
p1
physics associatedGroupsWithLeft p1
Definition: test_associatedgroupswithleft_full.fcl:24
geo::PlaneGeo::WirePitch
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
geo::GeometryCore::HasWire
bool HasWire(geo::WireID const &wireid) const
Returns whether we have the specified wire.
Definition: GeometryCore.h:3386
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 128 of file ProjectionMatchingAlg.cxx.

135 {
136  double max_d = fTrkValidationDist2D;
137  double d2, max_d2 = max_d * max_d;
138  unsigned int nAll = 0, nPassed = 0;
139  unsigned int testPlane = adcImage.Plane();
140 
141  std::vector<unsigned int> trkTPCs = trk.TPCs();
142  std::vector<unsigned int> trkCryos = trk.Cryos();
143  std::map<std::pair<unsigned int, unsigned int>, std::pair<TVector2, TVector2>> ranges;
144  std::map<std::pair<unsigned int, unsigned int>, double> wirePitch;
145  for (auto const& [t, c] : views::cartesian_product(trkTPCs, trkCryos)) {
146  ranges[{t, c}] = trk.WireDriftRange(detProp, testPlane, t, c);
147  wirePitch[{t, c}] = fGeom->TPC(t, c).Plane(testPlane).WirePitch();
148  }
149 
150  unsigned int tpc, cryo;
151  std::map<std::pair<unsigned int, unsigned int>, std::vector<pma::Vector2D>> all_close_points;
152 
153  for (auto const& h :
154  hits | views::transform(to_element) | views::filter(hits_on_plane{testPlane})) {
155  tpc = h.WireID().TPC;
156  cryo = h.WireID().Cryostat;
157  std::pair<unsigned int, unsigned int> tpc_cryo(tpc, cryo);
158  std::pair<TVector2, TVector2> rect = ranges[tpc_cryo];
159 
160  if ((h.WireID().Wire > rect.first.X() - 10) && // check only hits in the rectangle around
161  (h.WireID().Wire < rect.second.X() + 10) && // the track projection, it is faster than
162  (h.PeakTime() > rect.first.Y() - 100) && // calculation of trk.Dist2(p2d, testPlane)
163  (h.PeakTime() < rect.second.Y() + 100)) {
164  TVector2 p2d(wirePitch[tpc_cryo] * h.WireID().Wire,
165  detProp.ConvertTicksToX(h.PeakTime(), testPlane, tpc, cryo));
166 
167  d2 = trk.Dist2(p2d, testPlane, tpc, cryo);
168 
169  if (d2 < max_d2) { all_close_points[tpc_cryo].emplace_back(p2d.X(), p2d.Y()); }
170  }
171  }
172 
173  // then check how points-close-to-the-track-projection are distributed along
174  // the track, namely: are there track sections crossing empty spaces, except
175  // dead wires?
176  pma::Vector3D p(
177  trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
178  for (auto const* seg : trk.Segments()) {
179  if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
180  {
181  p = seg->End();
182  continue;
183  }
184  pma::Vector3D p0 = seg->Start();
185  pma::Vector3D p1 = seg->End();
186 
187  pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());
188 
189  tpc = seg->TPC();
190  cryo = seg->Cryo();
191 
192  pma::Vector3D dc = step * seg->GetDirection3D();
193 
194  auto const& points = all_close_points[{tpc, cryo}];
195 
196  double f = pma::GetSegmentProjVector(p, p0, p1);
197 
198  double wirepitch = fGeom->TPC(tpc, cryo).Plane(testPlane).WirePitch();
199  while ((f < 1.0) && node->SameTPC(p)) {
200  pma::Vector2D p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
201  geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
202 
203  int widx = (int)p2d.X();
204  int didx = (int)detProp.ConvertXToTicks(p2d.Y(), testPlane, tpc, cryo);
205 
206  if (fGeom->HasWire(wireID)) {
207  raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
208  if (channelStatus.IsGood(ch)) {
209  bool is_close = false;
210  float max_adc = adcImage.poolMax(widx, didx, 2);
211 
212  if (points.size()) {
213  p2d.SetX(wirepitch * p2d.X());
214  for (const auto& h : points) {
215  d2 = pma::Dist2(p2d, h);
216  if (d2 < max_d2) {
217  is_close = true;
218  nPassed++;
219  break;
220  }
221  }
222  }
223  nAll++;
224 
225  // now fill the calibration histograms
226  if (is_close) {
227  if (histoPassing) histoPassing->Fill(max_adc);
228  }
229  else {
230  if (histoRejected) histoRejected->Fill(max_adc);
231  }
232  }
233  //else mf::LogVerbatim("ProjectionMatchingAlg")
234  // << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
235  }
236 
237  p += dc;
238  f = pma::GetSegmentProjVector(p, p0, p1);
239  }
240  p = seg->End();
241  }
242 
243  if (nAll > 0) {
244  double v = nPassed / (double)nAll;
245  mf::LogVerbatim("ProjectionMatchingAlg") << " trk fraction ok: " << v;
246  return v;
247  }
248 
249  return 1.0;
250 }
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1217
lar::to_element
constexpr to_element_t to_element
Definition: ToElement.h:24
pma::Track3D::Dist2
double Dist2(const TVector2 &p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:2576
icarus::waveform_operations::details::transform
static constexpr Sample_t transform(Sample_t sample)
Definition: WaveformOperations.h:145
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
pma::ProjectionMatchingAlg::fGeom
geo::GeometryCore const * fGeom
Definition: ProjectionMatchingAlg.h:338
pma::Track3D::front
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:101
p
pdgs p
Definition: selectors.fcl:22
pma::Hit3D::Point3D
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
pma::Track3D::TPCs
std::vector< unsigned int > TPCs() const
Definition: PmaTrack3D.cxx:453
img::DataProviderAlg::Plane
unsigned int Plane() const
Definition: DataProviderAlg.h:183
pma::Vector3D
recob::tracking::Vector_t Vector3D
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.h:34
pma::Element3D::TPC
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
pma::Track3D::Cryos
std::vector< unsigned int > Cryos() const
Definition: PmaTrack3D.cxx:472
h
while getopts h
Definition: getStatistics.sh:13
pma::Track3D::Segments
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:329
detinfo::DetectorPropertiesData::ConvertXToTicks
double ConvertXToTicks(double X, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:99
lariov::ChannelStatusProvider::IsGood
virtual bool IsGood(raw::ChannelID_t channel) const
Returns whether the specified channel is physical and good.
Definition: ChannelStatusProvider.h:84
pma::Node3D::SameTPC
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:102
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
filter
physics filters filter
Definition: g4_3drift_faketrigger_filter.fcl:8
geo::GeometryCore::PlaneWireToChannel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Definition: GeometryCore.cxx:1315
pma::GetSegmentProjVector
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:120
pma::Vector2D
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.h:33
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
img::DataProviderAlg::poolMax
float poolMax(int wire, int drift, size_t r=0) const
Pool max value in a patch around the wire/drift pixel.
Definition: DataProviderAlg.cxx:146
p1
physics associatedGroupsWithLeft p1
Definition: test_associatedgroupswithleft_full.fcl:24
geo::PlaneGeo::WirePitch
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
pma::Track3D::WireDriftRange
std::pair< TVector2, TVector2 > WireDriftRange(detinfo::DetectorPropertiesData const &detProp, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:491
geo::GeometryCore::HasWire
bool HasWire(geo::WireID const &wireid) const
Returns whether we have the specified wire.
Definition: GeometryCore.h:3386

Member Data Documentation

double const pma::ProjectionMatchingAlg::fFineTuningEps
private

Definition at line 326 of file ProjectionMatchingAlg.h.

geo::GeometryCore const* pma::ProjectionMatchingAlg::fGeom
private

Definition at line 338 of file ProjectionMatchingAlg.h.

double const pma::ProjectionMatchingAlg::fHitTestingDist2D
private

Definition at line 331 of file ProjectionMatchingAlg.h.

double const pma::ProjectionMatchingAlg::fMinTwoViewFraction
private

Definition at line 334 of file ProjectionMatchingAlg.h.

double const pma::ProjectionMatchingAlg::fOptimizationEps
private

Definition at line 322 of file ProjectionMatchingAlg.h.

double const pma::ProjectionMatchingAlg::fTrkValidationDist2D
private

Definition at line 329 of file ProjectionMatchingAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5