trkf::SpacePointAlg Class Reference

#include "SpacePointAlg.h"

Classes
struct	HitMCInfo

Public Member Functions
	SpacePointAlg (const fhicl::ParameterSet &pset)
bool	filter () const noexcept
bool	merge () const noexcept
double	maxDT () const noexcept
double	maxS () const noexcept
int	minViews () const noexcept
bool	enableU () const noexcept
bool	enableV () const noexcept
bool	enableW () const noexcept
void	update (detinfo::DetectorPropertiesData const &detProp) const
double	correctedTime (detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
double	separation (const art::PtrVector< recob::Hit > &hits) const
bool	compatible (detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
void	fillSpacePoint (detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
void	fillSpacePoints (detinfo::DetectorPropertiesData const &detProp, std::vector< recob::SpacePoint > &spts, std::multimap< double, KHitTrack > const &trackMap) const
	Fill a collection of space points. More...
void	fillComplexSpacePoint (detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
void	makeSpacePoints (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
void	makeMCTruthSpacePoints (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
const art::PtrVector< recob::Hit > &	getAssociatedHits (const recob::SpacePoint &spt) const
void	clearHitMap () const
int	numHitMap () const

Private Member Functions
void	makeSpacePoints (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts, bool useMC) const

Private Attributes
double	fMaxDT
	Maximum time difference between planes. More...
double	fMaxS
	Maximum space separation between wires. More...
int	fMinViews
	Mininum number of views per space point. More...
bool	fEnableU
	Enable flag (U). More...
bool	fEnableV
	Enable flag (V). More...
bool	fEnableW
	Enable flag (W). More...
bool	fFilter
	Filter flag. More...
bool	fMerge
	Merge flag. More...
bool	fPreferColl
	Sort by collection wire. More...
double	fTickOffsetU
	Tick offset for plane U. More...
double	fTickOffsetV
	Tick offset for plane V. More...
double	fTickOffsetW
	Tick offset for plane W. More...
std::map< const recob::Hit *, HitMCInfo >	fHitMCMap
std::map< int, art::PtrVector< recob::Hit > >	fSptHitMap

Detailed Description

Definition at line 80 of file SpacePointAlg.h.

Constructor & Destructor Documentation

trkf::SpacePointAlg::SpacePointAlg

(

const fhicl::ParameterSet &

pset

)

Definition at line 47 of file SpacePointAlg.cxx.

References fEnableU, fEnableV, fEnableW, fFilter, fMaxDT, fMaxS, fMerge, fMinViews, fPreferColl, fTickOffsetU, fTickOffsetV, fTickOffsetW, and fhicl::ParameterSet::get().

    : fMaxDT{pset.get<double>("MaxDT")}
    , fMaxS{pset.get<double>("MaxS")}
    , fMinViews{pset.get<int>("MinViews")}
    , fEnableU{pset.get<bool>("EnableU")}
    , fEnableV{pset.get<bool>("EnableV")}
    , fEnableW{pset.get<bool>("EnableW")}
    , fFilter{pset.get<bool>("Filter")}
    , fMerge{pset.get<bool>("Merge")}
    , fPreferColl{pset.get<bool>("PreferColl")}
    , fTickOffsetU{pset.get<double>("TickOffsetU", 0.)}
    , fTickOffsetV{pset.get<double>("TickOffsetV", 0.)}
    , fTickOffsetW{pset.get<double>("TickOffsetW", 0.)}
  {
    // Only allow one of fFilter and fMerge to be true.

    if (fFilter && fMerge)
      throw cet::exception("SpacePointAlg") << "Filter and Merge flags are both true.\n";

    // Report.

    std::cout << "SpacePointAlg configured with the following parameters:\n"
              << "  MaxDT = " << fMaxDT << "\n"
              << "  MaxS = " << fMaxS << "\n"
              << "  MinViews = " << fMinViews << "\n"
              << "  EnableU = " << fEnableU << "\n"
              << "  EnableV = " << fEnableV << "\n"
              << "  EnableW = " << fEnableW << "\n"
              << "  Filter = " << fFilter << "\n"
              << "  Merge = " << fMerge << "\n"
              << "  PreferColl = " << fPreferColl << "\n"
              << "  TickOffsetU = " << fTickOffsetU << "\n"
              << "  TickOffsetV = " << fTickOffsetV << "\n"
              << "  TickOffsetW = " << fTickOffsetW << std::endl;
  }

Member Function Documentation

void trkf::SpacePointAlg::clearHitMap ( ) const

inline

Definition at line 149 of file SpacePointAlg.h.

Referenced by trkf::SeedFinderAlgorithm::FindSeeds(), cluster::ClusterMatchAlg::Match_SpacePoint(), trkf::Track3DKalmanHit::prepareForInput(), trkf::TrackKalmanCheater::produce(), and trkf::Track3DKalmanHit::produce().

{ fSptHitMap.clear(); }

bool trkf::SpacePointAlg::compatible	(	detinfo::DetectorPropertiesData const &	detProp,
		const art::PtrVector< recob::Hit > &	hits,
		bool	useMC = false
	)		const

Definition at line 260 of file SpacePointAlg.cxx.

References util::abs(), geo::CryostatID::Cryostat, larg4::dist(), fHitMCMap, fMaxDT, fMaxS, Get, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), trkf::SpacePointAlg::HitMCInfo::pchit, recob::Hit::PeakTime(), geo::PlaneID::Plane, art::PtrVector< T >::size(), t1, t2, geo::TPCID::TPC, trkf::SpacePointAlg::HitMCInfo::trackIDs, recob::Hit::View(), and recob::Hit::WireID().

Referenced by fillSpacePoints(), and makeSpacePoints().

  {
    art::ServiceHandle<geo::Geometry const> geom;

    int nhits = hits.size();

    // Fewer than two or more than three hits can never be compatible.

    bool result = nhits >= 2 && nhits <= 3;
    bool mc_ok = true;
    unsigned int tpc = 0;
    unsigned int cstat = 0;

    if (result) {

      // First do pairwise tests.
      // Do double loop over hits.

      for (int ihit1 = 0; result && ihit1 < nhits - 1; ++ihit1) {
        const recob::Hit& hit1 = *(hits[ihit1]);
        geo::WireID hit1WireID = hit1.WireID();
        geo::View_t view1 = hit1.View();

        double t1 = hit1.PeakTime() -
                    detProp.GetXTicksOffset(hit1WireID.Plane, hit1WireID.TPC, hit1WireID.Cryostat);

        // If using mc information, get a collection of track ids for hit 1.
        // If not using mc information, this section of code will trigger the
        // insertion of a single invalid HitMCInfo object into fHitMCMap.

        const HitMCInfo& mcinfo1 = fHitMCMap[(useMC ? &hit1 : 0)];
        const std::vector<int>& tid1 = mcinfo1.trackIDs;
        bool only_neg1 = tid1.size() > 0 && tid1.back() < 0;

        // Loop over second hit.

        for (int ihit2 = ihit1 + 1; result && ihit2 < nhits; ++ihit2) {
          const recob::Hit& hit2 = *(hits[ihit2]);
          geo::WireID hit2WireID = hit2.WireID();
          geo::View_t view2 = hit2.View();

          // Test for same tpc and different views.

          result = result && hit1WireID.TPC == hit2WireID.TPC && view1 != view2 &&
                   hit1WireID.Cryostat == hit2WireID.Cryostat;
          if (result) {

            // Remember which tpc and cryostat we are in.

            tpc = hit1WireID.TPC;
            cstat = hit1WireID.Cryostat;

            double t2 = hit2.PeakTime() - detProp.GetXTicksOffset(
                                            hit2WireID.Plane, hit2WireID.TPC, hit2WireID.Cryostat);

            // Test maximum time difference.

            result = result && std::abs(t1 - t2) <= fMaxDT;

            // Test mc truth.

            if (result && useMC) {

              // Test whether hits have a common parent track id.

              const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
              std::vector<int> tid2 = mcinfo2.trackIDs;
              bool only_neg2 = tid2.size() > 0 && tid2.back() < 0;
              std::vector<int>::iterator it = std::set_intersection(
                tid1.begin(), tid1.end(), tid2.begin(), tid2.end(), tid2.begin());
              tid2.resize(it - tid2.begin());

              // Hits are compatible if they have parents in common.
              // If the only parent id in common is negative (-999),
              // then hits are compatible only if both hits have only
              // negative parent tracks.

              bool only_neg3 = tid2.size() > 0 && tid2.back() < 0;
              mc_ok = tid2.size() > 0 && (!only_neg3 || (only_neg1 && only_neg2));
              result = result && mc_ok;

              // If we are still OK, check that either hit is
              // the nearest neighbor of the other.

              if (result) {
                result = mcinfo1.pchit[hit2WireID.Plane] == &hit2 ||
                         mcinfo2.pchit[hit1WireID.Plane] == &hit1;
              }
            }
          }
        }
      }

      // If there are exactly three hits, and they pass pairwise tests, check
      // for spatial compatibility.

      if (result && nhits == 3) {

        // Loop over hits.

        double dist[3] = {0., 0., 0.};
        double sinth[3] = {0., 0., 0.};
        double costh[3] = {0., 0., 0.};

        auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

        for (int i = 0; i < 3; ++i) {

          // Get tpc, plane, wire.

          const recob::Hit& hit = *(hits[i]);
          geo::WireID hitWireID = hit.WireID();

          const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());
          if ((hitWireID.TPC != tpc) || (hitWireID.Cryostat != cstat))
            throw cet::exception("SpacePointAlg") << "compatible(): geometry mismatch\n";

          // Get angles and distance of wire.

          double const hl = wgeom.HalfL();
          auto const xyz = wgeom.GetCenter();
          auto const xyz1 = wgeom.GetEnd();
          double s = (xyz1.Y() - xyz.Y()) / hl;
          double c = (xyz1.Z() - xyz.Z()) / hl;
          sinth[hit.WireID().Plane] = s;
          costh[hit.WireID().Plane] = c;
          dist[hit.WireID().Plane] = xyz.Z() * s - xyz.Y() * c;
        }

        // Do space cut.

        double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0] +
                    (sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1] +
                    (sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);

        result = result && std::abs(S) < fMaxS;
      }
    }

    // Done.

    return result;
  }

double trkf::SpacePointAlg::correctedTime	(	detinfo::DetectorPropertiesData const &	detProp,
		const recob::Hit &	hit
	)		const

Definition at line 161 of file SpacePointAlg.cxx.

References fTickOffsetU, fTickOffsetV, fTickOffsetW, detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::kU, geo::kV, geo::kW, recob::Hit::PeakTime(), recob::Hit::View(), and recob::Hit::WireID().

Referenced by trkf::SpacePointAna::analyze().

  {
    // Correct time for trigger offset and plane-dependent time offsets.

    double t = hit.PeakTime() - detProp.GetXTicksOffset(hit.WireID());
    if (hit.View() == geo::kU)
      t -= fTickOffsetU;
    else if (hit.View() == geo::kV)
      t -= fTickOffsetV;
    else if (hit.View() == geo::kW)
      t -= fTickOffsetW;

    return t;
  }

bool trkf::SpacePointAlg::enableU ( ) const

inlinenoexcept

Definition at line 91 of file SpacePointAlg.h.

Referenced by trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits(), trkf::SpacePointFinder::produce(), and trkf::SpacePointCheater::produce().

{ return fEnableU; }

bool trkf::SpacePointAlg::enableV ( ) const

inlinenoexcept

Definition at line 92 of file SpacePointAlg.h.

Referenced by trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits(), trkf::SpacePointFinder::produce(), and trkf::SpacePointCheater::produce().

{ return fEnableV; }

bool trkf::SpacePointAlg::enableW ( ) const

inlinenoexcept

Definition at line 93 of file SpacePointAlg.h.

References hits().

Referenced by trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits(), trkf::SpacePointFinder::produce(), and trkf::SpacePointCheater::produce().

{ return fEnableW; }

void trkf::SpacePointAlg::fillComplexSpacePoint	(	detinfo::DetectorPropertiesData const &	detProp,
		const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	sptv,
		int	sptid
	)		const

Definition at line 608 of file SpacePointAlg.cxx.

References art::PtrVector< T >::begin(), geo::CryostatID::Cryostat, art::PtrVector< T >::end(), art::PtrVector< T >::front(), fSptHitMap, Get, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), detinfo::DetectorPropertiesData::GetXTicksCoefficient(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), hits(), recob::Hit::PeakTime(), geo::PlaneID::Plane, recob::Hit::SigmaPeakTime(), art::PtrVector< T >::size(), ssc, lar::to_element, geo::TPCID::TPC, w, weight, and recob::Hit::WireID().

Referenced by makeSpacePoints().

  {
    art::ServiceHandle<geo::Geometry const> geom;
    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

    // Calculate time pitch.

    double timePitch = detProp.GetXTicksCoefficient(); // cm / tick

    // Figure out which tpc we are in.

    unsigned int tpc0 = 0;
    unsigned int cstat0 = 0;
    int nhits = hits.size();
    if (nhits > 0) {
      tpc0 = hits.front()->WireID().TPC;
      cstat0 = hits.front()->WireID().Cryostat;
    }

    // Remember associated hits internally.

    if (fSptHitMap.count(sptid) != 0)
      throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): hit already present!\n";
    fSptHitMap[sptid] = hits;

    // Do a preliminary scan of hits.
    // Determine weight given to hits in each view.

    unsigned int nplanes = wireReadoutGeom.Nplanes(geo::TPCID(cstat0, tpc0));
    std::vector<int> numhits(nplanes, 0);
    std::vector<double> weight(nplanes, 0.);

    for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
         ++ihit) {

      const recob::Hit& hit = **ihit;
      geo::WireID hitWireID = hit.WireID();
      // kept as assertions for performance reasons
      assert(hitWireID.Cryostat == cstat0);
      assert(hitWireID.TPC == tpc0);
      assert(hitWireID.Plane < nplanes);
      ++numhits[hitWireID.Plane];
    }

    for (unsigned int plane = 0; plane < nplanes; ++plane) {
      double np = numhits[plane];
      if (np > 0.) weight[plane] = 1. / (np * np * np);
    }

    // Calculate position and error matrix.

    double xyz[3] = {0., 0., 0.};
    double errxyz[6] = {0., 0., 0., 0., 0., 0.};

    // Calculate x using drift times.
    // Loop over all hits and calculate the weighted average drift time.
    // Also calculate time variance and chisquare.

    double sumt2w = 0.;
    double sumtw = 0.;
    double sumw = 0.;

    for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
         ++ihit) {

      const recob::Hit& hit = **ihit;
      geo::WireID hitWireID = hit.WireID();

      // Correct time for trigger offset and view-dependent time offsets.

      double t0 = detProp.GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
      double t = hit.PeakTime() - t0;
      double et = hit.SigmaPeakTime();
      double w = weight[hitWireID.Plane] / (et * et);

      sumt2w += w * t * t;
      sumtw += w * t;
      sumw += w;
    }

    double drift_time = 0.;
    double var_time = 0.;
    double chisq = 0.;
    if (sumw != 0.) {
      drift_time = sumtw / sumw;
      var_time = sumt2w / sumw - drift_time * drift_time;
      if (var_time < 0.) var_time = 0.;
      chisq = sumt2w - sumtw * drift_time;
      if (chisq < 0.) chisq = 0.;
    }
    xyz[0] = drift_time * timePitch;
    errxyz[0] = var_time * timePitch * timePitch;

    // Calculate y, z using wires (need at least two hits).

    if (nhits >= 2) {

      // Calculate y and z by chisquare minimization of wire coordinates.

      double sus = 0.; // sum w_i u_i sin_th_i
      double suc = 0.; // sum w_i u_i cos_th_i
      double sc2 = 0.; // sum w_i cos2_th_i
      double ss2 = 0.; // sum w_i sin2_th_i
      double ssc = 0.; // sum w_i sin_th_i cos_th_i

      // Loop over points.

      auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
      for (auto const& hit : hits | ranges::views::transform(to_element)) {
        geo::WireID hitWireID = hit.WireID();
        const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());

        // Calculate angle and wire coordinate in this view.

        double const hl = wgeom.HalfL();
        auto const cen = wgeom.GetCenter();
        auto const cen1 = wgeom.GetEnd();
        double s = (cen1.Y() - cen.Y()) / hl;
        double c = (cen1.Z() - cen.Z()) / hl;
        double u = cen.Z() * s - cen.Y() * c;
        double eu = wireReadoutGeom.Plane(hitWireID).WirePitch() / std::sqrt(12.);
        double w = weight[hitWireID.Plane] / (eu * eu);

        // Summations

        sus += w * u * s;
        suc += w * u * c;
        sc2 += w * c * c;
        ss2 += w * s * s;
        ssc += w * s * c;
      }

      // Calculate y,z

      double denom = sc2 * ss2 - ssc * ssc;
      if (denom != 0.) {
        xyz[1] = (-suc * ss2 + sus * ssc) / denom;
        xyz[2] = (sus * sc2 - suc * ssc) / denom;
        errxyz[2] = ss2 / denom;
        errxyz[4] = ssc / denom;
        errxyz[5] = sc2 / denom;
      }

      // Make space point.

      recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
      sptv.push_back(spt);
    }
  }

