#include "TripletFinder.h"

Public Member Functions
	TripletFinder (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &xhits, const std::vector< art::Ptr< recob::Hit >> &uhits, const std::vector< art::Ptr< recob::Hit >> &vhits, const std::vector< raw::ChannelID_t > &xbad, const std::vector< raw::ChannelID_t > &ubad, const std::vector< raw::ChannelID_t > &vbad, double distThresh, double distThreshDrift, double xhitOffset, int maxTriplets=0)

std::vector< HitTriplet >	Triplets ()

std::vector< HitTriplet >	TripletsTwoView ()
	Only search for XU intersections. More...

Protected Member Functions
void	FillHitMap (const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits, std::map< geo::TPCID, std::vector< HitOrChan >> &out)
	Helper for constructor. More...

void	FillBadMap (const std::vector< raw::ChannelID_t > &bads, std::map< geo::TPCID, std::vector< raw::ChannelID_t >> &out)
	Helper for constructor. More...

bool	CloseDrift (double xa, double xb) const

bool	CloseSpace (geo::WireIDIntersection ra, geo::WireIDIntersection rb) const

std::vector< ChannelDoublet >	DoubletsXU (geo::TPCID tpc)

std::vector< ChannelDoublet >	DoubletsXV (geo::TPCID tpc)

std::vector< ChannelDoublet >	DoubletHelper (geo::TPCID tpc, const std::vector< HitOrChan > &ahits, const std::vector< HitOrChan > &bhits, const std::vector< raw::ChannelID_t > &bbads) const

Protected Attributes
const geo::GeometryCore *	geom

const geo::WireReadoutGeom *	wireReadoutGeom

double	fDistThresh

double	fDistThreshDrift

double	fXHitOffset

int	fMaxTriplets

std::map< geo::TPCID, std::vector< HitOrChan > >	fX_by_tpc

std::map< geo::TPCID, std::vector< HitOrChan > >	fU_by_tpc

std::map< geo::TPCID, std::vector< HitOrChan > >	fV_by_tpc

std::map< geo::TPCID, std::vector< raw::ChannelID_t > >	fXbad_by_tpc

std::map< geo::TPCID, std::vector< raw::ChannelID_t > >	fUbad_by_tpc

std::map< geo::TPCID, std::vector< raw::ChannelID_t > >	fVbad_by_tpc

Detailed Description

Definition at line 49 of file TripletFinder.h.

Constructor & Destructor Documentation

reco3d::TripletFinder::TripletFinder	(	const detinfo::DetectorPropertiesData &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	xhits,
		const std::vector< art::Ptr< recob::Hit >> &	uhits,
		const std::vector< art::Ptr< recob::Hit >> &	vhits,
		const std::vector< raw::ChannelID_t > &	xbad,
		const std::vector< raw::ChannelID_t > &	ubad,
		const std::vector< raw::ChannelID_t > &	vbad,
		double	distThresh,
		double	distThreshDrift,
		double	xhitOffset,
		int	maxTriplets = `0`
	)

Definition at line 12 of file TripletFinder.cxx.

References fDistThresh, fDistThreshDrift, FillBadMap(), FillHitMap(), fMaxTriplets, fU_by_tpc, fUbad_by_tpc, fV_by_tpc, fVbad_by_tpc, fX_by_tpc, fXbad_by_tpc, fXHitOffset, and Get.

     : wireReadoutGeom{&art::ServiceHandle<geo::WireReadout const>()->Get()}
     , fDistThresh(distThresh)
     , fDistThreshDrift(distThreshDrift)
     , fXHitOffset(xhitOffset)
     , fMaxTriplets(maxTriplets)
   {
     FillHitMap(detProp, xhits, fX_by_tpc);
     FillHitMap(detProp, uhits, fU_by_tpc);
     FillHitMap(detProp, vhits, fV_by_tpc);
 
     FillBadMap(xbad, fXbad_by_tpc);
     FillBadMap(ubad, fUbad_by_tpc);
     FillBadMap(vbad, fVbad_by_tpc);
   }

Member Function Documentation

bool reco3d::TripletFinder::CloseDrift	(	double	xa,
		double	xb
	)		const

protected

Definition at line 114 of file TripletFinder.cxx.

References fDistThreshDrift.

Referenced by DoubletHelper(), and Triplets().

   {
     return fabs(xa - xb) < fDistThreshDrift;
   }

bool reco3d::TripletFinder::CloseSpace	(	geo::WireIDIntersection	ra,
		geo::WireIDIntersection	rb
	)		const

protected

Definition at line 120 of file TripletFinder.cxx.

References fDistThresh, geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by Triplets().

   {
     const TVector3 pa(ra.y, ra.z, 0);
     const TVector3 pb(rb.y, rb.z, 0);
 
     return (pa - pb).Mag() < fDistThresh;
   }

std::vector< ChannelDoublet > reco3d::TripletFinder::DoubletHelper	(	geo::TPCID	tpc,
		const std::vector< HitOrChan > &	ahits,
		const std::vector< HitOrChan > &	bhits,
		const std::vector< raw::ChannelID_t > &	bbads
	)		const

