#include "DisambigAlg.h"

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

void	RunDisambig (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::Handle< std::vector< recob::Hit >> GausHits)

void	TrivialDisambig (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
	Make the easiest and safest disambiguations in apa. More...

void	Crawl (unsigned int apa)
	Extend what we disambiguation we do have in apa. More...

unsigned int	FindChanTimeEndPts (detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
	Basic endpoint-hit finder per apa. More...

void	UseEndPts (detinfo::DetectorPropertiesData const &detProp, unsigned int apa)

unsigned int	CompareViews (detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
	Compare U and V to see if one says something about the other. More...

void	AssessDisambigSoFar (unsigned int apa)
	See how much disambiguation has been done in this apa so far. More...

Public Attributes
std::map< unsigned int, double >	fUeffSoFar

std::map< unsigned int, double >	fVeffSoFar

std::map< unsigned int, unsigned int >	fnUSoFar

std::map< unsigned int, unsigned int >	fnVSoFar

std::map< unsigned int, unsigned int >	fnDUSoFar

std::map< unsigned int, unsigned int >	fnDVSoFar

std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > >	fDisambigHits
	The final list of hits to pass back to be made. More...

Private Member Functions
void	MakeDisambigHit (art::Ptr< recob::Hit > const &hit, geo::WireID, unsigned int apa)
	Makes a disambiguated hit while keeping track of what has already been disambiguated. More...

unsigned int	MakeCloseHits (int ext, geo::WireID wid, double Dmin, double Dmax)
	Having disambiguated a time range on a wireID, extend to neighboring channels. More...

bool	HitsOverlapInTime (detinfo::DetectorPropertiesData const &detProp, recob::Hit const &hitA, recob::Hit const &hitB)

bool	HitsReasonablyMatch (art::Ptr< recob::Hit > hitA, art::Ptr< recob::Hit > hitB)

Private Attributes
apa::APAGeometryAlg	fAPAGeo

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

art::ServiceHandle< cheat::BackTrackerService const >	bt_serv

std::map< raw::ChannelID_t, std::vector< art::Ptr< recob::Hit > > >	fChannelToHits

std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > >	fAPAToUVHits

std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > >	fAPAToZHits

std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > >	fAPAToHits

std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > >	fAPAToEndPHits
	\ todo: Channel/APA to hits can be done in a unified way More...

std::map< unsigned int, std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > >	fAPAToDHits
	Hold the disambiguations per APA. More...

std::map< std::pair< double, double >, geo::WireID >	fChanTimeToWid
	If a hit is disambiguated, map its chan and peak time to the chosen wireID. More...

std::map< unsigned int, std::map< std::pair< double, double >, bool > >	fHasBeenDisambiged
	Convenient way to keep track of disambiguation so far. More...

bool	fCrawl
	\ todo: Write function that compares hits more detailedly More...

bool	fUseEndP

bool	fCompareViews

unsigned int	fNChanJumps
	Number of channels the crawl can jump over. More...

double	fCloseHitsRadius

double	fMaxEndPDegRange

Detailed Description

Definition at line 34 of file DisambigAlg.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 40 of file DisambigAlg.cxx.

References fCloseHitsRadius, fCompareViews, fCrawl, fMaxEndPDegRange, fNChanJumps, fUseEndP, and fhicl::ParameterSet::get().

     : fAPAGeo(p.get<fhicl::ParameterSet>("APAGeometryAlg"))
   {
     fCrawl = p.get<bool>("Crawl");
     fUseEndP = p.get<bool>("UseEndP");
     fCompareViews = p.get<bool>("CompareViews");
     fCloseHitsRadius = p.get<double>("CloseHitsRadius");
     fMaxEndPDegRange = p.get<double>("MaxEndPDegRange");
     fNChanJumps = p.get<unsigned int>("NChanJumps");
   }

Member Function Documentation

void apa::DisambigAlg::AssessDisambigSoFar ( unsigned int apa )

See how much disambiguation has been done in this apa so far.

Definition at line 633 of file DisambigAlg.cxx.

References fAPAToDHits, fAPAToUVHits, fnDUSoFar, fnDVSoFar, fnUSoFar, fnVSoFar, fUeffSoFar, fVeffSoFar, geo::kU, geo::kV, and recob::Hit::View().

Referenced by RunDisambig().

   {
     unsigned int nU(0), nV(0);
     for (size_t h = 0; h < fAPAToUVHits[apa].size(); h++) {
       art::Ptr<recob::Hit> hit = fAPAToUVHits[apa][h];
       if (hit->View() == geo::kU)
         nU++;
       else if (hit->View() == geo::kV)
         nV++;
     }
 
     unsigned int nDU(0), nDV(0);
     for (size_t h = 0; h < fAPAToDHits[apa].size(); h++) {
       art::Ptr<recob::Hit> hit = fAPAToDHits[apa][h].first;
       if (hit->View() == geo::kU)
         nDU++;
       else if (hit->View() == geo::kV)
         nDV++;
     }
 
     fUeffSoFar[apa] = (nDU * 1.) / (nU * 1.);
     fVeffSoFar[apa] = (nDV * 1.) / (nV * 1.);
     fnUSoFar[apa] = nU;
     fnVSoFar[apa] = nV;
     fnDUSoFar[apa] = nDU;
     fnDVSoFar[apa] = nDV;
   }

unsigned int apa::DisambigAlg::CompareViews	(	detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	apa
	)

Compare U and V to see if one says something about the other.

Definition at line 663 of file DisambigAlg.cxx.

References geo::GeometryCore::ChannelToWire(), d, fAPAToUVHits, fChanTimeToWid, fHasBeenDisambiged, geom, HitsOverlapInTime(), MakeDisambigHit(), lar::to_element, w, and geo::GeometryCore::WireIDsIntersect().

Referenced by RunDisambig().

   {
     unsigned int nDisambiguations(0);
 
     // loop through all hits that are still ambiguous
     for (auto const& ambighitPtr : fAPAToUVHits[apa]) {
       auto const& ambighit = *ambighitPtr;
       raw::ChannelID_t ambigchan = ambighit.Channel();
       std::pair<double, double> ambigChanTime(ambigchan * 1., ambighit.PeakTime());
       if (fHasBeenDisambiged[apa][ambigChanTime]) continue;
       geo::View_t view = ambighit.View();
       std::vector<geo::WireID> ambigwids = geom->ChannelToWire(ambigchan);
       std::vector<unsigned int> widDcounts(ambigwids.size(), 0);
       std::vector<unsigned int> widAcounts(ambigwids.size(), 0);
 
       // loop through hits in the other view which are close in time
       for (auto const& hit : fAPAToUVHits[apa] | transform(to_element)) {
         if (hit.View() == view || !this->HitsOverlapInTime(detProp, ambighit, hit)) continue;
 
         // An other-view-hit overlaps in time, see what
         // wids of the ambiguous hit's channels it overlaps
         raw::ChannelID_t chan = hit.Channel();
         std::vector<geo::WireID> wids = geom->ChannelToWire(chan);
         std::pair<double, double> ChanTime(chan * 1., hit.PeakTime());
         geo::WireIDIntersection widIntersect; // only so we can use the function
         if (fHasBeenDisambiged[apa][ChanTime]) {
           for (size_t a = 0; a < ambigwids.size(); a++)
             if (ambigwids[a].TPC == fChanTimeToWid[ChanTime].TPC &&
                 geom->WireIDsIntersect(ambigwids[a], fChanTimeToWid[ChanTime], widIntersect))
               widDcounts[a]++;
         }
         else {
           // still might be able to glean disambiguation
           // from the ambiguous hits at this time
           for (size_t a = 0; a < ambigwids.size(); a++)
             for (size_t w = 0; w < wids.size(); w++)
               if (ambigwids[a].TPC == wids[w].TPC &&
                   geom->WireIDsIntersect(ambigwids[a], wids[w], widIntersect))
                 widAcounts[a]++;
         }
       } // end loop through close-time hits
 
       // For now, just make a hit if either ambig or disambig hits
       // unanimously intersect a single wireID
       unsigned int Dcount(0), Acount(0);
       for (size_t d = 0; d < widDcounts.size(); d++)
         Dcount += widDcounts[d];
       for (size_t a = 0; a < widAcounts.size(); a++)
         Acount += widAcounts[a];
       for (size_t d = 0; d < widDcounts.size(); d++) {
         if (Dcount == widDcounts[d] && Dcount > 0 && Acount == 0) {
           this->MakeDisambigHit(ambighitPtr, ambigwids[d], apa);
           nDisambiguations++;
         }
       }
     } // end loop through still ambiguous hits
 
     return nDisambiguations;
   }

void apa::DisambigAlg::Crawl ( unsigned int apa )

Extend what we disambiguation we do have in apa.

\ todo: Evaluate how aggressive we can be here. How far should we jump? In what cases should we quit out?

Definition at line 370 of file DisambigAlg.cxx.

References fAPAToUVHits, fChanTimeToWid, fHasBeenDisambiged, fNChanJumps, hits(), and MakeCloseHits().

Referenced by RunDisambig(), and UseEndPts().

   {
 
     std::vector<art::Ptr<recob::Hit>> hits = fAPAToUVHits[apa];
 
     // repeat this method until stable
     unsigned int nExtended(1);
     while (nExtended > 0) {
       nExtended = 0;
 
       // Look for any disambiguated hit ...
       for (size_t h = 0; h < hits.size(); h++) {
         std::pair<double, double> ChanTime(hits[h]->Channel() * 1., hits[h]->PeakTime() * 1.);
         if (!fHasBeenDisambiged[apa][ChanTime]) continue;
         double stD = hits[h]->PeakTimePlusRMS(-1.);
         double etD = hits[h]->PeakTimePlusRMS(+1.);
         double hitWindow = etD - stD;
         geo::WireID Dwid = fChanTimeToWid[ChanTime];
 
         // ... and if any neighboring-channel hits are close enough in time,
         // extend the disambiguation to the neighboring wire.
         unsigned int extensions = 0;
         for (unsigned int ext = 1; ext < fNChanJumps + 1; ext++) {
           unsigned int N(0);
           double timeExt = hitWindow * ext;
           N += this->MakeCloseHits((int)(-ext), Dwid, stD - 5 - timeExt, etD + 5 + timeExt);
           N += this->MakeCloseHits((int)(ext), Dwid, stD - 5 - timeExt, etD + 5 + timeExt);
           extensions += N;
         }
         nExtended += extensions;
 
       } // end UV hits loop
 
     } // end while still disambiguating***
 
     // *** nested while loops allow disambigauation to hop the channel-wrap-boundary
   }

unsigned int apa::DisambigAlg::FindChanTimeEndPts	(	detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	apa
	)

Basic endpoint-hit finder per apa.

\ todo: Clean up and break down into two functions. \ todo: Make the conditions more robust to some spotty hits around a potential endpoint.

Definition at line 411 of file DisambigAlg.cxx.

References recob::Hit::Channel(), apa::APAGeometryAlg::ChannelsInView(), detinfo::DetectorPropertiesData::ConvertTicksToX(), geo::CryostatID::Cryostat, fAPAGeo, fAPAToEndPHits, fAPAToHits, fCloseHitsRadius, apa::APAGeometryAlg::FirstChannelInView(), fMaxEndPDegRange, geom, geo::kV, geo::kZ, recob::Hit::PeakTime(), util::pi(), util::size(), geo::TPCID::TPC, recob::Hit::View(), recob::Hit::WireID(), and geo::GeometryCore::WirePitch().

Referenced by RunDisambig().

   {
 
     double pi = 3.14159265;
     double fMaxEndPRadRange = fMaxEndPDegRange / 180. * (2 * pi);
 
     for (size_t h = 0; h < fAPAToHits[apa].size(); h++) {
       art::Ptr<recob::Hit> centhit = fAPAToHits[apa][h];
       geo::View_t view = centhit->View();
       unsigned int plane = 0;
       if (view == geo::kV) { plane = 1; }
       else if (view == geo::kZ)
         plane = 2;
       std::vector<double> ChanTimeCenter(2, 0.);
       unsigned int relchan = centhit->Channel() - fAPAGeo.FirstChannelInView(centhit->Channel());
       ChanTimeCenter[0] = relchan * geom->WirePitch(view);
       ChanTimeCenter[1] = detProp.ConvertTicksToX(centhit->PeakTime(),
                                                   plane,
                                                   apa * 2, // tpc doesnt matter
                                                   centhit->WireID().Cryostat);
       //std::vector< art::Ptr<recob::Hit> > CloseHits;
       std::vector<std::vector<double>> CloseHitsChanTime;
       std::vector<double> FurthestCloseChanTime(2, 0.); //double maxDist = 0.;
       std::vector<double> ClosestChanTime(2, 0.);
       double minDist = fCloseHitsRadius + 1.;
       double ChanDistRange = fAPAGeo.ChannelsInView(view) * geom->WirePitch(view);
 
       for (size_t c = 0; c < fAPAToHits[apa].size(); c++) {
         art::Ptr<recob::Hit> closehit = fAPAToHits[apa][c];
         if (view != closehit->View()) continue;
         if (view == geo::kZ && centhit->WireID().TPC != closehit->WireID().TPC) continue;
         unsigned int plane = 0;
         if (view == geo::kV) { plane = 1; }
         else if (view == geo::kZ)
           plane = 2;
         std::vector<double> ChanTimeClose(2, 0.);
         unsigned int relchanclose =
           closehit->Channel() - fAPAGeo.FirstChannelInView(closehit->Channel());
         ChanTimeClose[0] = relchanclose * geom->WirePitch(view);
         ChanTimeClose[1] = detProp.ConvertTicksToX(closehit->PeakTime(),
                                                    plane,
                                                    apa * 2, // tpc doesnt matter
                                                    closehit->WireID().Cryostat);
         if (ChanTimeClose == ChanTimeCenter) continue; // move on if the same one
 
         double ChanDist = ChanTimeClose[0] - ChanTimeCenter[0];
         if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
 
         double distance = std::hypot(ChanDist, ChanTimeClose[1] - ChanTimeCenter[1]);
 
         if (distance <= fCloseHitsRadius) CloseHitsChanTime.push_back(ChanTimeClose);
 
         if (distance < minDist) {
           ClosestChanTime = ChanTimeClose;
           minDist = distance;
         }
 
       } // end close-by hit loop
 
       if (CloseHitsChanTime.size() < 5) continue; // quick fix, to-be improved
 
       double minRad(2 * pi + 1.), maxRad(0.);
       bool CloseToNegPi(false), CloseToPosPi(false);
       for (size_t i = 0; i < CloseHitsChanTime.size(); i++) {
         std::vector<double> ThisChanTime(CloseHitsChanTime[i]);
         //if(ThisChanTime==ClosestChanTime) continue;
         double ChanDist = ThisChanTime[0] - ChanTimeCenter[0];
         if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
         double hitrad = std::atan2(ThisChanTime[1] - ChanTimeCenter[1], ChanDist);
         if (hitrad > maxRad) maxRad = hitrad;
         if (hitrad < minRad) minRad = hitrad;
         if (hitrad + fMaxEndPRadRange > pi)
           CloseToPosPi = true;
         else if (hitrad - fMaxEndPRadRange < -pi)
           CloseToNegPi = true;
       }
 
       // activity at this boundary automatically kills the test, move boundary and redo
       if (CloseToPosPi && CloseToNegPi) {
         for (size_t i = 0; i < CloseHitsChanTime.size(); i++) {
           std::vector<double> ThisChanTime(CloseHitsChanTime[i]);
           //if(ThisChanTime==ClosestChanTime) continue;
           double ChanDist = ThisChanTime[0] - ChanTimeCenter[0];
           if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
           double hitrad = std::atan2(ThisChanTime[1] - ChanTimeCenter[1], ChanDist);
           if (hitrad > 0) hitrad = pi - hitrad; // reflec across pi/2 line
           if (hitrad < 0) hitrad = -pi - hitrad;
           if (hitrad > maxRad) maxRad = hitrad;
           if (hitrad < minRad) minRad = hitrad;
         }
       }
 
       if (maxRad - minRad < fMaxEndPRadRange) fAPAToEndPHits[apa].push_back(centhit);
 
     } // end UV hit loop
 
     if (fAPAToEndPHits[apa].size() == 0) return 0;
     mf::LogVerbatim("FindChanTimeEndPts") << "          Found " << fAPAToEndPHits[apa].size()
                                           << " endpoint hits in apa " << apa << std::endl;
     for (size_t ep = 0; ep < fAPAToEndPHits[apa].size(); ep++) {
       art::Ptr<recob::Hit> epHit = fAPAToEndPHits[apa][ep];
       mf::LogVerbatim("FindChanTimeEndPts") << "           endP on channel " << epHit->Channel()
                                             << " at time " << epHit->PeakTime() << std::endl;
     }
 
     return fAPAToEndPHits[apa].size();
   }

bool apa::DisambigAlg::HitsOverlapInTime	(	detinfo::DetectorPropertiesData const &	detProp,
		recob::Hit const &	hitA,
		recob::Hit const &	hitB
	)

private

Definition at line 193 of file DisambigAlg.cxx.

References geo::kU, geo::kV, geo::kZ, recob::Hit::PeakTimeMinusRMS(), recob::Hit::PeakTimePlusRMS(), detinfo::DetectorPropertiesData::TimeOffsetU(), detinfo::DetectorPropertiesData::TimeOffsetV(), detinfo::DetectorPropertiesData::TimeOffsetZ(), and recob::Hit::View().

Referenced by CompareViews(), TrivialDisambig(), and UseEndPts().

   {
     double AsT = hitA.PeakTimeMinusRMS();
     double AeT = hitA.PeakTimePlusRMS();
     double BsT = hitB.PeakTimeMinusRMS();
     double BeT = hitB.PeakTimePlusRMS();
 
     if (hitA.View() == geo::kU) {
       AsT -= detProp.TimeOffsetU();
       AeT -= detProp.TimeOffsetU();
     }
     else if (hitA.View() == geo::kV) {
       AsT -= detProp.TimeOffsetV();
       AeT -= detProp.TimeOffsetV();
     }
     else if (hitA.View() == geo::kZ) {
       AsT -= detProp.TimeOffsetZ();
       AeT -= detProp.TimeOffsetZ();
     }
 
     if (hitB.View() == geo::kU) {
       BsT += detProp.TimeOffsetU();
       BeT -= detProp.TimeOffsetU();
     }
     else if (hitB.View() == geo::kV) {
       BsT -= detProp.TimeOffsetV();
       BeT -= detProp.TimeOffsetV();
     }
     else if (hitA.View() == geo::kZ) { // FIXME: Shouldn't this be hitB, BsT, and BeT?
       AsT -= detProp.TimeOffsetZ();
       AeT -= detProp.TimeOffsetZ();
     }
 
     return (AsT <= BsT && BsT <= AeT) || (AsT <= BeT && BeT <= AeT) || (BsT <= AsT && AsT <= BeT) ||
            (BsT <= AeT && AeT <= BeT);
   }

bool apa::DisambigAlg::HitsReasonablyMatch	(	art::Ptr< recob::Hit >	hitA,
		art::Ptr< recob::Hit >	hitB
	)

private

unsigned int apa::DisambigAlg::MakeCloseHits	(	int	ext,
		geo::WireID	wid,
		double	Dmin,
		double	Dmax
	)

private

Having disambiguated a time range on a wireID, extend to neighboring channels.

Definition at line 308 of file DisambigAlg.cxx.

References recob::Hit::Channel(), apa::APAGeometryAlg::ChannelsInView(), apa::APAGeometryAlg::ChannelToAPA(), geo::GeometryCore::ChannelToWire(), fAPAGeo, fChannelToHits, fHasBeenDisambiged, apa::APAGeometryAlg::FirstChannelInView(), geom, geo::kZ, MakeDisambigHit(), recob::Hit::PeakTime(), recob::Hit::PeakTimeMinusRMS(), recob::Hit::PeakTimePlusRMS(), geo::GeometryCore::PlaneWireToChannel(), geo::TPCID::TPC, geo::GeometryCore::View(), w, and geo::WireID::Wire.

Referenced by Crawl().

   {
     // Function to look, on a channel *ext* channels away from a
     // disambiguated hit channel, for hits with time windows touching
     // range *Dmin to Dmax*. If found, make such a hit to have a
     // wireID adjacent to supplied *wid*.  Returns number of NEW hits
     // made.
 
     raw::ChannelID_t Dchan = geom->PlaneWireToChannel(Dwid);
     geo::View_t view = geom->View(Dchan);
     if (view == geo::kZ)
       throw cet::exception("MakeCloseHits") << "Function not meant for non-wrapped channels.\n";
 
     // Account for wrapping
     raw::ChannelID_t firstChan = fAPAGeo.FirstChannelInView(view, Dchan);
     unsigned int ChanPerView = fAPAGeo.ChannelsInView(view);
     int tempchan = Dchan + ext; // need sign for the set of channels starting with channel 0
     if (tempchan < (int)firstChan) tempchan += ChanPerView;
     if (tempchan > (int)(firstChan + ChanPerView - 1)) tempchan -= ChanPerView;
     raw::ChannelID_t chan = (raw::ChannelID_t)(tempchan);
 
     // There may just be no hits
     if (fChannelToHits.count(chan) == 0) return 0;
 
     // There are close channel hits, so for each
     unsigned int apa(0), cryo(0);
     fAPAGeo.ChannelToAPA(chan, apa, cryo);
     unsigned int MakeCount(0);
     for (size_t i = 0; i < fChannelToHits[chan].size(); i++) {
       art::Ptr<recob::Hit> closeHit = fChannelToHits[chan][i];
       double st = closeHit->PeakTimeMinusRMS();
       double et = closeHit->PeakTimePlusRMS();
       std::vector<geo::WireID> wids = geom->ChannelToWire(chan);
 
       if (!(Dmin <= st && st <= Dmax) && !(Dmin <= et && et <= Dmax)) continue;
 
       // Found hit with window overlapping given range,
       // now find the only reasonable wireID.
       for (size_t w = 0; w < wids.size(); w++) {
         if (wids[w].TPC != Dwid.TPC) continue;
         if ((int)(wids[w].Wire) - (int)(Dwid.Wire) != ext) continue;
 
         // In this case, we have a unique wireID.
         // Check to see if it has already been made - if so, do not incriment count
         std::pair<double, double> ChanTime(closeHit->Channel() * 1., closeHit->PeakTime() * 1.);
         if (!fHasBeenDisambiged[apa][ChanTime]) {
           this->MakeDisambigHit(closeHit, wids[w], apa);
           MakeCount++;
           //std::cout << "     Close hit found on channel " << chan << ", time " << st<<"-"<<et << "... \n";
           //std::cout << " ... giving it wireID ("<< Dwid.Cryostat <<"," << Dwid.TPC
           //      << "," << Dwid.Plane << "," << Dwid.Wire << ")\n";
         }
         break;
       } // end find right wireID
 
     } // end loop through all hits on chan
 
     return MakeCount;
   }

void apa::DisambigAlg::MakeDisambigHit	(	art::Ptr< recob::Hit > const &	hit,
		geo::WireID	wid,
		unsigned int	apa
	)

private

Makes a disambiguated hit while keeping track of what has already been disambiguated.

Definition at line 174 of file DisambigAlg.cxx.

References recob::Hit::Channel(), fAPAToDHits, fChanTimeToWid, fHasBeenDisambiged, geo::CryostatID::isValid, and recob::Hit::PeakTime().

Referenced by CompareViews(), MakeCloseHits(), TrivialDisambig(), and UseEndPts().

   {
     std::pair<double, double> ChanTime(hit->Channel() * 1., hit->PeakTime() * 1.);
     if (fHasBeenDisambiged[apa][ChanTime]) return;
 
     if (!wid.isValid) {
       mf::LogWarning("InvalidWireID") << "wid is invalid, hit not being made\n";
       return;
     }
 
     fAPAToDHits[apa].emplace_back(hit, wid);
     fHasBeenDisambiged[apa][ChanTime] = true;
     fChanTimeToWid[ChanTime] = wid;
   }

void apa::DisambigAlg::RunDisambig	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		art::Handle< std::vector< recob::Hit >>	GausHits
	)

Definition at line 53 of file DisambigAlg.cxx.

References AssessDisambigSoFar(), bt_serv, recob::Hit::Channel(), apa::APAGeometryAlg::ChannelToAPA(), CompareViews(), Crawl(), fAPAGeo, fAPAToDHits, fAPAToEndPHits, fAPAToHits, fAPAToUVHits, fAPAToZHits, fChannelToHits, fChanTimeToWid, fCompareViews, fCrawl, fDisambigHits, fHasBeenDisambiged, art::fill_ptr_vector(), FindChanTimeEndPts(), fnDUSoFar, fnDVSoFar, fnUSoFar, fnVSoFar, fUeffSoFar, fUseEndP, fVeffSoFar, cheat::BackTrackerService::HitToXYZ(), geo::kU, geo::kV, geo::kZ, recob::Hit::PeakTime(), TrivialDisambig(), UseEndPts(), and recob::Hit::View().

Referenced by apa::APAHitFinder::produce().

   {
     fUeffSoFar.clear();
     fVeffSoFar.clear();
     fnUSoFar.clear();
     fnVSoFar.clear();
     fnDUSoFar.clear();
     fnDVSoFar.clear();
     fChannelToHits.clear();
     fAPAToUVHits.clear();
     fAPAToZHits.clear();
     fAPAToHits.clear();
     fAPAToEndPHits.clear();
     fAPAToDHits.clear();
     fDisambigHits.clear();
     fChanTimeToWid.clear();
     fHasBeenDisambiged.clear();
 
     std::vector<art::Ptr<recob::Hit>> ChHits;
     art::fill_ptr_vector(ChHits, ChannelHits);
 
     fHasBeenDisambiged.clear();
     unsigned int skipNoise(0);
     // Map hits by channel/APA, initialize the disambiguation status map
     for (size_t h = 0; h < ChHits.size(); h++) {
       art::Ptr<recob::Hit> const& hit = ChHits[h];
 
       // **temporary** option to skip noise hits
       try {
         bt_serv->HitToXYZ(clockData, hit);
       }
       catch (...) {
         skipNoise++;
         continue;
       }
 
       geo::View_t view = hit->View();
       unsigned int apa(0), cryo(0);
       fAPAGeo.ChannelToAPA(hit->Channel(), apa, cryo);
       fAPAToHits[apa].push_back(hit);
       if (view == geo::kZ) {
         fAPAToZHits[apa].push_back(hit);
         continue;
       }
       else if (view == geo::kU || view == geo::kV) {
         std::pair<double, double> ChanTime(hit->Channel() * 1., hit->PeakTime() * 1.);
         this->fHasBeenDisambiged[apa][ChanTime] = false;
         fChannelToHits[hit->Channel()].push_back(hit);
         fAPAToUVHits[apa].push_back(hit);
       }
     }
 
     if (skipNoise > 0)
       mf::LogWarning("DisambigAlg")
         << "\nSkipped " << skipNoise << " induction noise hits using the BackTrackerService.\n"
         << "This is only to temporarily deal with the excessive amount of noise due to the bad "
            "deconvolution.\n";
 
     mf::LogVerbatim("RunDisambig") << "\n~~~~~~~~~~~ Running Disambiguation ~~~~~~~~~~~\n";
 
     std::map<unsigned int, std::vector<art::Ptr<recob::Hit>>>::iterator APA_it;
     for (APA_it = fAPAToUVHits.begin(); APA_it != fAPAToUVHits.end(); APA_it++) {
       unsigned int apa = APA_it->first;
 
       mf::LogVerbatim("RunDisambig") << "APA " << apa << ":";
 
       fUeffSoFar[apa] = 0.;
       fVeffSoFar[apa] = 0.;
       fnUSoFar[apa] = 0;
       fnVSoFar[apa] = 0;
       fnDUSoFar[apa] = 0;
       fnDVSoFar[apa] = 0;
 
       // Always run this...
       this->TrivialDisambig(clockData, detProp, apa);
       this->AssessDisambigSoFar(apa);
       mf::LogVerbatim("RunDisambig")
         << "  Trivial Disambig -->  " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U,  "
         << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
 
       // ... and pick the rest with the configurations.
       if (fCrawl) {
         this->Crawl(apa);
         this->AssessDisambigSoFar(apa);
         mf::LogVerbatim("RunDisambig")
           << "  Crawl            -->  " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U,  "
           << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
       }
 
       if (fUseEndP) {
         this->FindChanTimeEndPts(detProp, apa);
         this->UseEndPts(detProp, apa); // does the crawl from inside
         this->AssessDisambigSoFar(apa);
         mf::LogVerbatim("RunDisambig")
           << "  Endpoint Crawl   -->  " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U,  "
           << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
       }
 
       if (fCompareViews) {
         unsigned int nDisambig(1);
         while (nDisambig > 0) {
           nDisambig = this->CompareViews(detProp, apa);
           this->Crawl(apa);
         }
         this->AssessDisambigSoFar(apa);
         mf::LogVerbatim("RunDisambig")
           << "  Compare Views    -->  " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U,  "
           << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
       }
 
       // For now just buld a simple list to get from the module
       for (size_t i = 0; i < fAPAToDHits[apa].size(); i++)
         fDisambigHits.push_back(fAPAToDHits[apa][i]);
 
     } // end loop through APA
   }

void apa::DisambigAlg::TrivialDisambig	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	apa
	)

Make the easiest and safest disambiguations in apa.

\ todo: Figure out why sometimes non-noise hits dont match any Z hits at all.

\ todo: Add mechanism to at least eliminate the wids that aren't even possible, for the benefit of future methods

Definition at line 234 of file DisambigAlg.cxx.

References bt_serv, geo::GeometryCore::ChannelToWire(), geo::CryostatID::Cryostat, d, fAPAToUVHits, fAPAToZHits, geom, geo::TPCGeo::GetCenter(), HitsOverlapInTime(), cheat::BackTrackerService::HitToXYZ(), MakeDisambigHit(), geo::GeometryCore::NearestChannel(), geo::GeometryCore::NTPC(), lar::to_element, geo::TPCID::TPC, geo::GeometryCore::TPC(), w, and geo::GeometryCore::WireEndPoints().

Referenced by RunDisambig().

   {
     // Loop through ambiguous hits (U/V) in this APA
     for (auto const& hitPtr : fAPAToUVHits[apa]) {
       auto const& hit = *hitPtr;
       raw::ChannelID_t chan = hit.Channel();
       unsigned int peakT = hit.PeakTime();
 
       std::vector<geo::WireID> hitwids = geom->ChannelToWire(chan);
       std::vector<bool> IsReasonableWid(hitwids.size(), false);
       unsigned short nPossibleWids(0);
       for (size_t w = 0; w < hitwids.size(); w++) {
         geo::WireID wid = hitwids[w];
 
         double xyzStart[3] = {0.};
         double xyzEnd[3] = {0.};
         geom->WireEndPoints(wid, xyzStart, xyzEnd);
         unsigned int side(wid.TPC % 2), cryo(wid.Cryostat);
         double zminPos(xyzStart[2]), zmaxPos(xyzEnd[2]);
 
         // get appropriate x and y with tpc center
         unsigned int tpc =
           2 * apa + side - cryo * geom->NTPC(); // apa number does not reset per cryo
         auto const tpcCenter = geom->TPC(geo::TPCID(cryo, tpc)).GetCenter();
 
         // get channel range
         auto Min = tpcCenter;
         Min.SetZ(zminPos);
         auto Max = tpcCenter;
         Max.SetZ(zmaxPos);
         raw::ChannelID_t ZminChan = geom->NearestChannel(Min, geo::PlaneID{cryo, tpc, 2});
         raw::ChannelID_t ZmaxChan = geom->NearestChannel(Max, geo::PlaneID{cryo, tpc, 2});
 
         for (auto const& zhit : fAPAToZHits[apa] | transform(to_element)) {
           raw::ChannelID_t chan = zhit.Channel();
           if (chan <= ZminChan || ZmaxChan <= chan) continue;
 
           if (this->HitsOverlapInTime(detProp, hit, zhit)) {
             IsReasonableWid[w] = true;
             nPossibleWids++;
             break;
           }
         }
 
       } // end hit chan-wid loop
 
       if (nPossibleWids == 0) {
         std::vector<double> xyz;
         try {
           xyz = bt_serv->HitToXYZ(clockData, hit);
         } // TEMPORARY
         catch (...) {
           continue;
         }
         mf::LogWarning("UniqueTimeSeg")
           << "U/V hit inconsistent with Z info; peak time is " << peakT << " in APA " << apa
           << " on channel " << hit.Channel();
       }
       else if (nPossibleWids == 1) {
         for (size_t d = 0; d < hitwids.size(); d++)
           if (IsReasonableWid[d]) this->MakeDisambigHit(hitPtr, hitwids[d], apa);
       }
       else if (nPossibleWids == 2) {
       }
 
     } // end ambig hits loop
   }

void apa::DisambigAlg::UseEndPts	(	detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	apa
	)

Try to associate endpoint hits and crawl from there

\ todo: This function could be made much cleaner and more compact

Definition at line 524 of file DisambigAlg.cxx.

References apa::APAGeometryAlg::APAChannelsIntersect(), recob::Hit::Channel(), Crawl(), fAPAGeo, fAPAToEndPHits, geom, geo::TPCGeo::GetCenter(), HitsOverlapInTime(), geo::kU, geo::kV, geo::kZ, MakeDisambigHit(), apa::APAGeometryAlg::NearestWireIDOnChan(), util::size(), apa::APAGeometryAlg::ThreeChanPos(), and geo::GeometryCore::TPC().

Referenced by RunDisambig().

   {
 
 
     if (fAPAToEndPHits[apa].size() == 0) {
       mf::LogVerbatim("UseEndPts") << "          APA " << apa << " has no endpoints.";
       return;
     }
     std::vector<art::Ptr<recob::Hit>> const& endPts = fAPAToEndPHits[apa];
 
     std::vector<std::vector<art::Ptr<recob::Hit>>> EndPMatch;
     unsigned short nZendPts(0);
 
     auto on_z_plane = [](art::Ptr<recob::Hit> const& hit) { return hit->View() == geo::kZ; };
     auto not_on_z_plane = [](art::Ptr<recob::Hit> const& hit) { return hit->View() != geo::kZ; };
     for (auto const& ZHitPtr : endPts | filter(on_z_plane)) {
       auto const& ZHit = *ZHitPtr;
       art::Ptr<recob::Hit> Uhit = ZHitPtr;
       art::Ptr<recob::Hit> Vhit = ZHitPtr;
       unsigned short Umatch(0), Vmatch(0);
       ++nZendPts;
 
       // look for U and V hits overlapping in time
       for (auto const& hitPtr : endPts | filter(not_on_z_plane)) {
         auto const& hit = *hitPtr;
         if (not HitsOverlapInTime(detProp, ZHit, hit)) continue;
 
         if (hit.View() == geo::kU) {
           Uhit = hitPtr;
           Umatch++;
         }
         else if (hit.View() == geo::kV) {
           Vhit = hitPtr;
           Vmatch++;
         }
       }
 
       unsigned int tpc(ZHit.WireID().TPC), cryo(ZHit.WireID().Cryostat);
       auto const tpcCenter = geom->TPC(geo::TPCID{cryo, tpc}).GetCenter();
 
       if (Umatch == 1 && Vmatch == 1) {
 
         std::vector<double> yzEndPt =
           fAPAGeo.ThreeChanPos(Uhit->Channel(), Vhit->Channel(), ZHit.Channel());
         double intersect[3] = {tpcCenter.X(), yzEndPt[0], yzEndPt[1]};
 
         geo::WireID Uwid = fAPAGeo.NearestWireIDOnChan(intersect, Uhit->Channel(), 0, tpc, cryo);
         geo::WireID Vwid = fAPAGeo.NearestWireIDOnChan(intersect, Vhit->Channel(), 1, tpc, cryo);
         this->MakeDisambigHit(Uhit, Uwid, apa);
         this->MakeDisambigHit(Vhit, Vwid, apa);
       }
       else if (Umatch == 1 && Vmatch != 1) {
 
         std::vector<geo::WireIDIntersection> widIntersects;
         fAPAGeo.APAChannelsIntersect(Uhit->Channel(), ZHit.Channel(), widIntersects);
         if (widIntersects.size() == 0)
           continue;
         else if (widIntersects.size() == 1) {
           double intersect[3] = {tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
           geo::WireID Uwid = fAPAGeo.NearestWireIDOnChan(intersect, Uhit->Channel(), 0, tpc, cryo);
           this->MakeDisambigHit(Uhit, Uwid, apa);
         }
         else {
           for (size_t i = 0; i < widIntersects.size(); i++) {
             // compare to V hit times, see if only one makes sense
           }
         }
       }
       else if (Umatch == 1 && Vmatch != 1) {
 
         std::vector<geo::WireIDIntersection> widIntersects;
         fAPAGeo.APAChannelsIntersect(Vhit->Channel(), ZHit.Channel(), widIntersects);
         if (widIntersects.size() == 0)
           continue;
         else if (widIntersects.size() == 1) {
           double intersect[3] = {tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
           geo::WireID Vwid = fAPAGeo.NearestWireIDOnChan(intersect, Vhit->Channel(), 0, tpc, cryo);
           this->MakeDisambigHit(Vhit, Vwid, apa);
         }
       }
     }
 
     if (nZendPts == 0 && endPts.size() == 2 &&
         this->HitsOverlapInTime(detProp, *endPts[0], *endPts[1])) {
       std::vector<geo::WireIDIntersection> widIntersects;
       fAPAGeo.APAChannelsIntersect(endPts[0]->Channel(), endPts[1]->Channel(), widIntersects);
       if (widIntersects.size() == 1) {
         unsigned int cryo = endPts[0]->WireID().Cryostat;
         unsigned int tpc = widIntersects[0].TPC;
         auto const tpcCenter = geom->TPC(geo::TPCID(cryo, tpc)).GetCenter();
         double intersect[3] = {tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
         unsigned int plane0(0), plane1(0);
         if (endPts[0]->View() == geo::kV) plane0 = 1;
         if (endPts[1]->View() == geo::kV) plane1 = 1;
         geo::WireID wid0 =
           fAPAGeo.NearestWireIDOnChan(intersect, endPts[0]->Channel(), plane0, tpc, cryo);
         this->MakeDisambigHit(endPts[0], wid0, apa);
         geo::WireID wid1 =
           fAPAGeo.NearestWireIDOnChan(intersect, endPts[1]->Channel(), plane1, tpc, cryo);
         this->MakeDisambigHit(endPts[1], wid1, apa);
       }
     }
 
     this->Crawl(apa);
   }