void trkf::SpacePointAlg::fillSpacePoint	(	detinfo::DetectorPropertiesData const &	detProp,
		const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	sptv,
		int	sptid
	)		const

Definition at line 410 of file SpacePointAlg.cxx.

References art::PtrVector< T >::begin(), geo::CryostatID::Cryostat, art::PtrVector< T >::end(), fSptHitMap, Get, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), detinfo::DetectorPropertiesData::GetXTicksCoefficient(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), hits(), recob::Hit::PeakTime(), geo::PlaneID::Plane, recob::Hit::SigmaPeakTime(), art::PtrVector< T >::size(), ssc, lar::to_element, geo::TPCID::TPC, w, and recob::Hit::WireID().

Referenced by fillSpacePoints(), and makeSpacePoints().

  {
    art::ServiceHandle<geo::Geometry const> geom;

    double timePitch = detProp.GetXTicksCoefficient();

    int nhits = hits.size();

    // Remember associated hits internally.

    if (fSptHitMap.find(sptid) != fSptHitMap.end())
      throw cet::exception("SpacePointAlg") << "fillSpacePoint(): hit already present!\n";
    fSptHitMap[sptid] = hits;

    // Calculate position and error matrix.

    double xyz[3] = {0., 0., 0.};
    double errxyz[6] = {0., 0., 0., 0., 0., 0.};

    // Calculate x using drift times.
    // Loop over all hits and calculate the weighted average drift time.
    // Also calculate time variance and chisquare.

    double sumt2w = 0.;
    double sumtw = 0.;
    double sumw = 0.;

    for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
         ++ihit) {

      const recob::Hit& hit = **ihit;
      geo::WireID hitWireID = hit.WireID();

      // Correct time for trigger offset and view-dependent time offsets.

      double t0 = detProp.GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
      double t = hit.PeakTime() - t0;
      double et = hit.SigmaPeakTime();
      double w = 1. / (et * et);

      sumt2w += w * t * t;
      sumtw += w * t;
      sumw += w;
    }

    double drift_time = 0.;
    double var_time = 0.;
    double chisq = 0.;
    if (sumw != 0.) {
      drift_time = sumtw / sumw;
      //var_time = sumt2w / sumw - drift_time * drift_time;
      var_time = 1. / sumw;
      if (var_time < 0.) var_time = 0.;
      chisq = sumt2w - sumtw * drift_time;
      if (chisq < 0.) chisq = 0.;
    }
    xyz[0] = drift_time * timePitch;
    errxyz[0] = var_time * timePitch * timePitch;

    // Calculate y, z using wires (need at least two hits).

    if (nhits >= 2) {

      // Calculate y and z by chisquare minimization of wire coordinates.

      double sus = 0.; // sum w_i u_i sin_th_i
      double suc = 0.; // sum w_i u_i cos_th_i
      double sc2 = 0.; // sum w_i cos2_th_i
      double ss2 = 0.; // sum w_i sin2_th_i
      double ssc = 0.; // sum w_i sin_th_i cos_th_i

      // Loop over points.

      auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
      for (auto const& hit : hits | ranges::views::transform(to_element)) {
        geo::WireID hitWireID = hit.WireID();
        const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());

        // Calculate angle and wire coordinate in this view.

        double const hl = wgeom.HalfL();
        auto const cen = wgeom.GetCenter();
        auto const cen1 = wgeom.GetEnd();
        double s = (cen1.Y() - cen.Y()) / hl;
        double c = (cen1.Z() - cen.Z()) / hl;
        double u = cen.Z() * s - cen.Y() * c;
        double eu = wireReadoutGeom.Plane(hitWireID).WirePitch() / std::sqrt(12.);
        double w = 1. / (eu * eu);

        // Summations

        sus += w * u * s;
        suc += w * u * c;
        sc2 += w * c * c;
        ss2 += w * s * s;
        ssc += w * s * c;
      }

      // Calculate y,z

      double denom = sc2 * ss2 - ssc * ssc;
      if (denom != 0.) {
        xyz[1] = (-suc * ss2 + sus * ssc) / denom;
        xyz[2] = (sus * sc2 - suc * ssc) / denom;
        errxyz[2] = ss2 / denom;
        errxyz[4] = ssc / denom;
        errxyz[5] = sc2 / denom;
      }

      // Make space point.

      recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
      sptv.push_back(spt);
    }
    return;
  }

