#include "APAGeometryAlg.h"

Public Member Functions
	APAGeometryAlg (fhicl::ParameterSet const &pset)

	APAGeometryAlg ()

virtual	~APAGeometryAlg ()

void	reconfigure (fhicl::ParameterSet const &p)

void	Init ()
	Initialize some chanel numbers to speed up other methods. More...

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

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

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

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

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

geo::WireID	NearestWireIDOnChan (const double WorldLoc[3], uint32_t chan, unsigned int const plane, unsigned int const tpc=0, unsigned int const cstat=0)

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

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

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

unsigned int	ChannelsInView (geo::View_t geoview)

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

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

uint32_t	FirstChannelInView (uint32_t chan)

unsigned int	ChannelsInAPAView (APAView_t apaview)

unsigned int	ChannelsPerAPA ()

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

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 41 of file APAGeometryAlg.h.

Constructor & Destructor Documentation

apa::APAGeometryAlg::APAGeometryAlg ( fhicl::ParameterSet const & pset )

Definition at line 39 of file APAGeometryAlg.cxx.

References Init(), and reconfigure().

   {
     this->reconfigure(pset); 
     this->Init();
   }

apa::APAGeometryAlg::APAGeometryAlg ( )

Definition at line 47 of file APAGeometryAlg.cxx.

References Init().

   {
     this->Init();
   }

apa::APAGeometryAlg::~APAGeometryAlg ( )

virtual

Definition at line 53 of file APAGeometryAlg.cxx.

54 {

55 }

Member Function Documentation

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

If the channels intersect, get all intersections.

Definition at line 478 of file APAGeometryAlg.cxx.

References ChannelToAPA(), geo::GeometryCore::ChannelToWire(), fGeom, geo::GeometryCore::View(), and geo::GeometryCore::WireIDsIntersect().

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

   {
 
     
     // Get the WireIDs and view for each channel, make sure views are different
     geo::WireIDIntersection    widIntersect;
     std::vector< geo::WireID > wids1 = fGeom->ChannelToWire( chan1 );
     std::vector< geo::WireID > wids2 = fGeom->ChannelToWire( chan2 );
     geo::View_t                view1 = fGeom->View( chan1 );
     geo::View_t                view2 = fGeom->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         ||
         this->ChannelToAPA(chan1) != this->ChannelToAPA(chan2)    ){
       throw cet::exception("APAChannelsIntersect") << "Comparing two channels in in different APAs: "
                                                    << "channel " << chan1 << " in Cryo " 
                                                    << wids1[0].Cryostat << ", APA " << this->ChannelToAPA(chan1)
                                                    << ", and channel " << chan2 << " in Cryo " 
                                                    << wids2[0].Cryostat << ", APA " << this->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;
 
 //      std::cout << "Checking: \n WireID 1 = (" 
 //                << wids1[i1].Cryostat << "," << wids1[i1].TPC << ","
 //                << wids1[i1].Plane    << "," << wids1[i1].Wire
 //                << ") \n WireID 2 = ("
 //                << wids2[i2].Cryostat << "," << wids2[i2].TPC << ","
 //                << wids2[i2].Plane    << "," << wids2[i2].Wire << ")" << std::endl;
 
         // Check if they even intersect; if they do, push back
         if( fGeom->WireIDsIntersect( wids1[i1], wids2[i2], widIntersect ) ){
 
 //        std::cout << "we have an intersect" << std::endl;
 
           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
     if( IntersectVector.size() == 0 ) return false;
     else return true;
 
   }

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

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

Definition at line 217 of file APAGeometryAlg.cxx.

References fChannelsPerAPA, fFirstZ1, fGeom, fLastV, fLastZ0, apa::kU, geo::kU, apa::kUnknown, apa::kV, geo::kV, geo::kZ, apa::kZ0, apa::kZ1, and geo::GeometryCore::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 = fGeom->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 )

Definition at line 156 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 )

Definition at line 136 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 ( )

inline

Definition at line 90 of file APAGeometryAlg.h.

90 { return fChannelsPerAPA; };

apa::APAGeometryAlg::fChannelsPerAPA

unsigned int fChannelsPerAPA

All APAs have this same number of channels.

Definition: APAGeometryAlg.h:97

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

Get number of the APA containing the given channel.

Definition at line 129 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
	)

Definition at line 113 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;
 
     return;
   }

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

Definition at line 243 of file APAGeometryAlg.cxx.

References geo::GeometryCore::ChannelToWire(), and fGeom.

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

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

Definition at line 253 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
	)

Definition at line 175 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
	)

Definition at line 205 of file APAGeometryAlg.cxx.

References ChannelToAPA(), and FirstChannelInView().

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

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

