#include "APAGeometryAlg.h"

Public Member Functions
	APAGeometryAlg ()

bool	APAChannelsIntersect (uint32_t chan1, uint32_t chan2, std::vector< geo::WireIDIntersection > &IntersectVector) const
	If the channels intersect, get all intersections. More...

bool	LineSegChanIntersect (geo::Point_t const &xyzStart, geo::Point_t const &xyzEnd, uint32_t chan, std::vector< geo::WireID > &widsCrossed, bool ExtendLine) const
	If a line given by start/end points intersects a channel. More...

std::vector< geo::WireID >	ChanSegsPerSide (uint32_t chan, unsigned int side) const

std::vector< geo::WireID >	ChanSegsPerSide (std::vector< geo::WireID > wids, unsigned int side) const

std::vector< double >	ThreeChanPos (uint32_t u, uint32_t v, uint32_t z) const
	Find the center of the 3 intersections, choose best if multiple. More...

geo::WireID	NearestWireIDOnChan (geo::Point_t const &WorldLoc, uint32_t chan, geo::PlaneID const &planeID) const

unsigned int	ChannelToAPA (uint32_t chan) const
	Get number of the APA containing the given channel. More...

void	ChannelToAPA (uint32_t chan, unsigned int &apa, unsigned int &cryo) const

APAView_t	APAView (uint32_t chan) const
	Get which of the 4 APA views the channel is in. More...

unsigned int	ChannelsInView (geo::View_t geoview) const

uint32_t	FirstChannelInView (geo::View_t geoview, unsigned int apa, unsigned int cryo) const

uint32_t	FirstChannelInView (geo::View_t geoview, uint32_t chan) const

uint32_t	FirstChannelInView (uint32_t chan) const

unsigned int	ChannelsInAPAView (APAView_t apaview) const

unsigned int	ChannelsPerAPA () const

Private Attributes
art::ServiceHandle< geo::Geometry const >	fGeom

geo::WireReadoutGeom const *	fWireReadoutGeom

unsigned int	fChannelsPerAPA
	All APAs have this same number of channels. More...

unsigned int	fAPAsPerCryo

uint32_t	fFirstU

uint32_t	fLastU

uint32_t	fFirstV

uint32_t	fLastV

uint32_t	fFirstZ0

uint32_t	fLastZ0

uint32_t	fFirstZ1

uint32_t	fLastZ1

double	fChannelRange [2]

Detailed Description

Definition at line 27 of file APAGeometryAlg.h.

Constructor & Destructor Documentation

apa::APAGeometryAlg::APAGeometryAlg ( )

Definition at line 32 of file APAGeometryAlg.cxx.

References geo::WireReadoutGeom::ChannelToWire(), fAPAsPerCryo, fChannelRange, fChannelsPerAPA, fFirstU, fFirstV, fFirstZ0, fFirstZ1, fGeom, fLastU, fLastV, fLastZ0, fLastZ1, fWireReadoutGeom, Get, geo::kU, geo::kV, geo::kZ, geo::WireReadoutGeom::Nchannels(), geo::GeometryCore::NTPC(), geo::WireReadoutGeom::Plane(), geo::TPCID::TPC, geo::WireReadoutGeom::View(), and geo::PlaneGeo::WirePitch().

     : fWireReadoutGeom{&art::ServiceHandle<geo::WireReadout const>()->Get()}
   {
     // find the number of channels per APA
     uint32_t channel = 0;
 
     fChannelsPerAPA = fWireReadoutGeom->Nchannels();
     for (channel = 0; channel < fWireReadoutGeom->Nchannels(); ++channel) {
       if (fWireReadoutGeom->ChannelToWire(channel)[0].TPC > 1) {
         fChannelsPerAPA = channel;
         break;
       }
     }
 
     // Step through channel c and find the view boundaries, until
     // outside of first APA - these help optimize ChannelToAPAView
     // (very dependent on the conventions implimented in the channel map)
     fFirstU = 0;
     uint32_t c = 1;
     geo::WireID wid = fWireReadoutGeom->ChannelToWire(c)[0];
     geo::WireID lastwid;
     while (wid.TPC < 2) {
       if (fWireReadoutGeom->View(c) == geo::kV && fWireReadoutGeom->View(c - 1) == geo::kU) {
         fLastU = c - 1;
         fFirstV = c;
       }
 
       if (fWireReadoutGeom->View(c) == geo::kZ && fWireReadoutGeom->View(c - 1) == geo::kV) {
         fLastV = c - 1;
         fFirstZ0 = c;
       }
 
       if (wid.TPC == lastwid.TPC + 1) {
         fLastZ0 = c - 1;
         fFirstZ1 = c;
       }
 
       lastwid = wid;
       c++;
       if (c >= fWireReadoutGeom->Nchannels()) break;
       wid = fWireReadoutGeom->ChannelToWire(c)[0]; // for the while condition
     }
 
     fLastZ1 = c - 1;
 
     if (fLastZ1 + 1 != fChannelsPerAPA)
       throw cet::exception("APAGeometryAlg") << "Channel boundaries are inconsistent.\n";
 
     // some other things that will be needed
     fAPAsPerCryo = fGeom->NTPC() / 2;
     constexpr geo::TPCID tpcid{0, 0};
     fChannelRange[0] =
       (fLastU - fFirstU + 1) * fWireReadoutGeom->Plane({tpcid, geo::kU}).WirePitch();
     fChannelRange[1] =
       (fLastV - fFirstV + 1) * fWireReadoutGeom->Plane({tpcid, geo::kV}).WirePitch();
   }