void trkf::SpacePointAlg::fillSpacePoints	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< recob::SpacePoint > &	spts,
		std::multimap< double, KHitTrack > const &	trackMap
	)		const

Fill a collection of space points.

Fill a collection of space points.

Arguments:

spts - Collection of space points to fill. sptalg - Space point algorithm object.

This method uses the hits contained in this track to construct space points.

This method does not have any knowledge of what constitutes a good space point, except that Hits are required to be consecutive when sorted by path distance, and space points are required to pass compatibility tests used by the space point algorithm object. This method will make space points from either two or three Hits (even for three-plane detectors), if the space point algorithm is configured to allow it.

Definition at line 548 of file SpacePointAlg.cxx.

References art::PtrVector< T >::clear(), compatible(), fillSpacePoint(), trkf::KHitWireX::getHit(), trkf::KHitTrack::getHit(), hits(), numHitMap(), art::PtrVector< T >::push_back(), art::PtrVector< T >::size(), and track.

Referenced by trkf::Track3DKalmanHit::createOutputs(), and trkf::TrackKalmanCheater::produce().

  {
    // Loop over KHitTracks.

    art::PtrVector<recob::Hit> hits;
    art::PtrVector<recob::Hit> compatible_hits;
    for (std::multimap<double, KHitTrack>::const_iterator it = trackMap.begin();
         it != trackMap.end();
         ++it) {
      const KHitTrack& track = (*it).second;

      // Extrack Hit from track.

      const std::shared_ptr<const KHitBase>& hit = track.getHit();
      const KHitWireX* phit = dynamic_cast<const KHitWireX*>(&*hit);
      if (phit != 0) {
        const art::Ptr<recob::Hit> prhit = phit->getHit();

        // Test this hit for compatibility.

        hits.push_back(prhit);
        bool ok = this->compatible(detProp, hits);
        if (!ok) {

          // The new hit is not compatible.  Make a space point out of
          // the last known compatible hits, provided there are at least
          // two.

          if (compatible_hits.size() >= 2) {
            this->fillSpacePoint(detProp, compatible_hits, spts, this->numHitMap());
            compatible_hits.clear();
          }

          // Forget about any previous hits.

          hits.clear();
          hits.push_back(prhit);
        }

        // Update the list of known compatible hits.

        compatible_hits = hits;
      }
    }

    // Maybe make one final space point.

    if (compatible_hits.size() >= 2) {
      this->fillSpacePoint(detProp, compatible_hits, spts, this->numHitMap());
    }
  }