Definition at line 195 of file APAGeometryAlg.cxx.

References ChannelToAPA(), fGeom, FirstChannelInView(), and geo::GeometryCore::View().

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

void apa::APAGeometryAlg::Init ( )

Initialize some chanel numbers to speed up other methods.

Definition at line 63 of file APAGeometryAlg.cxx.

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

Referenced by APAGeometryAlg().

   {
 
     // find the number of channels per APA
     uint32_t channel = 0;
     while( fGeom->ChannelToWire(channel+1)[0].TPC < 2 ) channel++;
     fChannelsPerAPA = channel + 1;
 
     // 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 = (fGeom->ChannelToWire(c))[0];
     geo::WireID lastwid;
     while( wid.TPC < 2 ){
 
       if( fGeom->View(c) == geo::kV && fGeom->View(c-1) == geo::kU ){
         fLastU = c-1;
         fFirstV = c;   }
       
       if( fGeom->View(c) == geo::kZ && fGeom->View(c-1) == geo::kV ){
         fLastV = c-1;
         fFirstZ0 = c;  }
       
       if( wid.TPC == lastwid.TPC + 1  ){
         fLastZ0 = c-1;
         fFirstZ1 = c;  }
 
       lastwid = wid;
       c++;
       wid = (fGeom->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(0)/2;
     fChannelRange[0] = (fLastU-fFirstU + 1)*fGeom->WirePitch(geo::kU);
     fChannelRange[1] = (fLastV-fFirstV + 1)*fGeom->WirePitch(geo::kV);
 
   }

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

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

Definition at line 305 of file APAGeometryAlg.cxx.

References geo::GeometryCore::ChannelToWire(), fGeom, geo::GeometryCore::NearestWire(), geo::GeometryCore::PositionToTPC(), geo::GeometryCore::ValueInRange(), and w.

   {
 
     // 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
     unsigned int tpc, cryo;
     double xyzMid[3];
     xyzMid[0] = (xyzStart[0]+xyzEnd[0])/2;
     xyzMid[1] = (xyzStart[1]+xyzEnd[1])/2;
     xyzMid[2] = (xyzStart[2]+xyzEnd[2])/2;
     fGeom->PositionToTPC(xyzMid, tpc, cryo);
 
     // Find the nearest wire number to the line segment endpoints
     std::vector<geo::WireID> wids = fGeom->ChannelToWire(chan);
     unsigned int startW = fGeom->NearestWire( xyzStart, wids[0].Plane, tpc, cryo );
     unsigned int endW   = fGeom->NearestWire( xyzEnd,   wids[0].Plane, tpc, cryo );
 
     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 != tpc ) continue;
       if( wids[w].Cryostat != cryo ) 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( fGeom->ValueInRange( wids[w].Wire*1., (startW-ext)*1., (endW+ext)*1. ) ) widsCrossed.push_back(wids[w]);
 
     }
 
     if( widsCrossed.size() > 0 ) return true;
     else return false;
 
   }

geo::WireID apa::APAGeometryAlg::NearestWireIDOnChan	(	const double	WorldLoc[3],
		uint32_t	chan,
		unsigned int const	plane,
		unsigned int const	tpc = `0`,
		unsigned int const	cstat = `0`
	)

Definition at line 270 of file APAGeometryAlg.cxx.

References geo::GeometryCore::ChannelToWire(), fChannelRange, fGeom, geo::kZ, geo::GeometryCore::NearestWireID(), geo::GeometryCore::View(), and geo::WireID::Wire.

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

   {
 
     std::vector<geo::WireID> cWids = fGeom->ChannelToWire( chan );
 
     if( cWids[0].Cryostat != cstat ) 
       throw cet::exception("APAGeometryAlg") << "Channel " << chan 
                                              << "not in cryostat " << cstat << "\n";
     if( std::floor( cWids[0].TPC / 2 ) != std::floor( tpc / 2 ) ) 
       throw cet::exception("APAGeometryAlg") << "Channel " << chan 
                                              << "not in APA " << std::floor(tpc/2) << "\n";
 
     // special case for vertical wires
     if(fGeom->View(chan)==geo::kZ) return fGeom->ChannelToWire(chan)[0];
 
     unsigned int xyzWire = fGeom->NearestWireID( WorldLoc, plane, tpc, cstat ).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 != tpc ) continue;
       if( std::abs((int)cWids[i].Wire - (int)xyzWire) < fChannelRange[plane]/2 ) wid=cWids[i];
     }
     
     return wid;
 
   }

void apa::APAGeometryAlg::reconfigure ( fhicl::ParameterSet const & p )

Definition at line 58 of file APAGeometryAlg.cxx.

Referenced by APAGeometryAlg().