Member Function Documentation

bool apa::APAGeometryAlg::APAChannelsIntersect	(	uint32_t	chan1,
		uint32_t	chan2,
		std::vector< geo::WireIDIntersection > &	IntersectVector
	)		const

If the channels intersect, get all intersections.

Definition at line 403 of file APAGeometryAlg.cxx.

References ChannelToAPA(), geo::WireReadoutGeom::ChannelToWire(), fWireReadoutGeom, geo::WireReadoutGeom::View(), and geo::WireReadoutGeom::WireIDsIntersect().

Referenced by ThreeChanPos(), and apa::DisambigAlg::UseEndPts().

   {
     // Get the WireIDs and view for each channel, make sure views are different
     std::vector<geo::WireID> wids1 = fWireReadoutGeom->ChannelToWire(chan1);
     std::vector<geo::WireID> wids2 = fWireReadoutGeom->ChannelToWire(chan2);
     geo::View_t view1 = fWireReadoutGeom->View(chan1);
     geo::View_t view2 = fWireReadoutGeom->View(chan2);
     if (view1 == view2) {
       mf::LogWarning("APAChannelsIntersect")
         << "Comparing two channels in the same view, return false";
       return false;
     }
     if (wids1[0].Cryostat != wids2[0].Cryostat || ChannelToAPA(chan1) != ChannelToAPA(chan2)) {
       throw cet::exception("APAChannelsIntersect")
         << "Comparing two channels in in different APAs: "
         << "channel " << chan1 << " in Cryo " << wids1[0].Cryostat << ", APA "
         << ChannelToAPA(chan1) << ", and channel " << chan2 << " in Cryo " << wids2[0].Cryostat
         << ", APA " << ChannelToAPA(chan2) << "\n";
       return false;
     }
 
     // Loop through wids1 and see if wids2 has any intersecting wires,
     // given that the WireIDs are in the same TPC
     for (unsigned int i1 = 0; i1 < wids1.size(); i1++) {
       for (unsigned int i2 = 0; i2 < wids2.size(); i2++) {
 
         // make sure it is reasonable to intersect
         if (wids1[i1].Plane == wids2[i2].Plane ||
             wids1[i1].TPC != wids2[i2].TPC || // not reasonable for a *WireID*
             wids1[i1].Cryostat != wids2[i2].Cryostat)
           continue;
 
         // Check if they even intersect; if they do, push back
         if (auto widIntersect = fWireReadoutGeom->WireIDsIntersect(wids1[i1], wids2[i2])) {
           IntersectVector.push_back(*widIntersect);
         }
       }
     }
 
     // Of all considered configurations, there are never more than
     // 4 intersections per channel pair
     if (IntersectVector.size() > 4) {
       mf::LogWarning("APAChannelsIntersect")
         << "Got " << IntersectVector.size() << " intersections for channels " << chan1 << " and "
         << chan2 << " - never expect more than 4, so far";
     }
 
     // With increasing IntersectVector index, the WireID vector indices of the
     // intersecting wireIDs increase.  This matches the direction in which ChannelToWire
     // builds its output WireID vector in the APA/35t Alg
     std::sort(IntersectVector.begin(), IntersectVector.end());
 
     // return true if any intersection points were found
     return not IntersectVector.empty();
   }

APAView_t apa::APAGeometryAlg::APAView ( uint32_t chan ) const

Get which of the 4 APA views the channel is in.