bool trkf::SpacePointAlg::filter ( ) const

inlinenoexcept

Definition at line 86 of file SpacePointAlg.h.

{ return fFilter; }

const art::PtrVector< recob::Hit > & trkf::SpacePointAlg::getAssociatedHits

(

const recob::SpacePoint &

spt

)

const

Definition at line 1441 of file SpacePointAlg.cxx.

References fSptHitMap, and recob::SpacePoint::ID().

Referenced by trkf::SpacePointAna::analyze(), trkf::Track3DKalmanHit::createOutputs(), trkf::SeedFinderAlgorithm::FindSeedAtEnd(), trkf::SeedFinderAlgorithm::FindSeeds(), makeSpacePoints(), trkf::TCTrack::produce(), trkf::SpacePointCheater::produce(), trkf::SpacePointFinder::produce(), and trkf::TrackKalmanCheater::produce().

  {
    // It is an error if no hits are associated with this space point (throw exception).

    std::map<int, art::PtrVector<recob::Hit>>::const_iterator it = fSptHitMap.find(spt.ID());
    if (it == fSptHitMap.end()) {
      mf::LogWarning("SpacePointAlg")
        << "Looking for ID " << spt.ID() << " from " << fSptHitMap.size() << std::endl;
      throw cet::exception("SpacePointAlg") << "No Hits associated with space point.\n";
    }
    return (*it).second;
  }

void trkf::SpacePointAlg::makeMCTruthSpacePoints	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	spts
	)		const

Definition at line 777 of file SpacePointAlg.cxx.

References makeSpacePoints().

Referenced by trkf::SpacePointAna::analyze(), and trkf::SpacePointCheater::produce().

  {
    makeSpacePoints(clockData, detProp, hits, spts, true);
  }

void trkf::SpacePointAlg::makeSpacePoints	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	spts
	)		const

Definition at line 765 of file SpacePointAlg.cxx.

Referenced by trkf::SpacePointAna::analyze(), trkf::SeedFinderAlgorithm::FindSeeds(), trkf::FeatureTracker::Get3DFeaturePoints(), makeMCTruthSpacePoints(), cluster::ClusterMatchAlg::Match_SpacePoint(), trkf::TCTrack::produce(), and trkf::SpacePointFinder::produce().

  {
    makeSpacePoints(clockData, detProp, hits, spts, false);
  }

void trkf::SpacePointAlg::makeSpacePoints ( detinfo::DetectorClocksData const &  clockData, detinfo::DetectorPropertiesData const &  detProp, const art::PtrVector< recob::Hit > &  hits, std::vector< recob::SpacePoint > &  spts, bool  useMC  ) const

private

Todo:

Why are we still checking on whether this is MC or not?

Todo:

Such checks should not be in reconstruction code.

Definition at line 789 of file SpacePointAlg.cxx.