59 {

60 }

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

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

Definition at line 357 of file APAGeometryAlg.cxx.

References APAChannelsIntersect(), geo::GeometryCore::ChannelToWire(), geo::GeometryCore::Cryostat(), fGeom, geo::WireGeo::GetCenter(), geo::TPCGeo::Length(), NearestWireIDOnChan(), geo::CryostatGeo::TPC(), geo::GeometryCore::WireIDsIntersect(), geo::GeometryCore::WireIDToWireGeo(), 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;
     this->APAChannelsIntersect( u, v, UVIntersects );
     std::vector< double > UVzToZ(UVIntersects.size());
     geo::WireID Zwid = fGeom->ChannelToWire(z)[0];
     unsigned int cryo = Zwid.Cryostat;
     unsigned int tpc = Zwid.TPC;
     std::vector<geo::WireID> Uwids = fGeom->ChannelToWire(u);
     std::vector<geo::WireID> Vwids = fGeom->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]);
     double Zcent[3] = {0.};
     fGeom->WireIDToWireGeo( Zwid ).GetCenter(Zcent);
 
     std::cout << "Zcent = " << Zcent[2] << ", UVintersects zpos = ";
     for(size_t uv=0; uv<UVIntersects.size(); uv++){
       std::cout << UVIntersects[uv].z << ", ";
     }
     std::cout << "\n";
 
 
     //std::cout << "U = " << u << " V = " << v << ", " << UVIntersects.size() << std::endl;
 
     if( UVIntersects.size() == 0 ){
       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];
       geo::WireIDIntersection UZInt, VZInt;
       bool checkUZ = fGeom->WireIDsIntersect( Uwid, Zwid, UZInt );
       bool checkVZ = fGeom->WireIDsIntersect( Vwid, Zwid, VZInt );
       if( !checkUZ && !checkVZ )
         throw cet::exception("NoChanIntersect") << "No channels intersect, bad return.\n";
       if( checkUZ && !checkVZ ){ yzCenter[0] = UZInt.y; yzCenter[1] = UZInt.z; }
       if( checkVZ && !checkUZ ){ yzCenter[0] = VZInt.y; yzCenter[1] = VZInt.z; }
       if( checkUZ &&  checkVZ ){ 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[2] );
       std::cout << UVzToZ[widI] << ", ";
     }
     std::cout<<"\n";
 
     unsigned int bestWidI = 0;
     double minZdiff = fGeom->Cryostat(cryo).TPC(tpc).Length(); // start it out at maximum z
     for( unsigned int widI = 0; widI < UVIntersects.size(); widI++ ){
 
       //std::cout << "widI = " << widI << std::endl; 
 
       if( UVIntersects[widI].TPC == tpc  &&  UVzToZ[widI] < minZdiff ){ 
         minZdiff = UVzToZ[widI];
         bestWidI = widI;
         //std::cout << "bestWidI = " << bestWidI << std::endl;
       }
     }
     geo::WireIDIntersection ChosenUVInt = UVIntersects[bestWidI];
     
     // Now having the UV intersection, get the UZ and VZ
     double UVInt[3] = {0.};
     UVInt[1] = ChosenUVInt.y; UVInt[2] = ChosenUVInt.z;
     geo::WireID Uwid = this->NearestWireIDOnChan( UVInt, u, 0, tpc, cryo );
     geo::WireID Vwid = this->NearestWireIDOnChan( UVInt, v, 1, tpc, cryo );
     geo::WireIDIntersection UZInt, VZInt;
     bool checkUZ = fGeom->WireIDsIntersect( Uwid, Zwid, UZInt );
     bool checkVZ = fGeom->WireIDsIntersect( Vwid, Zwid, VZInt );
 
     std::cout << "UZint.z = " << UZInt.z << " (" << checkUZ << "), VZint.z = " << VZInt.z << " (" << checkVZ << ")\n";
     
     // find the center
     std::vector<double> yzCenter(2,0.);
 
 
     if( !checkUZ || !checkVZ ){
       //throw cet::exception("ThreeChanPos")  << "WireIDs were expected to intersect.\n";
 
       std::cout << "ChosenUVint.y = " << ChosenUVInt.y << "ChosenUVint.z = " << ChosenUVInt.z << std::endl;
       
       //temporary case
       yzCenter[0] = ChosenUVInt.y;
       yzCenter[1] = ChosenUVInt.z;
 
 
     } else {
 
     yzCenter[0] =  (ChosenUVInt.y + UZInt.y + VZInt.y)/3; 
     yzCenter[1] =  (ChosenUVInt.z + UZInt.z + VZInt.z)/3;
 
     }
 
 
     return yzCenter;
 
   } 