Definition at line 168 of file APAGeometryAlg.cxx.

References fChannelsPerAPA, fFirstZ1, fLastV, fLastZ0, fWireReadoutGeom, apa::kU, geo::kU, apa::kUnknown, apa::kV, geo::kV, geo::kZ, apa::kZ0, apa::kZ1, and geo::WireReadoutGeom::View().

   {
     // it seems trivial to do this for U and V, but this gives a side to
     // geo::kZ, unlike Geometry::View(c), as is often needed in disambiguation
 
     geo::View_t view = fWireReadoutGeom->View(chan);
     switch (view) {
     default: break;
     case geo::kU: return kU;
     case geo::kV: return kV;
     case geo::kZ:
       unsigned int modchan = chan % fChannelsPerAPA;
       // Channel mapping number in the order of U, V, Z0, then Z1
       if (modchan > fLastV && modchan < fFirstZ1) return kZ0;
       if (modchan > fLastZ0) return kZ1;
     }
 
     return kUnknown;
   }

unsigned int apa::APAGeometryAlg::ChannelsInAPAView ( APAView_t apaview ) const

Definition at line 124 of file APAGeometryAlg.cxx.

References fFirstU, fFirstV, fFirstZ0, fFirstZ1, fLastZ1, apa::kU, apa::kV, apa::kZ0, and apa::kZ1.

Referenced by ChannelsInView().

   {
     switch (apaview) {
     default: return 0;
     case kU: return fFirstV - fFirstU;
     case kV: return fFirstZ0 - fFirstV;
     case kZ0: return fFirstZ1 - fFirstZ0;
     case kZ1: return fLastZ1 - fFirstZ1 + 1;
     }
   }

unsigned int apa::APAGeometryAlg::ChannelsInView ( geo::View_t geoview ) const

Definition at line 108 of file APAGeometryAlg.cxx.

References ChannelsInAPAView(), apa::kU, geo::kU, apa::kV, geo::kV, geo::kZ, apa::kZ0, and apa::kZ1.

Referenced by apa::DisambigAlg::FindChanTimeEndPts(), and apa::DisambigAlg::MakeCloseHits().

   {
     switch (geoview) {
     default: return 0;
     case geo::kU: return ChannelsInAPAView(kU);
     case geo::kV: return ChannelsInAPAView(kV);
     case geo::kZ:
       if (ChannelsInAPAView(kZ0) == ChannelsInAPAView(kZ1))
         return ChannelsInAPAView(kZ0);
       else
         throw cet::exception("ChannelsInView")
           << "Both Z sides should have the same amount of channels\n";
     }
   }

unsigned int apa::APAGeometryAlg::ChannelsPerAPA ( ) const

inline

Definition at line 63 of file APAGeometryAlg.h.

References fChannelsPerAPA.

63 { return fChannelsPerAPA; }

apa::APAGeometryAlg::fChannelsPerAPA

unsigned int fChannelsPerAPA

All APAs have this same number of channels.

Definition: APAGeometryAlg.h:69

unsigned int apa::APAGeometryAlg::ChannelToAPA ( uint32_t chan ) const

Get number of the APA containing the given channel.

Definition at line 102 of file APAGeometryAlg.cxx.

References fChannelsPerAPA.

Referenced by APAChannelsIntersect(), FirstChannelInView(), apa::DisambigAlg::MakeCloseHits(), and apa::DisambigAlg::RunDisambig().

   {
     return chan / fChannelsPerAPA;
   }

void apa::APAGeometryAlg::ChannelToAPA	(	uint32_t	chan,
		unsigned int &	apa,
		unsigned int &	cryo
	)		const

Definition at line 90 of file APAGeometryAlg.cxx.

References fAPAsPerCryo, and fChannelsPerAPA.

   {
     cryo = chan / (fAPAsPerCryo * fChannelsPerAPA);
 
     // Number apa uniquely across cryostats so that
     // apa to recob::Object maps are easy to work with.
     // If we decide to reset APA number per cryo, uncomment:
     //chan -= cryo*fAPAsPerCryo*fChannelsPerAPA;
     apa = chan / fChannelsPerAPA;
   }

std::vector< geo::WireID > apa::APAGeometryAlg::ChanSegsPerSide	(	uint32_t	chan,
		unsigned int	side
	)		const

Definition at line 189 of file APAGeometryAlg.cxx.