References util::abs(), util::begin(), art::PtrVector< T >::begin(), c1, c2, recob::SpacePoint::Chisq(), art::PtrVector< T >::clear(), compatible(), geo::CryostatID::Cryostat, geo::GeometryCore::Cryostat(), tca::debug, trkf::SpacePointAlg::HitMCInfo::dist2, util::empty(), art::PtrVector< T >::end(), art::PtrVector< T >::erase(), fEnableU, fEnableV, fEnableW, fFilter, fHitMCMap, fillComplexSpacePoint(), fillSpacePoint(), fMaxDT, fMaxS, fMerge, fMinViews, fPreferColl, fSptHitMap, Get, getAssociatedHits(), geo::WireGeo::GetEnd(), geo::WireGeo::GetStart(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), cheat::BackTrackerService::HitToAvgSimIDEs(), mf::isDebugEnabled(), geo::kCollection, geo::kU, geo::kV, geo::kZ, geo::GeometryCore::Ncryostats(), geo::CryostatGeo::NTPC(), trkf::SpacePointAlg::HitMCInfo::pchit, recob::Hit::PeakTime(), geo::PlaneID::Plane, art::PtrVector< T >::push_back(), art::PtrVector< T >::reserve(), cheat::BackTrackerService::SimIDEsToXYZ(), art::PtrVector< T >::size(), util::size(), t1, t2, geo::TPCID::TPC, trkf::SpacePointAlg::HitMCInfo::trackIDs, update(), recob::Hit::View(), geo::WireID::Wire, recob::Hit::WireID(), x, and trkf::SpacePointAlg::HitMCInfo::xyz.

  {
    art::ServiceHandle<geo::Geometry const> geom;

    fSptHitMap.clear();

    // Print diagnostic information.

    update(detProp);

    // First make sure result vector is empty.

    spts.erase(spts.begin(), spts.end());

    // Statistics.

    int n2 = 0;     // Number of two-hit space points.
    int n3 = 0;     // Number of three-hit space points.
    int n2filt = 0; // Number of two-hit space points after filtering/merging.
    int n3filt = 0; // Number of three-hit space pointe after filtering/merging.

    // Sort hits into maps indexed by [cryostat][tpc][plane][wire].
    // If using mc information, also generate maps of sim::IDEs and mc
    // position indexed by hit.

    std::vector<std::vector<std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>>> hitmap;
    fHitMCMap.clear();

    unsigned int ncstat = geom->Ncryostats();
    hitmap.resize(ncstat);

    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

    for (auto const& cryoid : geom->Iterate<geo::CryostatID>()) {
      unsigned int cstat = cryoid.Cryostat;
      unsigned int ntpc = geom->Cryostat(cryoid).NTPC();
      hitmap[cstat].resize(ntpc);
      for (unsigned int tpc = 0; tpc < ntpc; ++tpc) {
        unsigned int nplane = wireReadoutGeom.Nplanes({cstat, tpc});
        hitmap[cstat][tpc].resize(nplane);
      }
    }

    for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
         ++ihit) {
      const art::Ptr<recob::Hit>& phit = *ihit;
      geo::View_t view = phit->View();
      if ((view == geo::kU && fEnableU) || (view == geo::kV && fEnableV) ||
          (view == geo::kZ && fEnableW)) {
        geo::WireID phitWireID = phit->WireID();
        hitmap[phitWireID.Cryostat][phitWireID.TPC][phitWireID.Plane].insert(
          std::make_pair(phitWireID.Wire, phit));
      }
    }

    // Fill mc information, including IDEs and closest neighbors
    // of each hit.
    if (useMC) {
      art::ServiceHandle<cheat::BackTrackerService const> bt_serv;

      // First loop over hits and fill track ids and mc position.
      for (auto const& id : wireReadoutGeom.Iterate<geo::PlaneID>()) {
        auto const [cstat, tpc, plane] = std::make_tuple(id.Cryostat, id.TPC, id.Plane);
        int nplane = wireReadoutGeom.Nplanes(id);
        for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit =
               hitmap[cstat][tpc][plane].begin();
             ihit != hitmap[cstat][tpc][plane].end();
             ++ihit) {
          const art::Ptr<recob::Hit>& phit = ihit->second;
          const recob::Hit& hit = *phit;
          HitMCInfo& mcinfo = fHitMCMap[&hit]; // Default HitMCInfo.

          // Fill default nearest neighbor information (i.e. none).

          mcinfo.pchit.resize(nplane, 0);
          mcinfo.dist2.resize(nplane, 1.e20);

          // Get sim::IDEs for this hit.

          std::vector<sim::IDE> ides = bt_serv->HitToAvgSimIDEs(clockData, phit);

          // Get sorted track ids. for this hit.

          mcinfo.trackIDs.reserve(ides.size());
          for (std::vector<sim::IDE>::const_iterator i = ides.begin(); i != ides.end(); ++i)
            mcinfo.trackIDs.push_back(i->trackID);
          sort(mcinfo.trackIDs.begin(), mcinfo.trackIDs.end());

          // Get position of ionization for this hit.

          try {
            mcinfo.xyz = bt_serv->SimIDEsToXYZ(ides);
          }
          catch (cet::exception& x) {
            mcinfo.xyz.clear();
          }
        } // end loop over ihit
      }   // end loop oer planes

      // Loop over hits again and fill nearest neighbor information for real.
      for (auto const& id : wireReadoutGeom.Iterate<geo::PlaneID>()) {
        auto const [cstat, tpc, plane] = std::make_tuple(id.Cryostat, id.TPC, id.Plane);
        int nplane = wireReadoutGeom.Nplanes(id);
        for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit =
               hitmap[cstat][tpc][plane].begin();
             ihit != hitmap[cstat][tpc][plane].end();
             ++ihit) {
          const art::Ptr<recob::Hit>& phit = ihit->second;
          const recob::Hit& hit = *phit;
          HitMCInfo& mcinfo = fHitMCMap[&hit];
          if (mcinfo.xyz.size() != 0) {
            assert(mcinfo.xyz.size() == 3);

            // Fill nearest neighbor information for this hit.

            for (int plane2 = 0; plane2 < nplane; ++plane2) {
              for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator jhit =
                     hitmap[cstat][tpc][plane2].begin();
                   jhit != hitmap[cstat][tpc][plane2].end();
                   ++jhit) {
                const art::Ptr<recob::Hit>& phit2 = jhit->second;
                const recob::Hit& hit2 = *phit2;
                const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];

                if (mcinfo2.xyz.size() != 0) {
                  assert(mcinfo2.xyz.size() == 3);
                  double dx = mcinfo.xyz[0] - mcinfo2.xyz[0];
                  double dy = mcinfo.xyz[1] - mcinfo2.xyz[1];
                  double dz = mcinfo.xyz[2] - mcinfo2.xyz[2];
                  double dist2 = dx * dx + dy * dy + dz * dz;
                  if (dist2 < mcinfo.dist2[plane2]) {
                    mcinfo.dist2[plane2] = dist2;
                    mcinfo.pchit[plane2] = &hit2;
                  }
                } // end if mcinfo2.xyz valid
              }   // end loop over jhit
            }     // end loop over plane2
          }       // end if mcinfo.xyz valid.
        }         // end loop over ihit
      }           // end loop over planes
    }             // end if MC

    // use mf::LogDebug instead of MF_LOG_DEBUG because we reuse it in many lines
    // insertions are protected by mf::isDebugEnabled()
    mf::LogDebug debug("SpacePointAlg");
    if (mf::isDebugEnabled()) {
      debug << "Total hits = " << hits.size() << "\n\n";
      for (auto const& id : geom->Iterate<geo::TPCID>()) {
        auto const [cstat, tpc] = std::make_tuple(id.Cryostat, id.TPC);
        int nplane = hitmap[cstat][tpc].size();
        for (int plane = 0; plane < nplane; ++plane) {
          debug << "TPC, Plane: " << tpc << ", " << plane
                << ", hits = " << hitmap[cstat][tpc][plane].size() << "\n";
        }
      } // end loop over TPCs
    }   // if debug

    // Make empty multimap from hit pointer on preferred
    // (most-populated or collection) plane to space points that
    // include that hit (used for sorting, filtering, and
    // merging).

    typedef const recob::Hit* sptkey_type;
    std::multimap<sptkey_type, recob::SpacePoint> sptmap;
    std::set<sptkey_type> sptkeys; // Keys of multimap.

    // Loop over TPCs.
    for (auto const& tpcid : geom->Iterate<geo::TPCID>()) {

      auto const [cstat, tpc] = std::make_tuple(tpcid.Cryostat, tpcid.TPC);
      // Sort maps in increasing order of number of hits.
      // This is so that we can do the outer loops over hits
      // over the views with fewer hits.
      //
      // If config parameter PreferColl is true, treat the colleciton
      // plane as if it had the most hits, regardless of how many
      // hits it actually has.  This will force space points to be
      // filtered and merged with respect to the collection plane
      // wires.  It will also force space points to be sorted by
      // collection plane wire.

      int nplane = hitmap[cstat][tpc].size();
      std::vector<int> index(nplane);

      for (int i = 0; i < nplane; ++i)
        index[i] = i;

      for (int i = 0; i < nplane - 1; ++i) {
        for (int j = i + 1; j < nplane; ++j) {
          bool icoll = fPreferColl && wireReadoutGeom.SignalType(geo::PlaneID(tpcid, index[i])) ==
                                        geo::kCollection;
          bool jcoll = fPreferColl && wireReadoutGeom.SignalType(geo::PlaneID(tpcid, index[j])) ==
                                        geo::kCollection;
          if ((hitmap[cstat][tpc][index[i]].size() > hitmap[cstat][tpc][index[j]].size() &&
               !jcoll) ||
              icoll) {
            int temp = index[i];
            index[i] = index[j];
            index[j] = temp;
          }
        }
      } // end loop over i

      // how many views with hits?
      // This will allow for the special case where we might have only 2 planes of information and
      // still want space points even if a three plane TPC
      std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>& hitsByPlaneVec =
        hitmap[cstat][tpc];
      int nViewsWithHits(0);

      for (int i = 0; i < nplane; i++) {
        if (hitsByPlaneVec[index[i]].size() > 0) nViewsWithHits++;
      }

      // If two-view space points are allowed, make a double loop
      // over hits and produce space points for compatible hit-pairs.

      if ((nViewsWithHits == 2 || nplane == 2) && fMinViews <= 2) {

        // Loop over pairs of views.
        for (int i = 0; i < nplane - 1; ++i) {
          unsigned int plane1 = index[i];

          if (hitmap[cstat][tpc][plane1].empty()) continue;

          for (int j = i + 1; j < nplane; ++j) {
            unsigned int plane2 = index[j];

            if (hitmap[cstat][tpc][plane2].empty()) continue;

            // Get angle, pitch, and offset of plane2 wires.
            geo::PlaneID const plane2_id{tpcid, plane2};
            const geo::WireGeo& wgeo2 = wireReadoutGeom.Plane(plane2_id).Wire(0);
            double const hl2 = wgeo2.HalfL();
            auto const xyz21 = wgeo2.GetStart();
            auto const xyz22 = wgeo2.GetEnd();
            double s2 = (xyz22.Y() - xyz21.Y()) / (2. * hl2);
            double c2 = (xyz22.Z() - xyz21.Z()) / (2. * hl2);
            double dist2 = -xyz21.Y() * c2 + xyz21.Z() * s2;
            double pitch2 = wireReadoutGeom.Plane(plane2_id).WirePitch();

            if (!fPreferColl &&
                hitmap[cstat][tpc][plane1].size() > hitmap[cstat][tpc][plane2].size())
              throw cet::exception("SpacePointAlg")
                << "makeSpacePoints(): hitmaps with incompatible size\n";

            // Loop over pairs of hits.

            art::PtrVector<recob::Hit> hitvec;
            hitvec.reserve(2);

            for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit1 =
                   hitmap[cstat][tpc][plane1].begin();
                 ihit1 != hitmap[cstat][tpc][plane1].end();
                 ++ihit1) {

              const art::Ptr<recob::Hit>& phit1 = ihit1->second;
              geo::WireID phit1WireID = phit1->WireID();
              const geo::WireGeo& wgeo = wireReadoutGeom.Wire(phit1WireID);

              // Get endpoint coordinates of this wire.
              // (kept as assertions for performance reasons)
              assert(phit1WireID.Cryostat == cstat);
              assert(phit1WireID.TPC == tpc);
              assert(phit1WireID.Plane == plane1);
              auto const xyz1 = wgeo.GetStart();
              auto const xyz2 = wgeo.GetEnd();

              // Find the plane2 wire numbers corresponding to the endpoints.

              double wire21 = (-xyz1.Y() * c2 + xyz1.Z() * s2 - dist2) / pitch2;
              double wire22 = (-xyz2.Y() * c2 + xyz2.Z() * s2 - dist2) / pitch2;

              int wmin = std::max(0., std::min(wire21, wire22));
              int wmax = std::max(0., std::max(wire21, wire22) + 1.);

              std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
                ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
                ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);

              for (; ihit2 != ihit2end; ++ihit2) {

                const art::Ptr<recob::Hit>& phit2 = ihit2->second;

                // Check current pair of hits for compatibility.
                // By construction, hits should always have compatible views
                // and times, but may not have compatible mc information.

                hitvec.clear();
                hitvec.push_back(phit1);
                hitvec.push_back(phit2);
                bool ok = compatible(detProp, hitvec, useMC);
                if (ok) {

                  // Add a space point.

                  ++n2;

                  // make a dummy vector of recob::SpacePoints
                  // as we are filtering or merging and don't want to
                  // add the created SpacePoint to the final collection just yet
                  // This dummy vector will hold just one recob::SpacePoint,
                  // which will go into the multimap and then the vector
                  // will go out of scope.

                  std::vector<recob::SpacePoint> sptv;
                  fillSpacePoint(detProp, hitvec, sptv, sptmap.size());
                  sptkey_type key = &*phit2;
                  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
                  sptkeys.insert(key);
                }
              }
            }
          }
        }
      } // end if fMinViews <= 2

      // If three-view space points are allowed, make a triple loop
      // over hits and produce space points for compatible triplets.

      if (nplane >= 3 && fMinViews <= 3) {

        // Loop over triplets of hits.

        art::PtrVector<recob::Hit> hitvec;
        hitvec.reserve(3);

        unsigned int plane1 = index[0];
        unsigned int plane2 = index[1];
        unsigned int plane3 = index[2];

        // Get angle, pitch, and offset of plane1 wires.

        geo::PlaneID const plane1_id{tpcid, plane1};
        const geo::WireGeo& wgeo1 = wireReadoutGeom.Plane(plane1_id).Wire(0);
        double const hl1 = wgeo1.HalfL();
        auto const xyz11 = wgeo1.GetStart();
        auto const xyz12 = wgeo1.GetEnd();
        double s1 = (xyz12.Y() - xyz11.Y()) / (2. * hl1);
        double c1 = (xyz12.Z() - xyz11.Z()) / (2. * hl1);
        double dist1 = -xyz11.Y() * c1 + xyz11.Z() * s1;
        double pitch1 = wireReadoutGeom.Plane(plane1_id).WirePitch();
        const double TicksOffset1 = detProp.GetXTicksOffset(plane1_id);

        // Get angle, pitch, and offset of plane2 wires.

        geo::PlaneID const plane2_id{tpcid, plane2};
        const geo::WireGeo& wgeo2 = wireReadoutGeom.Plane(plane2_id).Wire(0);
        double const hl2 = wgeo2.HalfL();
        auto const xyz21 = wgeo2.GetStart();
        auto const xyz22 = wgeo2.GetEnd();
        double s2 = (xyz22.Y() - xyz21.Y()) / (2. * hl2);
        double c2 = (xyz22.Z() - xyz21.Z()) / (2. * hl2);
        double dist2 = -xyz21.Y() * c2 + xyz21.Z() * s2;
        double pitch2 = wireReadoutGeom.Plane(plane2_id).WirePitch();
        const double TicksOffset2 = detProp.GetXTicksOffset(plane2_id);

        // Get angle, pitch, and offset of plane3 wires.

        geo::PlaneID const plane3_id{tpcid, plane3};
        const geo::WireGeo& wgeo3 = wireReadoutGeom.Plane(plane3_id).Wire(0);
        double const hl3 = wgeo3.HalfL();
        auto const xyz31 = wgeo3.GetStart();
        auto const xyz32 = wgeo3.GetEnd();
        double s3 = (xyz32.Y() - xyz31.Y()) / (2. * hl3);
        double c3 = (xyz32.Z() - xyz31.Z()) / (2. * hl3);
        double dist3 = -xyz31.Y() * c3 + xyz31.Z() * s3;
        double pitch3 = wireReadoutGeom.Plane(plane3_id).WirePitch();
        const double TicksOffset3 = detProp.GetXTicksOffset(plane3_id);

        // Get sine of angle differences.

        double s12 = s1 * c2 - s2 * c1; // sin(theta1 - theta2).
        double s23 = s2 * c3 - s3 * c2; // sin(theta2 - theta3).
        double s31 = s3 * c1 - s1 * c3; // sin(theta3 - theta1).

        // Loop over hits in plane1.

        std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
          ihit1 = hitmap[cstat][tpc][plane1].begin(),
          ihit1end = hitmap[cstat][tpc][plane1].end();
        for (; ihit1 != ihit1end; ++ihit1) {

          unsigned int wire1 = ihit1->first;
          const art::Ptr<recob::Hit>& phit1 = ihit1->second;
          geo::WireID phit1WireID = phit1->WireID();
          const geo::WireGeo& wgeo = wireReadoutGeom.Wire(phit1WireID);

          // Get endpoint coordinates of this wire from plane1.
          // (kept as assertions for performance reasons)
          assert(phit1WireID.Cryostat == cstat);
          assert(phit1WireID.TPC == tpc);
          assert(phit1WireID.Plane == plane1);
          assert(phit1WireID.Wire == wire1);
          auto const xyz1 = wgeo.GetStart();
          auto const xyz2 = wgeo.GetEnd();

          // Get corrected time and oblique coordinate of first hit.

          double t1 = phit1->PeakTime() - TicksOffset1;
          double u1 = wire1 * pitch1 + dist1;

          // Find the plane2 wire numbers corresponding to the endpoints.

          double wire21 = (-xyz1.Y() * c2 + xyz1.Z() * s2 - dist2) / pitch2;
          double wire22 = (-xyz2.Y() * c2 + xyz2.Z() * s2 - dist2) / pitch2;

          int wmin = std::max(0., std::min(wire21, wire22));
          int wmax = std::max(0., std::max(wire21, wire22) + 1.);

          std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
            ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
            ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);

          for (; ihit2 != ihit2end; ++ihit2) {

            int wire2 = ihit2->first;
            const art::Ptr<recob::Hit>& phit2 = ihit2->second;

            // Get corrected time of second hit.

            double t2 = phit2->PeakTime() - TicksOffset2;

            // Check maximum time difference with first hit.

            bool dt12ok = std::abs(t1 - t2) <= fMaxDT;
            if (dt12ok) {

              // Test first two hits for compatibility before looping
              // over third hit.

              hitvec.clear();
              hitvec.push_back(phit1);
              hitvec.push_back(phit2);
              bool h12ok = compatible(detProp, hitvec, useMC);
              if (h12ok) {

                // Get oblique coordinate of second hit.

                double u2 = wire2 * pitch2 + dist2;

                // Predict plane3 oblique coordinate and wire number.

                double u3pred = (-u1 * s23 - u2 * s31) / s12;
                double w3pred = (u3pred - dist3) / pitch3;
                double w3delta = std::abs(fMaxS / (s12 * pitch3));
                int w3min = std::max(0., std::ceil(w3pred - w3delta));
                int w3max = std::max(0., std::floor(w3pred + w3delta));

                std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
                  ihit3 = hitmap[cstat][tpc][plane3].lower_bound(w3min),
                  ihit3end = hitmap[cstat][tpc][plane3].upper_bound(w3max);

                for (; ihit3 != ihit3end; ++ihit3) {

                  int wire3 = ihit3->first;
                  const art::Ptr<recob::Hit>& phit3 = ihit3->second;

                  // Get corrected time of third hit.

                  double t3 = phit3->PeakTime() - TicksOffset3;

                  // Check time difference of third hit compared to first two hits.

                  bool dt123ok = std::abs(t1 - t3) <= fMaxDT && std::abs(t2 - t3) <= fMaxDT;
                  if (dt123ok) {

                    // Get oblique coordinate of third hit and check spatial separation.

                    double u3 = wire3 * pitch3 + dist3;
                    double S = s23 * u1 + s31 * u2 + s12 * u3;
                    bool sok = std::abs(S) <= fMaxS;
                    if (sok) {

                      // Test triplet for compatibility.

                      hitvec.clear();
                      hitvec.push_back(phit1);
                      hitvec.push_back(phit2);
                      hitvec.push_back(phit3);
                      bool h123ok = compatible(detProp, hitvec, useMC);
                      if (h123ok) {

                        // Add a space point.

                        ++n3;

                        // make a dummy vector of recob::SpacePoints
                        // as we are filtering or merging and don't want to
                        // add the created SpacePoint to the final collection just yet
                        // This dummy vector will hold just one recob::SpacePoint,
                        // which will go into the multimap and then the vector
                        // will go out of scope.

                        std::vector<recob::SpacePoint> sptv;
                        fillSpacePoint(detProp, hitvec, sptv, sptmap.size() - 1);
                        sptkey_type key = &*phit3;
                        sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
                        sptkeys.insert(key);
                      }
                    }
                  }
                }
              }
            }
          }
        }
      } // end if fMinViews <= 3

      // Do Filtering.

      if (fFilter) {

        // Transfer (some) space points from sptmap to spts.

        spts.reserve(spts.size() + sptkeys.size());

        // Loop over keys of space point map.
        // Space points that have the same key are candidates for filtering.

        for (std::set<sptkey_type>::const_iterator i = sptkeys.begin(); i != sptkeys.end(); ++i) {
          sptkey_type key = *i;

          // Loop over space points corresponding to the current key.
          // Choose the single best space point from among this group.

          double best_chisq = 0.;
          const recob::SpacePoint* best_spt = 0;

          for (std::multimap<sptkey_type, recob::SpacePoint>::const_iterator j =
                 sptmap.lower_bound(key);
               j != sptmap.upper_bound(key);
               ++j) {
            const recob::SpacePoint& spt = j->second;
            if (best_spt == 0 || spt.Chisq() < best_chisq) {
              best_spt = &spt;
              best_chisq = spt.Chisq();
            }
          }

          // Transfer best filtered space point to result vector.

          if (!best_spt)
            throw cet::exception("SpacePointAlg") << "makeSpacePoints(): no best point\n";
          spts.push_back(*best_spt);
          if (fMinViews <= 2)
            ++n2filt;
          else
            ++n3filt;
        }
      } // end if filtering

      // Do merging.

      else if (fMerge) {

        // Transfer merged space points from sptmap to spts.

        spts.reserve(spts.size() + sptkeys.size());

        // Loop over keys of space point map.
        // Space points that have the same key are candidates for merging.

        for (std::set<sptkey_type>::const_iterator i = sptkeys.begin(); i != sptkeys.end(); ++i) {
          sptkey_type key = *i;

          // Loop over space points corresponding to the current key.
          // Make a collection of hits that is the union of the hits
          // from each candidate space point.

          std::multimap<sptkey_type, recob::SpacePoint>::const_iterator jSPT =
                                                                          sptmap.lower_bound(key),
                                                                        jSPTend =
                                                                          sptmap.upper_bound(key);

          art::PtrVector<recob::Hit> merged_hits;
          for (; jSPT != jSPTend; ++jSPT) {
            const recob::SpacePoint& spt = jSPT->second;

            // Loop over hits from this space points.
            // Add each hit to the collection of all hits.

            const art::PtrVector<recob::Hit>& spt_hits = getAssociatedHits(spt);
            merged_hits.reserve(merged_hits.size() +
                                spt_hits.size()); // better than nothing, but not ideal
            for (art::PtrVector<recob::Hit>::const_iterator k = spt_hits.begin();
                 k != spt_hits.end();
                 ++k) {
              const art::Ptr<recob::Hit>& hit = *k;
              merged_hits.push_back(hit);
            }
          }

          // Remove duplicates.

          std::sort(merged_hits.begin(), merged_hits.end());
          art::PtrVector<recob::Hit>::iterator it =
            std::unique(merged_hits.begin(), merged_hits.end());
          merged_hits.erase(it, merged_hits.end());

          // Construct a complex space points using merged hits.

          fillComplexSpacePoint(detProp, merged_hits, spts, sptmap.size() + spts.size() - 1);

          if (fMinViews <= 2)
            ++n2filt;
          else
            ++n3filt;
        }
      } // end if merging

      // No filter, no merge.

      else {

        // Transfer all space points from sptmap to spts.

        spts.reserve(spts.size() + sptkeys.size());

        // Loop over space points.

        for (std::multimap<sptkey_type, recob::SpacePoint>::const_iterator j = sptmap.begin();
             j != sptmap.end();
             ++j) {
          const recob::SpacePoint& spt = j->second;
          spts.push_back(spt);
        }

        // Update statistics.

        n2filt = n2;
        n3filt = n3;
      }
    } // end loop over tpcs

    if (mf::isDebugEnabled()) {
      debug << "\n2-hit space points = " << n2 << "\n"
            << "3-hit space points = " << n3 << "\n"
            << "2-hit filtered/merged space points = " << n2filt << "\n"
            << "3-hit filtered/merged space points = " << n3filt;
    } // if debug
  }