protected

Definition at line 289 of file TripletFinder.cxx.

References reco3d::HitOrChan::chan, CloseDrift(), pt, wireReadoutGeom, and reco3d::HitOrChan::xpos.

Referenced by DoubletsXU(), and DoubletsXV().

   {
     std::vector<ChannelDoublet> ret;
 
     IntersectionCache isect(wireReadoutGeom, tpc);
 
     auto b_begin = bhits.begin();
 
     for (const HitOrChan& a : ahits) {
       // Bad channels are easy because there's no timing constraint
       for (raw::ChannelID_t b : bbads) {
         if (auto pt = isect(a.chan, b)) { ret.emplace_back(a, b, *pt); }
       }
 
       while (b_begin != bhits.end() && b_begin->xpos < a.xpos && !CloseDrift(b_begin->xpos, a.xpos))
         ++b_begin;
 
       for (auto bit = b_begin; bit != bhits.end(); ++bit) {
         const HitOrChan& b = *bit;
 
         if (b.xpos > a.xpos && !CloseDrift(b.xpos, a.xpos)) break;
 
         auto pt = isect(a.chan, b.chan);
         if (!pt) continue;
 
         ret.emplace_back(a, b, *pt);
       } // end for b
     }   // end for a
 
     return ret;
   }

std::vector< ChannelDoublet > reco3d::TripletFinder::DoubletsXU ( geo::TPCID tpc )

protected

Definition at line 261 of file TripletFinder.cxx.

References DoubletHelper(), fU_by_tpc, fUbad_by_tpc, fX_by_tpc, and fXbad_by_tpc.

Referenced by Triplets(), and TripletsTwoView().

   {
     std::vector<ChannelDoublet> ret =
       DoubletHelper(tpc, fX_by_tpc[tpc], fU_by_tpc[tpc], fUbad_by_tpc[tpc]);
 
     // Find X(bad)+U(good) doublets, have to flip them for the final result
     for (auto it : DoubletHelper(tpc, fU_by_tpc[tpc], {}, fXbad_by_tpc[tpc])) {
       ret.push_back({it.b, it.a, it.pt});
     }
 
     return ret;
   }

std::vector< ChannelDoublet > reco3d::TripletFinder::DoubletsXV ( geo::TPCID tpc )

protected

Definition at line 275 of file TripletFinder.cxx.

References DoubletHelper(), fV_by_tpc, fVbad_by_tpc, fX_by_tpc, and fXbad_by_tpc.

Referenced by Triplets().

   {
     std::vector<ChannelDoublet> ret =
       DoubletHelper(tpc, fX_by_tpc[tpc], fV_by_tpc[tpc], fVbad_by_tpc[tpc]);
 
     // Find X(bad)+V(good) doublets, have to flip them for the final result
     for (auto it : DoubletHelper(tpc, fV_by_tpc[tpc], {}, fXbad_by_tpc[tpc])) {
       ret.push_back({it.b, it.a, it.pt});
     }
 
     return ret;
   }

void reco3d::TripletFinder::FillBadMap	(	const std::vector< raw::ChannelID_t > &	bads,
		std::map< geo::TPCID, std::vector< raw::ChannelID_t >> &	out
	)

protected

Helper for constructor.

Definition at line 62 of file TripletFinder.cxx.

References geo::WireReadoutGeom::ChannelToROP(), geo::WireReadoutGeom::ROPtoTPCs(), and wireReadoutGeom.

Referenced by TripletFinder().

   {
     for (raw::ChannelID_t chan : bads) {
       for (geo::TPCID tpc : wireReadoutGeom->ROPtoTPCs(wireReadoutGeom->ChannelToROP(chan))) {
         out[tpc].push_back(chan);
       }
     }
   }

void reco3d::TripletFinder::FillHitMap	(	const detinfo::DetectorPropertiesData &	clockData,
		const std::vector< art::Ptr< recob::Hit >> &	hits,
		std::map< geo::TPCID, std::vector< HitOrChan >> &	out
	)

protected

Helper for constructor.

Definition at line 39 of file TripletFinder.cxx.

References geo::WireReadoutGeom::ChannelToROP(), geo::WireReadoutGeom::ChannelToWire(), detinfo::DetectorPropertiesData::ConvertTicksToX(), fXHitOffset, hits(), geo::kCollection, geo::WireReadoutGeom::ROPtoTPCs(), geo::WireReadoutGeom::SignalType(), and wireReadoutGeom.

Referenced by TripletFinder().

   {
     for (const art::Ptr<recob::Hit>& hit : hits) {
       for (geo::TPCID tpc :
            wireReadoutGeom->ROPtoTPCs(wireReadoutGeom->ChannelToROP(hit->Channel()))) {
         double xpos = 0;
         for (geo::WireID wire : wireReadoutGeom->ChannelToWire(hit->Channel())) {
           if (geo::TPCID(wire) == tpc) {
             xpos = detProp.ConvertTicksToX(hit->PeakTime(), wire);
             if (wireReadoutGeom->SignalType(wire) == geo::kCollection) xpos += fXHitOffset;
           }
         }
 
         out[tpc].emplace_back(hit.get(), xpos);
       }
     }
     for (auto& it : out)
       std::sort(it.second.begin(), it.second.end(), [](auto a, auto b) { return a.xpos < b.xpos; });
   }