References geo::WireReadoutGeom::ChannelToWire(), and fWireReadoutGeom.

   {
     std::vector<geo::WireID> wids = fWireReadoutGeom->ChannelToWire(chan);
     return ChanSegsPerSide(wids, side);
   }

std::vector< geo::WireID > apa::APAGeometryAlg::ChanSegsPerSide	(	std::vector< geo::WireID >	wids,
		unsigned int	side
	)		const

Definition at line 196 of file APAGeometryAlg.cxx.

   {
     // Given a vector of wireIDs and an APA side, return
     // the wireIDs the the tpc side where tpc%2 = side
 
     std::vector<geo::WireID> thisSide;
 
     for (size_t i = 0; i < wids.size(); i++)
       if (wids[i].TPC % 2 == side) thisSide.push_back(wids[i]);
 
     return thisSide;
   }

uint32_t apa::APAGeometryAlg::FirstChannelInView	(	geo::View_t	geoview,
		unsigned int	apa,
		unsigned int	cryo
	)		const

Definition at line 136 of file APAGeometryAlg.cxx.

References fAPAsPerCryo, fChannelsPerAPA, fFirstU, fFirstV, fFirstZ0, geo::kU, geo::kV, and geo::kZ.

Referenced by apa::DisambigAlg::FindChanTimeEndPts(), FirstChannelInView(), and apa::DisambigAlg::MakeCloseHits().

   {
     switch (geoview) {
     default: return 0 + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
     case geo::kU: return fFirstU + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
     case geo::kV: return fFirstV + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
     case geo::kZ:
       //TODO: would need tpc number for the rest of this
       return fFirstZ0 + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
     }
   }

uint32_t apa::APAGeometryAlg::FirstChannelInView	(	geo::View_t	geoview,
		uint32_t	chan
	)		const

Definition at line 160 of file APAGeometryAlg.cxx.

References ChannelToAPA(), and FirstChannelInView().

   {
     unsigned int apa, cryo;
     ChannelToAPA(chan, apa, cryo);
     return FirstChannelInView(geoview, apa, cryo);
   }

uint32_t apa::APAGeometryAlg::FirstChannelInView ( uint32_t chan ) const

Definition at line 151 of file APAGeometryAlg.cxx.

References ChannelToAPA(), FirstChannelInView(), fWireReadoutGeom, and geo::WireReadoutGeom::View().

   {
     geo::View_t geoview = fWireReadoutGeom->View(chan);
     unsigned int apa, cryo;
     ChannelToAPA(chan, apa, cryo);
     return FirstChannelInView(geoview, chan);
   }

bool apa::APAGeometryAlg::LineSegChanIntersect	(	geo::Point_t const &	xyzStart,
		geo::Point_t const &	xyzEnd,
		uint32_t	chan,
		std::vector< geo::WireID > &	widsCrossed,
		bool	ExtendLine = `true`
	)		const

If a line given by start/end points intersects a channel.

Definition at line 244 of file APAGeometryAlg.cxx.

References geo::WireReadoutGeom::ChannelToWire(), fGeom, fWireReadoutGeom, geo::vect::middlePoint(), geo::PlaneGeo::NearestWireID(), geo::WireReadoutGeom::Plane(), geo::GeometryCore::PositionToTPCID(), lar::util::ValueInRange(), w, and geo::WireID::Wire.

   {
     // This assumes a smooth wire numbering, and that the line seg is contained in a tpc.
     // Meant for use with the approximate line calculated by matching cluster endpoints in
     // disambiguation.
 
     // Find tpc, use midpoint in case start/end is on a boundary
     geo::Point_t const xyzMid = geo::vect::middlePoint({xyzStart, xyzEnd});
     auto const& tpcid = fGeom->PositionToTPCID(xyzMid);
 
     // Find the nearest wire number to the line segment endpoints
     std::vector<geo::WireID> wids = fWireReadoutGeom->ChannelToWire(chan);
     geo::PlaneID const planeID{tpcid, wids[0].Plane};
     auto const& plane = fWireReadoutGeom->Plane(planeID);
     unsigned int startW = plane.NearestWireID(xyzStart).Wire;
     unsigned int endW = plane.NearestWireID(xyzEnd).Wire;
 
     if (startW > endW) std::swap(startW, endW);
 
     // Loop through wireIDs and check for intersection, if in the right TPC
     for (size_t w = 0; w < wids.size(); w++) {
       if (wids[w].TPC != tpcid.TPC) continue;
       if (wids[w].Cryostat != tpcid.Cryostat)
         throw cet::exception("LineSegChanIntersect")
           << "Channel and line not in the same crostat.\n";
 
       // If the current wire id wire number is inbetween the start/end
       // point wires, the line segment intersects the wireID at some point.
 
       // TODO: for now, extend range, but that is application specific. fix asap
       // The longer we make the range, the more conservative it is, so it is safe
       // to extend the range a bit to get hits at the ends of the line
       unsigned int ext = 0;
       if (ExtendLine) ext = 10;
 
       if (lar::util::ValueInRange(wids[w].Wire * 1., (startW - ext) * 1., (endW + ext) * 1.))
         widsCrossed.push_back(wids[w]);
     }
 
     return not widsCrossed.empty();
   }