double trkf::SpacePointAlg::maxDT ( ) const

inlinenoexcept

Definition at line 88 of file SpacePointAlg.h.

Referenced by cluster::ClusterMatchAlg::Match_SpacePoint().

{ return fMaxDT; }

double trkf::SpacePointAlg::maxS ( ) const

inlinenoexcept

Definition at line 89 of file SpacePointAlg.h.

{ return fMaxS; }

bool trkf::SpacePointAlg::merge ( ) const

inlinenoexcept

Definition at line 87 of file SpacePointAlg.h.

Referenced by trkf::SpacePointAna::analyze(), and trkf::SpacePointAna::bookHistograms().

{ return fMerge; }

int trkf::SpacePointAlg::minViews ( ) const

inlinenoexcept

Definition at line 90 of file SpacePointAlg.h.

Referenced by trkf::SpacePointCheater::produce(), and trkf::SpacePointFinder::produce().

{ return fMinViews; }

int trkf::SpacePointAlg::numHitMap ( ) const

inline

Definition at line 152 of file SpacePointAlg.h.

Referenced by fillSpacePoints().

{ return fSptHitMap.size(); }

double trkf::SpacePointAlg::separation

(

const art::PtrVector< recob::Hit > &

hits

)

const

Definition at line 179 of file SpacePointAlg.cxx.