std::vector< HitTriplet > reco3d::TripletFinder::Triplets ( )

Definition at line 143 of file TripletFinder.cxx.

References reco3d::HitOrChan::chan, CloseDrift(), CloseSpace(), DoubletsXU(), DoubletsXV(), geo::TPCID::first(), fMaxTriplets, fX_by_tpc, reco3d::HitOrChan::hit, reco3d::LessThanXHit(), pt, reco3d::SameXHit(), wireReadoutGeom, x, and reco3d::HitOrChan::xpos.

Referenced by reco3d::SpacePointSolver::produce().

   {
     std::vector<HitTriplet> ret;
 
     for (const auto& it : fX_by_tpc) {
       const geo::TPCID& tpc = it.first;
 
       std::vector<ChannelDoublet> xus = DoubletsXU(tpc);
       std::vector<ChannelDoublet> xvs = DoubletsXV(tpc);
 
       // Cache to prevent repeating the same questions
       IntersectionCache isectUV(wireReadoutGeom, tpc);
 
       // For the efficient looping below to work we need to sort the doublet
       // lists so the X hits occur in the same order.
       std::sort(xus.begin(), xus.end(), LessThanXHit);
       std::sort(xvs.begin(), xvs.end(), LessThanXHit);
 
       auto xvit_begin = xvs.begin();
 
       int nxuv = 0;
       for (const ChannelDoublet& xu : xus) {
         const HitOrChan& x = xu.a;
         const HitOrChan& u = xu.b;
 
         // Catch up until we're looking at the same X hit in XV
         while (xvit_begin != xvs.end() && LessThanXHit(*xvit_begin, xu))
           ++xvit_begin;
 
         // Loop through all those matching hits
         for (auto xvit = xvit_begin; xvit != xvs.end() && SameXHit(*xvit, xu); ++xvit) {
           if (fMaxTriplets > 0 and nxuv > fMaxTriplets) break;
           const HitOrChan& v = xvit->b;
 
           // Only allow one bad channel per triplet
           if (!x.hit && !u.hit) continue;
           if (!x.hit && !v.hit) continue;
           if (!u.hit && !v.hit) continue;
 
           if (u.hit && v.hit && !CloseDrift(u.xpos, v.xpos)) continue;
 
           auto maybe_ptUV = isectUV(u.chan, v.chan);
           if (!maybe_ptUV) continue;
 
           auto const& ptUV = *maybe_ptUV;
           if (!CloseSpace(xu.pt, xvit->pt) || !CloseSpace(xu.pt, ptUV) ||
               !CloseSpace(xvit->pt, ptUV))
             continue;
 
           double xavg = 0;
           int nx = 0;
           if (x.hit) {
             xavg += x.xpos;
             ++nx;
           }
           if (u.hit) {
             xavg += u.xpos;
             ++nx;
           }
           if (v.hit) {
             xavg += v.xpos;
             ++nx;
           }
           xavg /= nx;
 
           const XYZ pt{
             xavg, (xu.pt.y + xvit->pt.y + ptUV.y) / 3, (xu.pt.z + xvit->pt.z + ptUV.z) / 3};
 
           ret.emplace_back(HitTriplet{x.hit, u.hit, v.hit, pt});
           ++nxuv;
           if (fMaxTriplets > 0 and nxuv > fMaxTriplets) break;
         } // end for xv
       }   // end for xu
 
       std::cout << tpc << " " << xus.size() << " XUs and " << xvs.size() << " XVs -> " << nxuv
                 << " XUVs" << std::endl;
 
     } // end for tpc
 
     std::cout << ret.size() << " XUVs total" << std::endl;
 
     return ret;
   }

std::vector< HitTriplet > reco3d::TripletFinder::TripletsTwoView ( )

Only search for XU intersections.

Definition at line 228 of file TripletFinder.cxx.

References DoubletsXU(), geo::TPCID::first(), fX_by_tpc, reco3d::HitOrChan::hit, pt, x, and reco3d::HitOrChan::xpos.

   {
     std::vector<HitTriplet> ret;
 
     for (const auto& it : fX_by_tpc) {
       const geo::TPCID& tpc = it.first;
 
       std::vector<ChannelDoublet> xus = DoubletsXU(tpc);
 
       for (const ChannelDoublet& xu : xus) {
         const HitOrChan& x = xu.a;
         const HitOrChan& u = xu.b;
 
         double xavg = x.xpos;
         int nx = 1;
         if (u.hit) {
           xavg += u.xpos;
           ++nx;
         }
         xavg /= nx;
 
         const XYZ pt{xavg, xu.pt.y, xu.pt.z};
 
         ret.emplace_back(HitTriplet{x.hit, u.hit, 0, pt});
       } // end for xu
     }   // end for tpc
 
     std::cout << ret.size() << " XUs total" << std::endl;
 
     return ret;
   }