geo::WireID apa::APAGeometryAlg::NearestWireIDOnChan	(	geo::Point_t const &	WorldLoc,
		uint32_t	chan,
		geo::PlaneID const &	planeID
	)		const

Definition at line 211 of file APAGeometryAlg.cxx.

References util::abs(), geo::WireReadoutGeom::ChannelToWire(), geo::CryostatID::Cryostat, fChannelRange, fWireReadoutGeom, geo::kZ, geo::PlaneGeo::NearestWireID(), geo::WireReadoutGeom::Plane(), geo::PlaneID::Plane, geo::TPCID::TPC, geo::WireReadoutGeom::View(), and geo::WireID::Wire.

Referenced by ThreeChanPos(), and apa::DisambigAlg::UseEndPts().

   {
     std::vector<geo::WireID> cWids = fWireReadoutGeom->ChannelToWire(chan);
 
     if (cWids[0].Cryostat != planeID.Cryostat)
       throw cet::exception("APAGeometryAlg")
         << "Channel " << chan << "not in cryostat " << planeID.Cryostat << "\n";
 
     if (std::floor(cWids[0].TPC / 2) != std::floor(planeID.TPC / 2))
       throw cet::exception("APAGeometryAlg")
         << "Channel " << chan << "not in APA " << std::floor(planeID.TPC / 2) << "\n";
 
     // special case for vertical wires
     if (fWireReadoutGeom->View(chan) == geo::kZ) {
       return fWireReadoutGeom->ChannelToWire(chan)[0];
     }
 
     unsigned int xyzWire = fWireReadoutGeom->Plane(planeID).NearestWireID(worldLoc).Wire;
 
     // The desired wire ID will be the only channel segment within half the channel range.
     geo::WireID wid;
     for (size_t i = 0; i < cWids.size(); i++) {
       if (cWids[i].TPC != planeID.TPC) continue;
       if (std::abs((int)cWids[i].Wire - (int)xyzWire) < fChannelRange[planeID.Plane] / 2)
         wid = cWids[i];
     }
 
     return wid;
   }

std::vector< double > apa::APAGeometryAlg::ThreeChanPos	(	uint32_t	u,
		uint32_t	v,
		uint32_t	z
	)		const

Find the center of the 3 intersections, choose best if multiple.

Definition at line 291 of file APAGeometryAlg.cxx.