References geo::CryostatID::Cryostat, larg4::dist(), Get, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), geo::WireGeo::HalfL(), geo::PlaneID::Plane, art::PtrVector< T >::size(), geo::TPCID::TPC, and recob::Hit::WireID().

Referenced by trkf::SpacePointAna::analyze().

  {
    // Trivial case - fewer than three hits.

    if (hits.size() < 3) return 0.;

    // Error case - more than three hits.

    if (hits.size() > 3) {
      mf::LogError("SpacePointAlg") << "Method separation called with more than three htis.";
      return 0.;
    }

    // Got exactly three hits.

    // Calculate angles and distance of each hit from origin.

    double dist[3] = {0., 0., 0.};
    double sinth[3] = {0., 0., 0.};
    double costh[3] = {0., 0., 0.};
    unsigned int cstats[3];
    unsigned int tpcs[3];
    unsigned int planes[3];

    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

    for (int i = 0; i < 3; ++i) {

      // Get tpc, plane, wire.

      const recob::Hit& hit = *(hits[i]);
      const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());
      cstats[i] = hit.WireID().Cryostat;
      tpcs[i] = hit.WireID().TPC;
      planes[i] = hit.WireID().Plane;

      // Check tpc and plane errors.

      for (int j = 0; j < i; ++j) {

        if (cstats[j] != hit.WireID().Cryostat) {
          mf::LogError("SpacePointAlg")
            << "Method separation called with hits from multiple cryostats..";
          return 0.;
        }

        if (tpcs[j] != hit.WireID().TPC) {
          mf::LogError("SpacePointAlg")
            << "Method separation called with hits from multiple tpcs..";
          return 0.;
        }

        if (planes[j] == hit.WireID().Plane) {
          mf::LogError("SpacePointAlg")
            << "Method separation called with hits from the same plane..";
          return 0.;
        }
      }

      // Get angles and distance of wire.

      double const hl = wgeom.HalfL();
      auto const xyz = wgeom.GetCenter();
      auto const xyz1 = wgeom.GetEnd();
      double s = (xyz1.Y() - xyz.Y()) / hl;
      double c = (xyz1.Z() - xyz.Z()) / hl;
      sinth[hit.WireID().Plane] = s;
      costh[hit.WireID().Plane] = c;
      dist[hit.WireID().Plane] = xyz.Z() * s - xyz.Y() * c;
    }

    double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0] +
                (sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1] +
                (sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
    return S;
  }

