apa::DisambigAlg Class Reference

#include "DisambigAlg.h"

Public Member Functions
	DisambigAlg (fhicl::ParameterSet const &pset)
virtual	~DisambigAlg ()
void	reconfigure (fhicl::ParameterSet const &p)
void	RunDisambig (art::Handle< std::vector< recob::Hit > > GausHits)
	Run disambiguation as currently configured. More...
void	TrivialDisambig (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 (unsigned int apa)
	Basic endpoint-hit finder per apa. More...
void	UseEndPts (unsigned int apa)
	Try to associate endpoint hits and crawl from there. More...
unsigned int	CompareViews (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 > 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 (art::Ptr< recob::Hit > hitA, art::Ptr< recob::Hit > hitB)
bool	HitsReasonablyMatch (art::Ptr< recob::Hit > hitA, art::Ptr< recob::Hit > hitB)

Private Attributes
apa::APAGeometryAlg	fAPAGeo
art::ServiceHandle< geo::Geometry >	geom
const detinfo::DetectorProperties *	detprop
art::ServiceHandle< cheat::BackTrackerService >	bt_serv
	For TEMPORARY monitering of potential problems. More...
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
	Distance (cm) away from a hit to look when checking if it's an endpoint. More...
double	fMaxEndPDegRange

Detailed Description

Definition at line 40 of file DisambigAlg.h.

Constructor & Destructor Documentation

apa::DisambigAlg::DisambigAlg

(

fhicl::ParameterSet const &

pset

)

Definition at line 45 of file DisambigAlg.cxx.

References reconfigure().

  : fAPAGeo(pset.get< fhicl::ParameterSet >("APAGeometryAlg"))
{
  this->reconfigure(pset); 
}

apa::DisambigAlg::~DisambigAlg ( )

virtual

Definition at line 52 of file DisambigAlg.cxx.

{
}

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 676 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

(

unsigned int

apa

)

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

Definition at line 706 of file DisambigAlg.cxx.

References recob::Hit::Channel(), geo::GeometryCore::ChannelToWire(), d, fAPAToUVHits, fChanTimeToWid, fHasBeenDisambiged, geom, HitsOverlapInTime(), MakeDisambigHit(), recob::Hit::PeakTime(), recob::Hit::View(), w, and geo::GeometryCore::WireIDsIntersect().

Referenced by RunDisambig().

{

  unsigned int nDisambiguations(0);

  // loop through all hits that are still ambiguous
  for(size_t h=0; h < fAPAToUVHits[apa].size(); h++){
    art::Ptr<recob::Hit>      ambighit  = fAPAToUVHits[apa][h];
    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(size_t i=0; i < fAPAToUVHits[apa].size(); i++){
      art::Ptr<recob::Hit> hit = fAPAToUVHits[apa][i];
      if(hit->View()==view || !this->HitsOverlapInTime(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(ambighit, ambigwids[d], apa);
        nDisambiguations++;  } }
    for(size_t a=0; a<widAcounts.size(); a++){
      if( Acount == widAcounts[a] && Acount==1 ){
        //this->MakeDisambigHit(ambighit, ambigwids[a], 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 399 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

(

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 444 of file DisambigAlg.cxx.

References recob::Hit::Channel(), apa::APAGeometryAlg::ChannelsInView(), detinfo::DetectorProperties::ConvertTicksToX(), geo::CryostatID::Cryostat, detprop, fAPAGeo, fAPAToEndPHits, fAPAToHits, fCloseHitsRadius, apa::APAGeometryAlg::FirstChannelInView(), fMaxEndPDegRange, geom, geo::kV, geo::kZ, recob::Hit::PeakTime(), util::pi(), 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);

  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();

  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::sqrt(   std::pow(ChanDist,2) 
                                   + std::pow(ChanTimeClose[1]-ChanTimeCenter[1],2)  );

      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 ( art::Ptr< recob::Hit >  hitA, art::Ptr< recob::Hit >  hitB  )

private

Definition at line 225 of file DisambigAlg.cxx.

References detprop, geo::kU, geo::kV, geo::kZ, recob::Hit::PeakTimeMinusRMS(), recob::Hit::PeakTimePlusRMS(), detinfo::DetectorProperties::TimeOffsetU(), detinfo::DetectorProperties::TimeOffsetV(), detinfo::DetectorProperties::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();

    const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
    
    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 ){ 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 326 of file DisambigAlg.cxx.

References recob::Hit::Channel(), apa::APAGeometryAlg::ChannelsInView(), apa::APAGeometryAlg::ChannelToAPA(), geo::GeometryCore::ChannelToWire(), geo::CryostatID::Cryostat, fAPAGeo, fChannelToHits, fHasBeenDisambiged, apa::APAGeometryAlg::FirstChannelInView(), geom, geo::kZ, MakeDisambigHit(), recob::Hit::PeakTime(), recob::Hit::PeakTimeMinusRMS(), recob::Hit::PeakTimePlusRMS(), geo::PlaneID::Plane, 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.Plane, Dwid.Wire, Dwid.TPC, Dwid.Cryostat);
  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 >  hit, geo::WireID  wid, unsigned int  apa  )

private

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

Definition at line 201 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; }
    
  std::pair<art::Ptr<recob::Hit>,geo::WireID> Dhit(hit, wid);
  fAPAToDHits[apa].push_back(Dhit);
  fHasBeenDisambiged[apa][ChanTime] = true;
  fChanTimeToWid[ChanTime] = wid;
  return;
  
}

void apa::DisambigAlg::reconfigure

(

fhicl::ParameterSet const &

p

)

Definition at line 57 of file DisambigAlg.cxx.

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

Referenced by DisambigAlg().

{

  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");

}

void apa::DisambigAlg::RunDisambig

(

art::Handle< std::vector< recob::Hit > >

GausHits

)

Run disambiguation as currently configured.

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

{

  // **tomporarily** here to look at performance without noise hits
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;

  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> hit = ChHits[h];

    // **temporary** option to skip noise hits
    try{ bt_serv->HitToXYZ(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(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(apa);
      this->UseEndPts(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 = 0;
        nDisambig = this->CompareViews(apa);
        this->Crawl(apa);
      }
      this->AssessDisambigSoFar(apa);
      mf::LogVerbatim("RunDisambig") << "  Compare Views    -->  " 
                                     << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U,  "
                                     << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
    }      


    //this->GatherLeftoverHits()

    // 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

(

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 250 of file DisambigAlg.cxx.

References bt_serv, recob::Hit::Channel(), geo::GeometryCore::ChannelToWire(), geo::CryostatID::Cryostat, geo::GeometryCore::Cryostat(), d, fAPAToUVHits, fAPAToZHits, geom, HitsOverlapInTime(), cheat::BackTrackerService::HitToXYZ(), geo::TPCGeo::LocalToWorld(), MakeDisambigHit(), geo::GeometryCore::NearestChannel(), geo::GeometryCore::NTPC(), recob::Hit::PeakTime(), geo::PlaneID::Plane, geo::CryostatGeo::TPC(), geo::TPCID::TPC, w, geo::WireID::Wire, geo::GeometryCore::WireEndPoints(), and z.

Referenced by RunDisambig().

{

  // Loop through ambiguous hits (U/V) in this APA
  for( size_t h=0; h<fAPAToUVHits[apa].size(); h++ ){
    const art::Ptr<recob::Hit> hit = fAPAToUVHits[apa][h];
    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.Cryostat, wid.TPC, wid.Plane, wid.Wire, 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
      TVector3 tpcCenter(0,0,0);
      unsigned int tpc = 2*apa + side - cryo*geom->NTPC(); // apa number does not reset per cryo
      tpcCenter = geom->Cryostat(cryo).TPC(tpc).LocalToWorld(tpcCenter);
      
      // get channel range
      TVector3 Min(tpcCenter); Min[2] = zminPos;
      TVector3 Max(tpcCenter); Max[2] = zmaxPos;
      raw::ChannelID_t ZminChan = geom->NearestChannel( Min, 2, tpc, cryo );
      raw::ChannelID_t ZmaxChan = geom->NearestChannel( Max, 2, tpc, cryo );
      
      for( size_t z=0; z < fAPAToZHits[apa].size(); z++ ){
        raw::ChannelID_t chan = fAPAToZHits[apa][z]->Channel();
        if( chan <= ZminChan || ZmaxChan <= chan ) continue;
        art::Ptr<recob::Hit> zhit = fAPAToZHits[apa][z];

//      try{ bt_serv->HitToXYZ(zhit); }
//      catch(...){ 
//        mf::LogWarning("DisambigAlg")<<"skipping a noise collection hit"; 
//        continue;
//      }

        if( this->HitsOverlapInTime(hit,zhit) ){
          IsReasonableWid[w] = true; nPossibleWids++; break; }
      } 

    } // end hit chan-wid loop

      

    if(nPossibleWids==0){
      std::vector<double> xyz;
      try{ xyz = bt_serv->HitToXYZ(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( hit, hitwids[d], apa );
    }
    else if(nPossibleWids==2){
    }
    

  } // end ambig hits loop  


}

void apa::DisambigAlg::UseEndPts

(

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 557 of file DisambigAlg.cxx.

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

Referenced by RunDisambig().

{


  if(fAPAToEndPHits[apa].size()==0){
    mf::LogVerbatim("UseEndPts") << "          APA " << apa << " has no endpoints.";
    return; }
  std::vector< art::Ptr<recob::Hit> > endPts = fAPAToEndPHits[apa];


  std::vector<std::vector< art::Ptr<recob::Hit> > > EndPMatch;
  unsigned short nZendPts(0);
  for(size_t ep=0; ep<endPts.size(); ep++){
    if(endPts[ep]->View()!=geo::kZ) continue;
    art::Ptr<recob::Hit> Zhit = endPts[ep];
    art::Ptr<recob::Hit> Uhit(Zhit), Vhit(Zhit);
    unsigned short Umatch(0), Vmatch(0);
    nZendPts++;
    
    // look for U and V hits overlapping in time
    for(size_t p=0; p<endPts.size(); p++){
      if(endPts[p]->View()==geo::kU){
        if( this->HitsOverlapInTime(Zhit, endPts[p]) ){
          Uhit = endPts[p];
          //endPts.erase(endPts.begin()+p);
          Umatch++; }
      } else if(endPts[p]->View()==geo::kV){
        if( this->HitsOverlapInTime(Zhit, endPts[p]) ){
          Vhit = endPts[p];
          //endPts.erase(endPts.begin()+p);
          Vmatch++; }
      }
    }


    TVector3 tpcCenter(0,0,0);
    unsigned int tpc(endPts[ep]->WireID().TPC), cryo(endPts[ep]->WireID().Cryostat);  
    tpcCenter = geom->Cryostat(cryo).TPC(tpc).LocalToWorld(tpcCenter);


    if(Umatch == 1 && Vmatch == 1){

      std::vector<double> yzEndPt 
        = fAPAGeo.ThreeChanPos(Uhit->Channel(), Vhit->Channel(), endPts[ep]->Channel());
      double intersect[3] = {tpcCenter[0], 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[0], 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[0], widIntersects[0].y, widIntersects[0].z};
        geo::WireID Vwid = fAPAGeo.NearestWireIDOnChan(intersect, Vhit->Channel(), 0, tpc, cryo);
        this->MakeDisambigHit(Vhit, Vwid, apa);
      } else {
        for(size_t i=0; i < widIntersects.size(); i++){
          // compare to V hit times, see if only one makes sense
        }
      }

    }
    
    //endPts.erase(endPts.begin()+ep);
  }


  if( nZendPts==0 && endPts.size()==2 && this->HitsOverlapInTime(endPts[0],endPts[1]) ){
    std::vector< geo::WireIDIntersection > widIntersects;
    fAPAGeo.APAChannelsIntersect(endPts[0]->Channel(), endPts[1]->Channel(), widIntersects);
    if( widIntersects.size() == 1 ){
      TVector3 tpcCenter(0,0,0);
      unsigned int cryo = endPts[0]->WireID().Cryostat;
      unsigned int tpc = widIntersects[0].TPC;
      tpcCenter = geom->Cryostat(cryo).TPC(tpc).LocalToWorld(tpcCenter);
      double intersect[3] = {tpcCenter[0], 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);
    } else if( widIntersects.size() > 1 ){
        for(size_t i=0; i < widIntersects.size(); i++){
          // compare to V hit times, see if only one makes sense
        }
      }
  }

 
  this->Crawl(apa);
 
}

Member Data Documentation

art::ServiceHandle<cheat::BackTrackerService> apa::DisambigAlg::bt_serv

private

For TEMPORARY monitering of potential problems.

Definition at line 76 of file DisambigAlg.h.

Referenced by RunDisambig(), and TrivialDisambig().

const detinfo::DetectorProperties* apa::DisambigAlg::detprop

private

Definition at line 75 of file DisambigAlg.h.

Referenced by FindChanTimeEndPts(), and HitsOverlapInTime().

apa::APAGeometryAlg apa::DisambigAlg::fAPAGeo

private

Definition at line 73 of file DisambigAlg.h.

Referenced by FindChanTimeEndPts(), MakeCloseHits(), RunDisambig(), and UseEndPts().

std::map< unsigned int, std::vector< std::pair<art::Ptr<recob::Hit>, geo::WireID> > > apa::DisambigAlg::fAPAToDHits

private

Hold the disambiguations per APA.

Definition at line 84 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), MakeDisambigHit(), and RunDisambig().

std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > apa::DisambigAlg::fAPAToEndPHits

private

\ todo: Channel/APA to hits can be done in a unified way

Definition at line 83 of file DisambigAlg.h.

Referenced by FindChanTimeEndPts(), RunDisambig(), and UseEndPts().

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

private

Definition at line 81 of file DisambigAlg.h.

Referenced by FindChanTimeEndPts(), and RunDisambig().

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

private

Definition at line 80 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), CompareViews(), Crawl(), RunDisambig(), and TrivialDisambig().

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

private

Definition at line 80 of file DisambigAlg.h.

Referenced by RunDisambig(), and TrivialDisambig().

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

private

Definition at line 79 of file DisambigAlg.h.

Referenced by MakeCloseHits(), and RunDisambig().

std::map<std::pair<double,double>, geo::WireID> apa::DisambigAlg::fChanTimeToWid

private

If a hit is disambiguated, map its chan and peak time to the chosen wireID.

Definition at line 90 of file DisambigAlg.h.

Referenced by CompareViews(), Crawl(), MakeDisambigHit(), and RunDisambig().

double apa::DisambigAlg::fCloseHitsRadius

private

Distance (cm) away from a hit to look when checking if it's an endpoint.

Definition at line 115 of file DisambigAlg.h.

Referenced by FindChanTimeEndPts(), and reconfigure().

bool apa::DisambigAlg::fCompareViews

private

Definition at line 113 of file DisambigAlg.h.

Referenced by reconfigure(), and RunDisambig().

bool apa::DisambigAlg::fCrawl

private

\ todo: Write function that compares hits more detailedly

Definition at line 111 of file DisambigAlg.h.

Referenced by reconfigure(), and RunDisambig().

std::vector< std::pair<art::Ptr<recob::Hit>, geo::WireID> > apa::DisambigAlg::fDisambigHits

The final list of hits to pass back to be made.

Definition at line 67 of file DisambigAlg.h.

Referenced by apa::APAHitFinder::produce(), and RunDisambig().

std::map< unsigned int, std::map<std::pair<double,double>, bool> > apa::DisambigAlg::fHasBeenDisambiged

private

Convenient way to keep track of disambiguation so far.

Definition at line 92 of file DisambigAlg.h.

Referenced by CompareViews(), Crawl(), MakeCloseHits(), MakeDisambigHit(), and RunDisambig().

double apa::DisambigAlg::fMaxEndPDegRange

private

Within the close hits radius, how spread can the majority of the activity be around a possible endpoint

Definition at line 116 of file DisambigAlg.h.

Referenced by FindChanTimeEndPts(), and reconfigure().

unsigned int apa::DisambigAlg::fNChanJumps

private

Number of channels the crawl can jump over.

Definition at line 114 of file DisambigAlg.h.

Referenced by Crawl(), and reconfigure().

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

Definition at line 64 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), and RunDisambig().

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

Definition at line 65 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), and RunDisambig().

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

Definition at line 62 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), and RunDisambig().

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

Definition at line 63 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), and RunDisambig().

std::map<unsigned int, double> apa::DisambigAlg::fUeffSoFar

Definition at line 60 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), and RunDisambig().

bool apa::DisambigAlg::fUseEndP

private

Definition at line 112 of file DisambigAlg.h.

Referenced by reconfigure(), and RunDisambig().

std::map<unsigned int, double> apa::DisambigAlg::fVeffSoFar

Definition at line 61 of file DisambigAlg.h.

Referenced by AssessDisambigSoFar(), and RunDisambig().

art::ServiceHandle<geo::Geometry> apa::DisambigAlg::geom

private

Definition at line 74 of file DisambigAlg.h.

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

---

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

• /cvmfs/larsoft.opensciencegrid.org/products/larreco/v07_10_02/source/larreco/RecoAlg/DisambigAlg.h
• /cvmfs/larsoft.opensciencegrid.org/products/larreco/v07_10_02/source/larreco/RecoAlg/DisambigAlg.cxx