References util::abs(), APAChannelsIntersect(), geo::WireReadoutGeom::ChannelToWire(), fGeom, fWireReadoutGeom, geo::WireGeo::GetCenter(), geo::TPCGeo::Length(), NearestWireIDOnChan(), geo::GeometryCore::TPC(), geo::WireReadoutGeom::Wire(), geo::WireReadoutGeom::WireIDsIntersect(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by apa::DisambigAlg::UseEndPts().

   {
     // Say we've associated a U, V, and Z channel -- perhaps by associating hits
     // or cluster endpoints -- these don't necessarily all intersect, but they
     // are hopefully pretty close. Find the center of the 3 intersections.
 
     // get data needed along the way
     std::vector<geo::WireIDIntersection> UVIntersects;
     APAChannelsIntersect(u, v, UVIntersects);
     std::vector<double> UVzToZ(UVIntersects.size());
     geo::WireID Zwid = fWireReadoutGeom->ChannelToWire(z)[0];
     unsigned int cryo = Zwid.Cryostat;
     unsigned int tpc = Zwid.TPC;
     std::vector<geo::WireID> Uwids = fWireReadoutGeom->ChannelToWire(u);
     std::vector<geo::WireID> Vwids = fWireReadoutGeom->ChannelToWire(v);
     std::vector<geo::WireID> UwidsInTPC, VwidsInTPC;
     for (size_t i = 0; i < Uwids.size(); i++)
       if (Uwids[i].TPC == tpc) UwidsInTPC.push_back(Uwids[i]);
     for (size_t i = 0; i < Vwids.size(); i++)
       if (Vwids[i].TPC == tpc) VwidsInTPC.push_back(Vwids[i]);
     auto const Zcent = fWireReadoutGeom->Wire(Zwid).GetCenter();
 
     std::cout << "Zcent = " << Zcent.Z() << ", UVintersects zpos = ";
     for (size_t uv = 0; uv < UVIntersects.size(); uv++) {
       std::cout << UVIntersects[uv].z << ", ";
     }
     std::cout << "\n";
 
 
     if (UVIntersects.empty()) {
       if (UwidsInTPC.size() > 1 || VwidsInTPC.size() > 1)
         throw cet::exception("ThreeChanPos") << "U/V channels don't intersect, bad return.\n";
 
       // Now assume there are only one of each u and v wireIDs on in this TPC
       mf::LogWarning("ThreeChanPos")
         << "No U/V intersect, exceptional channels. See if U or V intersects Z\n";
       std::vector<double> yzCenter(2, 0.);
       geo::WireID Uwid = UwidsInTPC[0];
       geo::WireID Vwid = VwidsInTPC[0];
       auto const UZInt = fWireReadoutGeom->WireIDsIntersect(Uwid, Zwid);
       auto const VZInt = fWireReadoutGeom->WireIDsIntersect(Vwid, Zwid);
       if (!UZInt && !VZInt)
         throw cet::exception("NoChanIntersect") << "No channels intersect, bad return.\n";
       if (UZInt && !VZInt) {
         yzCenter[0] = UZInt->y;
         yzCenter[1] = UZInt->z;
       }
       if (VZInt && !UZInt) {
         yzCenter[0] = VZInt->y;
         yzCenter[1] = VZInt->z;
       }
       if (UZInt && VZInt) {
         yzCenter[0] = (VZInt->y + UZInt->y) / 2;
         yzCenter[1] = (VZInt->z + UZInt->z) / 2;
       }
       return yzCenter;
     }
 
 
     // In case the uv channels intersect twice on the same side, choose the best case.
     // Note: this will not happen for APAs with UV angle at about 36, but will for 45
     std::cout << "UVzToZ = ";
     for (size_t widI = 0; widI < UVIntersects.size(); widI++) {
       UVzToZ[widI] = std::abs(UVIntersects[widI].z - Zcent.Z());
       std::cout << UVzToZ[widI] << ", ";
     }
     std::cout << "\n";
 
     unsigned int bestWidI = 0;
     double minZdiff = fGeom->TPC(Zwid).Length(); // start it out at maximum z
     for (unsigned int widI = 0; widI < UVIntersects.size(); widI++) {
       if (UVIntersects[widI].TPC == tpc && UVzToZ[widI] < minZdiff) {
         minZdiff = UVzToZ[widI];
         bestWidI = widI;
       }
     }
     geo::WireIDIntersection ChosenUVInt = UVIntersects[bestWidI];
 
     // Now having the UV intersection, get the UZ and VZ
     geo::Point_t const UVInt{0., ChosenUVInt.y, ChosenUVInt.z};
     geo::WireID Uwid = NearestWireIDOnChan(UVInt, u, {cryo, tpc, 0});
     geo::WireID Vwid = NearestWireIDOnChan(UVInt, v, {cryo, tpc, 1});
     auto const UZInt = fWireReadoutGeom->WireIDsIntersect(Uwid, Zwid);
     auto const VZInt = fWireReadoutGeom->WireIDsIntersect(Vwid, Zwid);
 
     // find the center
     std::vector<double> yzCenter(2, 0.);
 
     std::cout << "UZint = " << std::boolalpha << static_cast<bool>(UZInt)
               << ", VZint = " << static_cast<bool>(VZInt) << '\n';
     if (!UZInt || !VZInt) {
       std::cout << "ChosenUVint.y = " << ChosenUVInt.y << "ChosenUVint.z = " << ChosenUVInt.z
                 << std::endl;
 
       //temporary case
       yzCenter[0] = ChosenUVInt.y;
       yzCenter[1] = ChosenUVInt.z;
     }
     else {
       std::cout << "UZint.z = " << UZInt->z << ", VZint.z = " << VZInt->z << '\n';
       yzCenter[0] = (ChosenUVInt.y + UZInt->y + VZInt->y) / 3;
       yzCenter[1] = (ChosenUVInt.z + UZInt->z + VZInt->z) / 3;
     }
 
     return yzCenter;
   }