void trkf::SpacePointAlg::update

(

detinfo::DetectorPropertiesData const &

detProp

)

const

Definition at line 86 of file SpacePointAlg.cxx.

References Get, detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::kCollection, geo::kHorizontal, geo::kInduction, geo::kU, geo::kV, geo::kVertical, and geo::kZ.

Referenced by makeSpacePoints().

  {
    // Generate info report on first call only.

    static bool first = true;
    bool report = first;
    first = false;

    // Get services.

    art::ServiceHandle<geo::Geometry const> geom;
    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

    // Calculate and print geometry information.

    if (report) mf::LogInfo("SpacePointAlg") << "Updating geometry constants.\n";

    for (auto const& plane : wireReadoutGeom.Iterate<geo::PlaneGeo>()) {

      // Fill view-dependent quantities.

      geo::View_t view = plane.View();
      std::string viewname = "?";
      if (view == geo::kU) { viewname = "U"; }
      else if (view == geo::kV) {
        viewname = "V";
      }
      else if (view == geo::kZ) {
        viewname = "Z";
      }
      else
        throw cet::exception("SpacePointAlg") << "Bad view = " << view << "\n";

      std::string sigtypename = "?";
      geo::SigType_t sigtype = wireReadoutGeom.SignalType(plane.ID());
      if (sigtype == geo::kInduction)
        sigtypename = "Induction";
      else if (sigtype == geo::kCollection)
        sigtypename = "Collection";
      else
        throw cet::exception("SpacePointAlg") << "Bad signal type = " << sigtype << "\n";

      std::string orientname = "?";
      geo::Orient_t orient = plane.Orientation();
      if (orient == geo::kVertical)
        orientname = "Vertical";
      else if (orient == geo::kHorizontal)
        orientname = "Horizontal";
      else
        throw cet::exception("SpacePointAlg") << "Bad orientation = " << orient << "\n";

      if (report) {
        auto const xyz = plane.GetCenter();
        auto const& planeid = plane.ID();
        mf::LogInfo("SpacePointAlg")
          << '\n'
          << planeid << '\n'
          << "  View: " << viewname << "\n"
          << "  SignalType: " << sigtypename << "\n"
          << "  Orientation: " << orientname << "\n"
          << "  Plane location: " << xyz.X() << "\n"
          << "  Plane pitch: " << wireReadoutGeom.Plane0Pitch(planeid.parentID(), planeid.Plane)
          << "\n"
          << "  Wire angle: " << plane.Wire(0).ThetaZ() << "\n"
          << "  Wire pitch: " << wireReadoutGeom.Plane(planeid).WirePitch() << "\n"
          << "  Time offset: " << detProp.GetXTicksOffset(planeid) << "\n";
      }

      if (orient != geo::kVertical)
        throw cet::exception("SpacePointAlg") << "Horizontal wire geometry not implemented.\n";
    } // end loop over planes
  }

Member Data Documentation

bool trkf::SpacePointAlg::fEnableU

private

Enable flag (U).

Definition at line 168 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fEnableV

private

Enable flag (V).

Definition at line 169 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fEnableW

private

Enable flag (W).

Definition at line 170 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fFilter

private

Filter flag.

Definition at line 171 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

std::map<const recob::Hit*, HitMCInfo> trkf::SpacePointAlg::fHitMCMap

mutableprivate

Definition at line 186 of file SpacePointAlg.h.

Referenced by compatible(), and makeSpacePoints().

double trkf::SpacePointAlg::fMaxDT

private

Maximum time difference between planes.

Definition at line 165 of file SpacePointAlg.h.

Referenced by compatible(), makeSpacePoints(), and SpacePointAlg().

double trkf::SpacePointAlg::fMaxS

private

Maximum space separation between wires.

Definition at line 166 of file SpacePointAlg.h.

Referenced by compatible(), makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fMerge

private

Merge flag.

Definition at line 172 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

int trkf::SpacePointAlg::fMinViews

private

Mininum number of views per space point.

Definition at line 167 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fPreferColl

private

Sort by collection wire.

Definition at line 173 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

std::map<int, art::PtrVector<recob::Hit> > trkf::SpacePointAlg::fSptHitMap

mutableprivate

Definition at line 187 of file SpacePointAlg.h.

Referenced by fillComplexSpacePoint(), fillSpacePoint(), getAssociatedHits(), and makeSpacePoints().

double trkf::SpacePointAlg::fTickOffsetU

private

Tick offset for plane U.

Definition at line 174 of file SpacePointAlg.h.

Referenced by correctedTime(), and SpacePointAlg().

double trkf::SpacePointAlg::fTickOffsetV

private

Tick offset for plane V.

Definition at line 175 of file SpacePointAlg.h.

Referenced by correctedTime(), and SpacePointAlg().

double trkf::SpacePointAlg::fTickOffsetW

private

Tick offset for plane W.

Definition at line 176 of file SpacePointAlg.h.

Referenced by correctedTime(), and SpacePointAlg().

---

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

• larreco/v10_01_07/source/larreco/RecoAlg/SpacePointAlg.h
• larreco/v10_01_07/source/larreco/RecoAlg/SpacePointAlg.cxx