Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
dt::DiscreteTracker Class Reference

Job module performing "discrete" track reconstruction. More...

Inheritance diagram for dt::DiscreteTracker:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 DiscreteTracker (const fhicl::ParameterSet &pset)
 
 ~DiscreteTracker ()
 
void produce (art::Event &evt)
 
void reconfigure (const fhicl::ParameterSet &pset)
 
void beginJob ()
 
void endJob ()
 
std::string workerType () const
 
void doBeginJob ()
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< TconsumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< TmayConsumeView (InputTag const &tag)
 

Protected Member Functions

bool HighActivity (const art::Ptr< rb::Cluster > &clust) const
 Cut to prevent choking on large air shower events. More...
 
bool TryExtendChainOne (Cand &cur, const Cand &prev, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
 
void ExtendChain (const View &v, const dt::View::ChunkMap &chunk_map, Chain &chain, Direction dir) const
 
Chain RecoOneTrack (const View &v, const dt::View::ChunkMap &chunk_map, std::vector< rb::Track > *seeds_out, bool &ok) const
 
void SpliceChains (std::vector< Chain > &chains) const
 
bool SpliceIsSane (const Chain &a, const Chain &b) const
 
std::vector< rb::TrackRecoDiscrete (const art::PtrVector< rb::CellHit > &cells, std::vector< rb::Track > *seeds_out, std::vector< rb::Track > *unspliced_out, std::vector< rb::Track > *chains_out) const
 
std::vector< SegmentPossibleNextSegs (const View &v, const dt::View::ChunkMap &chunk_map, const Cand &cand, Direction dir, std::vector< Chunk > &chunks_out, bool shw, bool exactHits) const
 
double DotScoreHelper (double m1, double m2) const
 
double DotScore (const Cand &c1, const Cand &c2, bool allowFlip=false) const
 
bool TryExtendCandOne (Cand &cand, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, const Cand &prev_dir, bool shw, bool exactHits) const
 
bool TryExtendCand (Cand &cand, const Cand &old_cand, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
 
Cand BestSeed (const dt::View &v, bool &ok) const
 
std::vector< ChainCheckAndFixSparseChain (Chain &chain) const
 
rb::VertexFindVertex (std::vector< rb::Track > xTrks, std::vector< rb::Track > yTrks) const
 Toy vertex finder. More...
 
void MergeParasiteTracks (std::vector< rb::Track > &tracks) const
 
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Protected Attributes

std::string fClusterInput
 Module name to look for Clusters under. More...
 
bool fConsiderBadChans
 Passed to View. More...
 
bool fObeyPreselection
 Check rb::IsFiltered? More...
 
double fMinDotProduct
 Joined Cands that only have one chunk in common must have a dot-product better than this. More...
 
int fVerticalTrackMinCells
 The number of adjacent cells in a single plane, otherwise not fit into any track, required to declare them a vertical track. More...
 
bool fOnlyOneTrack
 Stop after the first track in each view. More...
 
TH1 * fAngle
 Opening angle of all cands considered for joining. More...
 
TH1 * fSinglesAngle
 For all cands with only one overlapping chunk. More...
 
FILE * fDebugFile
 File that track details are written to. More...
 

Detailed Description

Job module performing "discrete" track reconstruction.

Definition at line 42 of file DiscreteTracker_module.cc.

Member Typedef Documentation

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 25 of file Producer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 18 of file EDProducer.h.

Constructor & Destructor Documentation

dt::DiscreteTracker::DiscreteTracker ( const fhicl::ParameterSet pset)
explicit

Definition at line 151 of file DiscreteTracker_module.cc.

References reconfigure().

152  : EDProducer(pset)
153  , fDebugFile(0)
154  {
155  produces<std::vector<rb::Track> >();
156  produces<std::vector<rb::Track> >("seeds");
157  produces<std::vector<rb::Track> >("unspliced");
158  produces<std::vector<rb::Track> >("chains");
159  produces<std::vector<rb::Vertex> >();
160  // Tracks are associated with the slice they came from
161  produces<art::Assns<rb::Track, rb::Cluster> >();
162 
163  reconfigure(pset);
164  }
FILE * fDebugFile
File that track details are written to.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
void reconfigure(const fhicl::ParameterSet &pset)
dt::DiscreteTracker::~DiscreteTracker ( )

Definition at line 167 of file DiscreteTracker_module.cc.

168  {
169  }

Member Function Documentation

void dt::DiscreteTracker::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 186 of file DiscreteTracker_module.cc.

References fAngle, and fSinglesAngle.

187  {
189  fAngle = tfs->make<TH1F>("angle", ";#Delta#theta", 100, 0, 1);
190  fSinglesAngle = tfs->make<TH1F>("singles_angle", ";#Delta#theta",
191  100, 0, 1);
192  }
TH1 * fAngle
Opening angle of all cands considered for joining.
TH1 * fSinglesAngle
For all cands with only one overlapping chunk.
Cand dt::DiscreteTracker::BestSeed ( const dt::View v,
bool &  ok 
) const
protected

Definition at line 750 of file DiscreteTracker_module.cc.

References dt::Cand::AddChunk(), dt::Chunk::AllHits(), ana::assert(), make_associated_cosmic_defs::chunk, make_goodruns_defs::chunks, art::PtrVector< T >::empty(), dt::View::GeoView(), dt::Chunk::GetSegs(), it, dt::kDownstream, dt::View::MakeChunks(), dt::View::NextPlane(), dt::Cand::NSegs(), NDAPDHVSetting::plane, dt::Cand::TryAddSeg(), dt::Chunk::Up(), and dt::Cand::Up().

Referenced by RecoOneTrack().

751  {
752  ok = true;
753 
754  // TODO, we could iterate over these instead
755  // const std::set<unsigned int> planes = v.GetPlanes();
756  const dt::View::ChunkMap chunks_map = v.MakeChunks();
757 
758  Cand best(v.GeoView(), true);
759 
760  if(chunks_map.empty()){
761  ok = false;
762  return best;
763  }
764 
765  std::list<Cand> live;
766 
767  // Go through planes
768  for(dt::View::ChunkMap::const_iterator it = chunks_map.begin(); it != chunks_map.end(); ++it){
769  const int plane = it->first;
770 
771  const std::vector<Chunk>& chunks = it->second;
772 
773  std::list<Cand> newLive;
774 
775  for(std::list<Cand>::iterator live_it = live.begin(); live_it != live.end(); ++live_it){
776 
777  if(plane != v.NextPlane(live_it->LastChunk().Plane())){
778  if(live_it->NSegs() > best.NSegs()) best = *live_it;
779  continue;
780  }
781 
782  // TODO: make a less wasteful version of TryAddSeg? Try to keep just
783  // the one list.
784  Cand curCand = *live_it; // Take a copy
785 
786  // Were any chunks any use?
787  bool any = false;
788 
789  for(unsigned int chunkIdx = 0; chunkIdx < chunks.size(); ++chunkIdx){
790  const Chunk& chunk = chunks[chunkIdx];
791  // TODO: potentially real data could have seeds that never had two
792  // adjacent hits?
793  if(chunk.AllHits().empty()) continue;
794 
795  // TODO: have chunks given back in two lists?
796  if(chunk.Up() != curCand.Up()) continue;
797 
798  const std::pair<Segment, Segment>& segs = chunk.GetSegs();
799 
800  if(curCand.TryAddSeg(segs.first, kDownstream)){
801  curCand.AddChunk(chunk, kDownstream);
802  if(curCand.TryAddSeg(segs.second, kDownstream)){
803  // Good, keep going
804  newLive.insert(newLive.end(), curCand);
805  curCand = *live_it; // New one to play with
806  }
807  else{
808  // Failed second half, first might be useful
809  if(curCand.NSegs() > best.NSegs()) best = curCand;
810  curCand = *live_it;
811  }
812  // Got at least the first seg in
813  any = true;
814  }
815  } // end for chunkIdx
816 
817  if(!any){
818  // Got stuck, but it might still be the best
819  if(curCand.NSegs() > best.NSegs()) best = curCand;
820  // No need to make a new curCand anymore
821  }
822  } // end for live_it
823 
824  live.swap(newLive);
825  newLive.clear();
826 
827  // Seed new cands from this chunk
828  for(unsigned int chunkIdx = 0; chunkIdx < chunks.size(); ++chunkIdx){
829  const Chunk& chunk = chunks[chunkIdx];
830 
831  const std::pair<Segment, Segment>& segs = chunk.GetSegs();
832 
833  Cand newcand(v.GeoView(), chunk.Up());
834  newcand.AddChunk(chunk, kDownstream);
835  bool startok = newcand.TryAddSeg(segs.first, kDownstream);
836  assert(startok);
837  startok = startok && newcand.TryAddSeg(segs.second, kDownstream);
838  assert(startok);
839  live.insert(live.end(), newcand);
840 
841  // Only the second seg
842  Cand newcand2(v.GeoView(), chunk.Up());
843  newcand2.AddChunk(chunk, kDownstream);
844  startok = startok && newcand2.TryAddSeg(segs.second, kDownstream);
845  live.insert(live.end(), newcand2);
846 
847  assert(startok);
848  } // end for chunkIdx
849  } // end for it
850 
851  // Maybe one that made it all the way to the end is the best?
852  for(std::list<Cand>::iterator live_it = live.begin(); live_it != live.end(); ++live_it){
853  if(live_it->NSegs() > best.NSegs()) best = *live_it;
854  }
855 
856  if(best.NSegs() == 0) ok = false;
857  // TODO: how can this even happen? But if it does we'll get stuck in an
858  // endless loop, so bail out of the rest of the event.
859  if(best.AllHits().empty()) ok = false;
860  return best;
861  }
set< int >::iterator it
geo::View_t GeoView() const
Definition: View.h:26
std::map< int, std::vector< Chunk > > ChunkMap
Definition: View.h:27
int NextPlane(unsigned int plane) const
Definition: View.cxx:111
ChunkMap MakeChunks() const
Definition: View.cxx:148
assert(nhit_max >=nhit_nbins)
std::vector< Chain > dt::DiscreteTracker::CheckAndFixSparseChain ( Chain chain) const
protected

Definition at line 864 of file DiscreteTracker_module.cc.

References dt::Chain::AllChunks(), dt::Chain::AllHits(), ana::assert(), productionTest::chain, make_goodruns_defs::chunks, dir, art::PtrVector< T >::empty(), febshutoff_auto::end, dt::Chain::ExtremalPlane(), fConsiderBadChans, dt::kDownstream, dt::kUpstream, dt::OtherDir(), file_size_ana::progress, dt::Chain::Truncate(), and registry_explorer::v.

Referenced by RecoDiscrete().

865  {
866  std::vector<Chain> remnants;
867 
868  const std::vector<Chunk> chunks = chain.AllChunks();
869  assert(!chunks.empty());
870 
871  std::set<int> planes;
872  for(unsigned int n = 0; n < chunks.size(); ++n)
873  planes.insert(chunks[n].Plane());
874 
875  bool progress;
876  do{
877  progress = false;
879  const Direction dir = OtherDir(end);
880 
881  int nTot = 0;
882  int nMiss = 0;
883 
884  int curPlane = chain.ExtremalPlane(end);
885 
886  if(chain.AllHits().empty()) break;
887  View v(chain.AllHits(), fConsiderBadChans);
888 
889  do{
890  ++nTot;
891 
892  const bool hitHere = planes.count(curPlane);
893  if(!hitHere) ++nMiss;
894 
895  // If we've accumulated more than 1/2 gaps, and have now come back
896  // into hits
897  if(hitHere && nMiss > 0){
898  // Want to allow HGH... and ban GH... and HGH
899  if(nTot == nMiss || (nMiss > 1 && nMiss*2 >= nTot-1)){
900  Chain remnant = chain;
901  remnant.Truncate(curPlane+((dir == kUpstream) ? +1 : -1), dir);
902  remnants.push_back(remnant);
903 
904  chain.Truncate(curPlane, end);
905  progress = true;
906  // Shouldn't look any further from this end, but try the other one
907  break;
908  }
909  }
910 
911  curPlane = v.AdjacentPlane(curPlane, dir);
912  } while(curPlane != -1);
913  } // end for end/dir
914  } while(progress);
915 
916  return remnants;
917  }
Definition: Chain.py:1
bool fConsiderBadChans
Passed to View.
bool progress
Insert implicit nodes #####.
std::void_t< T > n
chain
Check that an output directory exists.
Direction OtherDir(Direction d)
Definition: Types.h:7
Direction
Definition: Types.h:5
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 55 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumes(), T, and getGoodRuns4SAM::tag.

56  {
57  return collector_.consumes<T, BT>(tag);
58  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
ProductToken< T > consumes(InputTag const &)
double T
Definition: Xdiff_gwt.C:5
ConsumesCollector& art::ModuleBase::consumesCollector ( )
protectedinherited
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 69 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumesMany(), and T.

70  {
71  collector_.consumesMany<T, BT>();
72  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
double T
Definition: Xdiff_gwt.C:5
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 62 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumesView(), T, and getGoodRuns4SAM::tag.

63  {
64  return collector_.consumesView<T, BT>(tag);
65  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
ViewToken< Element > consumesView(InputTag const &)
double T
Definition: Xdiff_gwt.C:5
void art::detail::Producer::doBeginJob ( )
inherited
bool art::detail::Producer::doBeginRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited
void art::detail::Producer::doEndJob ( )
inherited
bool art::detail::Producer::doEndRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited
bool art::detail::Producer::doEvent ( EventPrincipal ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited
double dt::DiscreteTracker::DotScore ( const Cand c1,
const Cand c2,
bool  allowFlip = false 
) const
protected

Definition at line 359 of file DiscreteTracker_module.cc.

References ana::assert(), demo5::c1, demo5::c2, DotScoreHelper(), std::max(), dt::Cand::MaxGradient(), dt::Cand::MaybeFlip(), dt::Cand::MinGradient(), and dt::Cand::NSegs().

Referenced by SpliceIsSane(), and TryExtendCandOne().

361  {
362  // Can't even determine a direction
363  if(c1.NSegs() < 2 || c2.NSegs() < 2) return -1;
364 
365  double bestDot = -2;
366 
367  for(bool flip1: {false, true}){
368  if(flip1 && !allowFlip) continue;
369 
370  bool ok = true;
371  const Cand c1f = flip1 ? c1.MaybeFlip(ok) : c1;
372  if(!ok) continue;
373 
374  for(bool flip2: {false, true}){
375  if(flip2 && !allowFlip) continue;
376 
377  bool ok = true;
378  const Cand c2f = flip2 ? c2.MaybeFlip(ok) : c2;
379  if(!ok) continue;
380 
381  const double m1lo = c1f.MinGradient();
382  const double m1hi = c1f.MaxGradient();
383  const double m2lo = c2f.MinGradient();
384  const double m2hi = c2f.MaxGradient();
385 
386  assert(m1lo <= m1hi);
387  assert(m2lo <= m2hi);
388 
389  // Overlapping ideas of direction, best is perfect
390  if((m1hi > m2lo && m1lo < m2hi) ||
391  (m1lo < m2hi && m1hi > m2lo)) return 1;
392 
393  bestDot = std::max(bestDot, DotScoreHelper(m1hi, m2lo));
394  bestDot = std::max(bestDot, DotScoreHelper(m1lo, m2hi));
395  } // end for flip2
396  } // end for flip1
397 
398  return bestDot;
399  }
T max(const caf::Proxy< T > &a, T b)
double DotScoreHelper(double m1, double m2) const
c2
Definition: demo5.py:33
assert(nhit_max >=nhit_nbins)
c1
Definition: demo5.py:24
double dt::DiscreteTracker::DotScoreHelper ( double  m1,
double  m2 
) const
protected

Definition at line 353 of file DiscreteTracker_module.cc.

References std::sqrt().

Referenced by DotScore().

354  {
355  return fabs(1+m1*m2)/sqrt((1+m1*m1)*(1+m2*m2));
356  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
static constexpr Double_t m2
Definition: Munits.h:145
void dt::DiscreteTracker::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 195 of file DiscreteTracker_module.cc.

References fclose(), and fDebugFile.

196  {
198  }
FILE * fDebugFile
File that track details are written to.
fclose(fg1)
void dt::DiscreteTracker::ExtendChain ( const View v,
const dt::View::ChunkMap chunk_map,
Chain chain,
Direction  dir 
) const
protected

Definition at line 533 of file DiscreteTracker_module.cc.

References dt::Chain::Add(), fillBadChanDBTables::cur, dir, dt::Chain::ExtremalCand(), TryExtendChainOne(), and registry_explorer::v.

Referenced by RecoOneTrack().

536  {
537  Cand cur = chain.ExtremalCand(dir);
538 
539  while(TryExtendChainOne(cur, chain.ExtremalCand(dir), v, chunk_map, dir, false)){
540  chain.Add(cur, dir);
541  }
542  }
chain
Check that an output directory exists.
TDirectory * dir
Definition: macro.C:5
bool TryExtendChainOne(Cand &cur, const Cand &prev, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
rb::Vertex * dt::DiscreteTracker::FindVertex ( std::vector< rb::Track xTrks,
std::vector< rb::Track yTrks 
) const
protected

Toy vertex finder.

Uses the intersection of the two tracks in each view whose intersection is proportionally nearest their ends.

Definition at line 925 of file DiscreteTracker_module.cc.

References dt::CompareByLength(), geo::kX, geo::kY, geo::LineIntersection(), util::pythag(), and POTSpillRate::view.

Referenced by produce().

927  {
928  if(xTrks.empty() || yTrks.empty()) return 0;
929 
930  // Longest tracks first. Only used for fallback at the end
931  std::sort(xTrks.rbegin(), xTrks.rend(), CompareByLength);
932  std::sort(yTrks.rbegin(), yTrks.rend(), CompareByLength);
933 
934  double bestVx, bestVy, bestVz1, bestVz2;
935  // I don't believe there's a path that doesn't set these, but need this to
936  // shut the compiler up in optimized builds.
937  bestVx = bestVy = bestVz1 = bestVz2 = 0;
938 
939  for(geo::View_t view: {geo::kX, geo::kY}){
940  double bestScore = 1e10;
941  const std::vector<rb::Track>& trks = (view == geo::kX) ? xTrks : yTrks;
942  for(unsigned int idx1 = 0; idx1 < trks.size()-1; ++idx1){
943  if(trks[idx1].ExtentPlane() < 3) continue;
944  for(bool end1: {false, true}){
945  for(unsigned int idx2 = idx1+1; idx2 < trks.size(); ++idx2){
946  if(trks[idx2].ExtentPlane() < 3) continue;
947  for(bool end2: {false, true}){
948 
949  // Start and direction defining point of these two lines
950  const TVector3 s1 = end1 ? trks[idx1].Stop() : trks[idx1].Start();
951  const TVector3 s2 = end2 ? trks[idx2].Stop() : trks[idx2].Start();
952  const TVector3 d1 = trks[idx1].TrajectoryPoint(end1 ? trks[idx1].NTrajectoryPoints()-2 : 1);
953  const TVector3 d2 = trks[idx2].TrajectoryPoint(end2 ? trks[idx2].NTrajectoryPoints()-2 : 1);
954 
955  double vxy, vz;
956  const bool ok = geo::LineIntersection(s1[view], s1.Z(), d1[view], d1.Z(),
957  s2[view], s2.Z(), d2[view], d2.Z(),
958  vxy, vz);
959  if(!ok) continue;
960 
961  const double dist1 = util::pythag(s1[view]-vxy, s1.Z()-vz);
962  const double dist2 = util::pythag(s2[view]-vxy, s2.Z()-vz);
963 
964  // Favour both long tracks and hits close to the ends
965  const double score = dist1/trks[idx1].TotalLength()+dist2/trks[idx2].TotalLength();
966  if(score < bestScore){
967  bestScore = score;
968  if(view == geo::kX){
969  bestVx = vxy;
970  bestVz1 = vz;
971  }
972  else{
973  bestVy = vxy;
974  bestVz2 = vz;
975  }
976  } // end if better
977  } // end for end2
978  } // end for idx2
979  } // end for end1
980  } // end for idx1
981  // Too far from the ends, fallback
982  if(bestScore > 6){
983  if(view == geo::kX){
984  bestVx = trks[0].Start().X();
985  bestVz1 = trks[0].Start().Z();
986  }
987  else{
988  bestVy = trks[0].Start().Y();
989  bestVz2 = trks[0].Start().Z();
990  }
991  }
992  } // end for view
993 
994  // Time is an afterthought
995  return new rb::Vertex(bestVx, bestVy, (bestVz1+bestVz2)/2, xTrks[0].MinTNS());
996  }
bool LineIntersection(double x0, double y0, double x1, double y1, double X0, double Y0, double X1, double Y1, double &x, double &y)
Find the intersection between two line segments.
Definition: Geo.cxx:184
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
bool CompareByLength(const rb::Track &a, const rb::Track &b)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
std::array<std::vector<ProductInfo>, NumBranchTypes> const& art::ModuleBase::getConsumables ( ) const
inherited
bool dt::DiscreteTracker::HighActivity ( const art::Ptr< rb::Cluster > &  clust) const
protected

Cut to prevent choking on large air shower events.

Adapted from the version in KalmanTrack

Definition at line 267 of file DiscreteTracker_module.cc.

References geom(), chaninfo::BadChanList::IsBad(), rb::Cluster::NCell(), geo::PlaneGeo::Ncells(), geo::GeometryBase::NPlanes(), and geo::GeometryBase::Plane().

Referenced by produce().

268  {
271 
272  // get the total number of non masked off channels by looping over all
273  // planes/cells in the detector
274  int numGoodChannels = 0;
275 
276  for(unsigned int iPlane = 0; iPlane < geom->NPlanes(); ++iPlane){
277  const geo::PlaneGeo* p = geom->Plane(iPlane);
278  for(unsigned int iCell = 0; iCell < p->Ncells(); ++iCell){
279  if(!bcl->IsBad(iPlane, iCell)) ++numGoodChannels;
280  }
281  }
282 
283  return clust->NCell() > 0.15*numGoodChannels;
284  }
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const char * p
Definition: xmltok.h:285
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const PlaneGeo * Plane(unsigned int i) const
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
void geom(int which=0)
Definition: geom.C:163
unsigned int NPlanes() const
bool IsBad(int plane, int cell)
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 76 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsume(), T, and getGoodRuns4SAM::tag.

77  {
78  return collector_.mayConsume<T, BT>(tag);
79  }
ProductToken< T > mayConsume(InputTag const &)
ConsumesCollector collector_
Definition: ModuleBase.h:50
double T
Definition: Xdiff_gwt.C:5
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 90 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsumeMany(), and T.

91  {
93  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
double T
Definition: Xdiff_gwt.C:5
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 83 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsumeView(), T, and getGoodRuns4SAM::tag.

84  {
85  return collector_.mayConsumeView<T, BT>(tag);
86  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
ViewToken< Element > mayConsumeView(InputTag const &)
double T
Definition: Xdiff_gwt.C:5
void dt::DiscreteTracker::MergeParasiteTracks ( std::vector< rb::Track > &  tracks) const
protected

Definition at line 999 of file DiscreteTracker_module.cc.

References a, b, rb::CellHit::Cell(), rb::Cluster::Cell(), rb::HitMap::CellExists(), DEFINE_ART_MODULE(), rb::Cluster::ExtentPlane(), MECModelEnuComparisons::i, cafExposure::match, rb::Cluster::NCell(), and rb::CellHit::Plane().

Referenced by RecoDiscrete().

1000  {
1001  // Worried this could be too aggressive around the vertex
1002 
1003  const unsigned int N = tracks.size();
1004  if(N < 2) return;
1005 
1006  // Longest first
1007  std::sort(tracks.begin(), tracks.end(),
1008  [](const rb::Track& a, const rb::Track& b){
1009  return a.ExtentPlane() > b.ExtentPlane();
1010  });
1011 
1012  for(unsigned int hostIdx = 0; hostIdx < N-1; ++hostIdx){
1013  const rb::HitMap hmap(&tracks[hostIdx]);
1014  for(unsigned int paraIdx = hostIdx+1; paraIdx < N; ++paraIdx){
1015  const rb::Track& para = tracks[paraIdx];
1016  bool match = true;
1017  for(unsigned int i = 0; i < para.NCell(); ++i){
1018  const art::Ptr<rb::CellHit> chit = para.Cell(i);
1019  if(!hmap.CellExists(chit->Plane(), chit->Cell()-1) &&
1020  !hmap.CellExists(chit->Plane(), chit->Cell()+1)){
1021  match = false;
1022  break;
1023  }
1024  } // end for i
1025  if(match){
1026  // Make sure there are no duplicates in the final host track
1027  rb::HitMap hmap(tracks[hostIdx].AllCells());
1028  for(unsigned int n = 0; n < para.NCell(); ++n){
1029  const art::Ptr<rb::CellHit> chit = para.Cell(n);
1030  if(!hmap.CellExists(chit->Plane(), chit->Cell()))
1031  tracks[hostIdx].Add(chit);
1032  }
1033  tracks.erase(tracks.begin()+paraIdx);
1034  MergeParasiteTracks(tracks);
1035  return;
1036  }
1037  } // end for paraIdx
1038  } // end for hostIdx
1039  }
void MergeParasiteTracks(std::vector< rb::Track > &tracks) const
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned short Plane() const
Definition: CellHit.h:39
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
unsigned short Cell() const
Definition: CellHit.h:40
std::void_t< T > n
const double a
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
const hit & b
Definition: hits.cxx:21
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
ModuleDescription const& art::ModuleBase::moduleDescription ( ) const
inherited
std::vector< Segment > dt::DiscreteTracker::PossibleNextSegs ( const View v,
const dt::View::ChunkMap chunk_map,
const Cand cand,
Direction  dir,
std::vector< Chunk > &  chunks_out,
bool  shw,
bool  exactHits 
) const
protected

Definition at line 288 of file DiscreteTracker_module.cc.

References dt::View::AdjacentPlane(), make_associated_cosmic_defs::chunk, make_goodruns_defs::chunks, dir, dt::Cand::ExtremalPlane(), dt::Cand::FirstChunk(), dt::Cand::FirstPlane(), dt::Cand::FirstSeg(), dt::Chunk::GetSegs(), dt::Chunk::IsWorthShowering(), it, dt::kDownstream, dt::kUpstream, dt::Cand::LastChunk(), dt::Cand::LastPlane(), dt::Cand::LastSeg(), dt::Segment::left, runNovaSAM::ret, dt::Chunk::SetShowerChunk(), dt::Chunk::SetUp(), dt::Chunk::Up(), and dt::Cand::Up().

Referenced by TryExtendCandOne().

293  {
294  std::vector<Segment> ret;
295 
296  // If the current chunk is only partially taken then the next thing to do
297  // must be to try and complete it.
298  if(dir == kUpstream && !cand.FirstSeg().left){
299  ret.push_back(cand.FirstChunk().GetSegs().first);
300  return ret;
301  }
302 
303  if(dir == kDownstream && cand.LastSeg().left){
304  ret.push_back(cand.LastChunk().GetSegs().second);
305  return ret;
306  }
307 
308  int nextPlane = v.AdjacentPlane(cand.ExtremalPlane(dir), dir);
309  if(nextPlane == -1){
310  if(!exactHits) return ret;
311  if(dir == kDownstream){
312  nextPlane = chunk_map.begin()->first;
313  if(nextPlane <= int(cand.LastPlane())) return ret;
314  }
315  else{
316  nextPlane = chunk_map.rbegin()->first;
317  if(nextPlane >= int(cand.FirstPlane())) return ret;
318  }
319  }
320 
321  dt::View::ChunkMap::const_iterator it = chunk_map.find(nextPlane);
322  if(exactHits){
323  while(it == chunk_map.end()){
324  nextPlane = v.AdjacentPlane(nextPlane, dir);
325  if(nextPlane == -1) return ret;
326  it = chunk_map.find(nextPlane);
327  }
328  }
329  if(it == chunk_map.end()) return ret;
330 
331  const std::vector<Chunk>& chunks = it->second;
332  const unsigned int chunkMax = chunks.size();
333 
334  for(unsigned int chunkIdx = 0; chunkIdx < chunkMax; ++chunkIdx){
335  Chunk chunk = chunks[chunkIdx];
336  if(exactHits) chunk.SetUp(cand.Up());
337 
338  if(chunk.Up() != cand.Up()) continue;
339  if(shw && !chunk.IsWorthShowering()) continue;
340  if(shw) chunk.SetShowerChunk();
341 
342  chunks_out.push_back(chunk);
343  if(dir == kUpstream)
344  ret.push_back(chunk.GetSegs().second);
345  else
346  ret.push_back(chunk.GetSegs().first);
347  }
348 
349  return ret;
350  }
set< int >::iterator it
TDirectory * dir
Definition: macro.C:5
void dt::DiscreteTracker::produce ( art::Event evt)
virtual

Splits hits according to slice and view for feeding to RecoDiscrete. Writes the returned tracks into the event.

Implements art::EDProducer.

Definition at line 201 of file DiscreteTracker_module.cc.

References rb::Cluster::Cell(), util::CreateAssn(), art::DataViewImpl::event(), fClusterInput, fDebugFile, FindVertex(), fObeyPreselection, art::DataViewImpl::getByLabel(), HighActivity(), rb::IsFiltered(), rb::Cluster::IsNoise(), it, geo::kX, geo::kY, rb::Cluster::NCell(), art::PtrVector< T >::push_back(), art::DataViewImpl::put(), RecoDiscrete(), and POTSpillRate::view.

202  {
203  std::unique_ptr<std::vector<rb::Track> > debug_seeds_col(new std::vector<rb::Track>);
204  std::unique_ptr<std::vector<rb::Track> > debug_unspliced_col(new std::vector<rb::Track>);
205  std::unique_ptr<std::vector<rb::Track> > debug_chains_col(new std::vector<rb::Track>);
206  std::unique_ptr<std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
207  std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
208  // Associate tracks with the slice they came from
209  std::unique_ptr<art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
210 
212  evt.getByLabel(fClusterInput, clusts);
213 
214  const int clustMax = clusts->size();
215  for(int clustIdx = 0; clustIdx < clustMax; ++clustIdx){
216  art::Ptr<rb::Cluster> clust(clusts, clustIdx);
217  if(fObeyPreselection && rb::IsFiltered(evt, clust)) continue;
218 
219  if(clust->IsNoise()) continue;
220 
221  if(HighActivity(clust)){
222  mf::LogWarning("DiscreteTracker") << "Skipping event due to very high activity";
223  continue;
224  }
225 
226  std::vector<rb::Track> tracks[2];
227  for(geo::View_t view: {geo::kX, geo::kY}){
228  const char viewC = (view == geo::kX) ? 'x' : 'y';
229 
231  const int Nx = clust->NCell(view);
232  for(int n = 0; n < Nx; ++n) cells.push_back(clust->Cell(view, n));
233 
234  fprintf(fDebugFile, "%d: %c: ", evt.event(), viewC);
235 
236  tracks[view] = RecoDiscrete(cells,
237  debug_seeds_col.get(),
238  debug_unspliced_col.get(),
239  debug_chains_col.get());
240 
241  for(unsigned int n = 0; n < tracks[view].size(); ++n)
242  fprintf(fDebugFile, "%u ", tracks[view][n].NCell());
243  fprintf(fDebugFile, "\n");
244 
245  for(std::vector<rb::Track>::const_iterator it = tracks[view].begin();
246  it != tracks[view].end(); ++it){
247  trackcol->push_back(*it);
248  // Associate most recently added track with the slice
249  util::CreateAssn(evt, *trackcol, clust, *assns);
250  } // end for it
251  } // end for view
252 
253  const rb::Vertex* vtx = FindVertex(tracks[geo::kX], tracks[geo::kY]);
254  if(vtx) vtxcol->push_back(*vtx);
255  delete vtx;
256  } // end for clustIdx
257 
258  evt.put(std::move(debug_seeds_col), "seeds");
259  evt.put(std::move(debug_unspliced_col), "unspliced");
260  evt.put(std::move(debug_chains_col), "chains");
261  evt.put(std::move(trackcol));
262  evt.put(std::move(vtxcol));
263  evt.put(std::move(assns));
264  }
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
EventNumber_t event() const
set< int >::iterator it
FILE * fDebugFile
File that track details are written to.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
bool HighActivity(const art::Ptr< rb::Cluster > &clust) const
Cut to prevent choking on large air shower events.
std::void_t< T > n
rb::Vertex * FindVertex(std::vector< rb::Track > xTrks, std::vector< rb::Track > yTrks) const
Toy vertex finder.
bool fObeyPreselection
Check rb::IsFiltered?
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
std::vector< rb::Track > RecoDiscrete(const art::PtrVector< rb::CellHit > &cells, std::vector< rb::Track > *seeds_out, std::vector< rb::Track > *unspliced_out, std::vector< rb::Track > *chains_out) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
std::string fClusterInput
Module name to look for Clusters under.
std::vector< rb::Track > dt::DiscreteTracker::RecoDiscrete ( const art::PtrVector< rb::CellHit > &  cells,
std::vector< rb::Track > *  seeds_out,
std::vector< rb::Track > *  unspliced_out,
std::vector< rb::Track > *  chains_out 
) const
protected

Definition at line 563 of file DiscreteTracker_module.cc.

References dt::Chain::AllHits(), productionTest::chain, makeDatasetsPage::chains, dt::View::Channel(), CheckAndFixSparseChain(), gen_hdf5record::debug, art::PtrVector< T >::empty(), fConsiderBadChans, dt::Chain::FirstCand(), fOnlyOneTrack, hits(), dt::Chain::LastCand(), dt::View::MakeChunks(), dt::Channel::MarkHitFound(), MergeParasiteTracks(), dt::Cand::NSegs(), RecoOneTrack(), runNovaSAM::ret, art::PtrVector< T >::size(), SpliceChains(), dt::Chain::ToTrack(), dt::Chain::TrimEnds(), and registry_explorer::v.

Referenced by produce().

567  {
568  std::vector<rb::Track> ret;
569  if(cells.empty()) return ret;
570 
571  std::vector<Chain> chains;
572 
573  dt::View v(cells, fConsiderBadChans);
574 
575  while(true){
576  // TODO, this should be shared
577  const dt::View::ChunkMap chunk_map = v.MakeChunks();
578 
579  bool ok;
580  Chain chain = RecoOneTrack(v, chunk_map, seeds_out, ok);
581 
582  if(!ok) break;
583 
584  chain.TrimEnds();
585  // TODO: how can this happen?
586  if(chain.FirstCand().NSegs() < 2) break;
587  if(chain.LastCand().NSegs() < 2) break;
588 
589  chains.push_back(chain);
590  unspliced_out->push_back(chain.ToTrack());
591 
592  // We're going to try and find another track, but first remove all the
593  // hits we've already explained.
594  art::PtrVector<rb::CellHit> hits = chain.AllHits();
595  for(unsigned int n = 0; n < hits.size(); ++n){
596  v.Channel(hits[n]->Plane(), hits[n]->Cell()).MarkHitFound();
597  }
598 
599  if(fOnlyOneTrack) break;
600  } // end while
601 
602  SpliceChains(chains);
603 
604  unsigned int N = chains.size();
605  for(unsigned int n = 0; n < N; ++n){
606  const std::vector<Chain> remnants = CheckAndFixSparseChain(chains[n]);
607  chains.insert(chains.end(), remnants.begin(), remnants.end());
608  }
609 
610  N = chains.size();
611  for(unsigned int n = 0; n < N; ++n){
612  // Comes out of Truncate() at least sometimes
613  if(chains[n].AllHits().empty()){
614  mf::LogWarning("DiscreteTracker") << "Track with no hits. Skipping.";
615  continue;
616  }
617 
618  if(chains[n].FirstSeg().zL > chains[n].LastSeg().zL){
619  mf::LogWarning("DiscreteTracker") << "Reversed track. Skipping.";
620  continue;
621  }
622 
623  ret.push_back(chains[n].ToTrack());
624  const std::vector<rb::Track> debug = chains[n].ToDebugTracks();
625  chains_out->insert(chains_out->end(), debug.begin(), debug.end());
626  }
627 
628  MergeParasiteTracks(ret);
629 
630  return ret;
631  }
void MergeParasiteTracks(std::vector< rb::Track > &tracks) const
Representation of the channels in one view.
Definition: View.h:22
Definition: Chain.py:1
bool fConsiderBadChans
Passed to View.
Chain RecoOneTrack(const View &v, const dt::View::ChunkMap &chunk_map, std::vector< rb::Track > *seeds_out, bool &ok) const
std::map< int, std::vector< Chunk > > ChunkMap
Definition: View.h:27
std::vector< Chain > CheckAndFixSparseChain(Chain &chain) const
void hits()
Definition: readHits.C:15
std::void_t< T > n
void SpliceChains(std::vector< Chain > &chains) const
chain
Check that an output directory exists.
bool empty() const
Definition: PtrVector.h:330
size_type size() const
Definition: PtrVector.h:302
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool fOnlyOneTrack
Stop after the first track in each view.
void dt::DiscreteTracker::reconfigure ( const fhicl::ParameterSet pset)

Definition at line 172 of file DiscreteTracker_module.cc.

References fclose(), fClusterInput, fConsiderBadChans, fDebugFile, fMinDotProduct, fObeyPreselection, fOnlyOneTrack, fVerticalTrackMinCells, fhicl::ParameterSet::get(), and string.

Referenced by DiscreteTracker().

173  {
174  fClusterInput = pset.get<std::string>("ClusterInput");
175  fConsiderBadChans = pset.get<bool>("ConsiderBadChans");
176  fMinDotProduct = pset.get<double>("MinDotProduct");
177  fVerticalTrackMinCells = pset.get<int>("VerticalTrackMinCells");
178  fOnlyOneTrack = pset.get<bool>("OnlyOneTrack");
179  fObeyPreselection = pset.get<bool>("ObeyPreselection");
180 
182  fDebugFile = fopen(pset.get<std::string>("DebugFile").c_str(), "w");
183  }
FILE * fDebugFile
File that track details are written to.
bool fConsiderBadChans
Passed to View.
double fMinDotProduct
Joined Cands that only have one chunk in common must have a dot-product better than this...
fclose(fg1)
T get(std::string const &key) const
Definition: ParameterSet.h:231
bool fObeyPreselection
Check rb::IsFiltered?
int fVerticalTrackMinCells
The number of adjacent cells in a single plane, otherwise not fit into any track, required to declare...
bool fOnlyOneTrack
Stop after the first track in each view.
std::string fClusterInput
Module name to look for Clusters under.
enum BeamMode string
Chain dt::DiscreteTracker::RecoOneTrack ( const View v,
const dt::View::ChunkMap chunk_map,
std::vector< rb::Track > *  seeds_out,
bool &  ok 
) const
protected

Definition at line 545 of file DiscreteTracker_module.cc.

References BestSeed(), productionTest::chain, ExtendChain(), dt::kDownstream, dt::kUpstream, seed, and dt::Cand::ToTrack().

Referenced by RecoDiscrete().

549  {
550  const Cand& seed = BestSeed(v, ok);
551  if(!ok) return Chain();
552  seeds_out->push_back(seed.ToTrack());
553 
554  Chain chain(seed);
555  ExtendChain(v, chunk_map, chain, kUpstream);
556  ExtendChain(v, chunk_map, chain, kDownstream);
557  ok = true;
558  return chain;
559  }
Definition: Chain.py:1
void ExtendChain(const View &v, const dt::View::ChunkMap &chunk_map, Chain &chain, Direction dir) const
unsigned int seed
Definition: runWimpSim.h:102
chain
Check that an output directory exists.
Cand BestSeed(const dt::View &v, bool &ok) const
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  )
inherited
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited
void dt::DiscreteTracker::SpliceChains ( std::vector< Chain > &  chains) const
protected

Definition at line 672 of file DiscreteTracker_module.cc.

References dt::Chain::Add(), dt::Chain::AllHits(), dt::Chain::AsChunkMap(), dt::SpliceDesc::childIdx, dir, dt::Chain::ExtremalCand(), dt::Chain::ExtremalPlane(), dt::Cand::ExtremalPlane(), fConsiderBadChans, dt::Chain::FirstPlane(), MECModelEnuComparisons::i, calib::j, dt::kDownstream, dt::kUpstream, dt::Chain::LastPlane(), dt::Chain::NHitPlanes(), dt::OtherDir(), file_size_ana::parent, dt::SpliceDesc::parentIdx, make_associated_cosmic_defs::sd, dt::SpliceDesc::spliced, SpliceIsSane(), TryExtendChainOne(), and registry_explorer::v.

Referenced by RecoDiscrete().

673  {
674  while(true){
675  std::vector<int> parentIdxs(chains.size());
676  for(unsigned int n = 0; n < chains.size(); ++n) parentIdxs[n] = n;
677 
678  std::vector<int> childIdxs(chains.size());
679  for(unsigned int n = 0; n < chains.size(); ++n) childIdxs[n] = n;
680 
681  // Keep the best splice first
682  std::multiset<SpliceDesc> splices;
683 
684  for(unsigned int i = 0; i < chains.size(); ++i){
685  const unsigned int parentIdx = parentIdxs[i];
686  const Chain& parent = chains[parentIdx];
687 
689  for(unsigned int j = 0; j < chains.size(); ++j){
690  const unsigned int childIdx = childIdxs[j];
691 
692  if(childIdx == parentIdx) continue;
693 
694  const Chain& child = chains[childIdx];
695 
696  int gap;
697  if(dir == kDownstream)
698  gap = child.FirstPlane() - parent.LastPlane();
699  else
700  gap = parent.FirstPlane() - child.LastPlane();
701 
702  // Have to get to correct side
703  if(gap < 0) continue;
704 
705  // Bail out early if we know we can't beat the best splice
706  if(!splices.empty() && gap > splices.begin()->gap) continue;
707 
708  if(!SpliceIsSane(parent, child)) continue;
709 
710  View v(child.AllHits(), fConsiderBadChans);
711  // Doing it this way preserves the "shower" flag on chunks
712  const dt::View::ChunkMap chunk_map = child.AsChunkMap();
713 
714  Cand join = parent.ExtremalCand(dir);
715  if(!TryExtendChainOne(join, parent.ExtremalCand(dir),
716  v, chunk_map,
717  dir, true)) continue;
718 
719  // Also need to check that the join is willing to be made in the
720  // other direction. Otherwise you can get a one-hit cand on one
721  // side, that's willing to match almost anything, in crazy ways.
722  View pview(parent.AllHits(), fConsiderBadChans);
723  const dt::View::ChunkMap pmap = parent.AsChunkMap();
724  Cand join2 = child.ExtremalCand(OtherDir(dir));
725  if(!TryExtendChainOne(join2, child.ExtremalCand(OtherDir(dir)),
726  pview, pmap, OtherDir(dir), true)) continue;
727 
728  Chain spliced = parent;
729  spliced.Add(join, dir);
730  // The splice can often subsume the whole child chain
731  if(child.ExtremalPlane(dir) != join.ExtremalPlane(dir)){
732  spliced.Add(child, dir);
733  }
734 
735  splices.insert(SpliceDesc(childIdx, parentIdx, spliced, gap,
736  child.NHitPlanes()+parent.NHitPlanes()));
737  } // end for j (childIdx)
738  } // end for dir
739  } // end for i (parentIdx)
740 
741  if(splices.empty()) return;
742 
743  const SpliceDesc sd = *splices.begin();
744  chains[sd.parentIdx] = sd.spliced;
745  chains.erase(chains.begin()+sd.childIdx);
746  }
747  }
Definition: Chain.py:1
bool fConsiderBadChans
Passed to View.
std::map< int, std::vector< Chunk > > ChunkMap
Definition: View.h:27
bool SpliceIsSane(const Chain &a, const Chain &b) const
std::void_t< T > n
const double j
Definition: BetheBloch.cxx:29
Direction OtherDir(Direction d)
Definition: Types.h:7
Direction
Definition: Types.h:5
TDirectory * dir
Definition: macro.C:5
bool TryExtendChainOne(Cand &cur, const Cand &prev, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
bool dt::DiscreteTracker::SpliceIsSane ( const Chain a,
const Chain b 
) const
protected

Definition at line 634 of file DiscreteTracker_module.cc.

References DotScore(), dt::Chain::FirstCand(), dt::Chain::FirstPlane(), fMinDotProduct, dt::Chain::LastCand(), and dt::Chain::LastPlane().

Referenced by SpliceChains().

635  {
636  // Must not overlap. Can't jump more than half the whole total length
637 
638  if(a.LastPlane() < b.FirstPlane()){
639  if(DotScore(a.LastCand(), b.FirstCand(), true) < fMinDotProduct) return false;
640 
641  return (2*(b.FirstPlane()-a.LastPlane()) < b.LastPlane()-a.FirstPlane());
642  }
643  if(b.LastPlane() < a.FirstPlane()){
644  if(DotScore(b.LastCand(), a.FirstCand(), true) < fMinDotProduct) return false;
645 
646  return (2*(a.FirstPlane()-b.LastPlane()) < a.LastPlane()-b.FirstPlane());
647  }
648 
649  return false;
650  }
double DotScore(const Cand &c1, const Cand &c2, bool allowFlip=false) const
double fMinDotProduct
Joined Cands that only have one chunk in common must have a dot-product better than this...
const double a
const hit & b
Definition: hits.cxx:21
bool dt::DiscreteTracker::TryExtendCand ( Cand cand,
const Cand old_cand,
const View v,
const dt::View::ChunkMap chunk_map,
Direction  dir,
bool  exactHits 
) const
protected

Definition at line 463 of file DiscreteTracker_module.cc.

References dt::Cand::ExtremalChunk(), dt::Cand::IsExtremalChunkComplete(), dt::Chunk::IsShowerChunk(), dt::Chunk::IsWorthShowering(), dt::Cand::MarkExtremalChunkShower(), dt::Cand::MaybeFlip(), dt::Cand::NSegs(), and TryExtendCandOne().

Referenced by TryExtendChainOne().

469  {
470  const Cand bak = cand;
471 
472  for(bool shwNext: {false, true}){
473  if(shwNext && !cand.IsExtremalChunkComplete(dir)) continue;
474  for(bool shwLast: {false, true}){
475  if(shwLast && cand.ExtremalChunk(dir).IsShowerChunk()) continue;
476  if(shwLast && !cand.ExtremalChunk(dir).IsWorthShowering()) continue;
477 
478  if(shwLast && cand.NSegs() <= 2) continue;
479 
480  if(shwLast) cand.MarkExtremalChunkShower(dir);
481  // Reset back by hand at end of loop. Careful about early exits
482 
483  bool any = false;
484  while(TryExtendCandOne(cand, v, chunk_map, dir, old_cand, shwNext, exactHits)) any = true;
485 
486  if(any) return true;
487 
488  // TODO: can probably hoist this higher with care
489  bool ok;
490  Cand flipped = cand.MaybeFlip(ok);
491  if(!ok){
492  if(shwLast) cand = bak;
493  continue;
494  }
495 
496  while(TryExtendCandOne(flipped, v, chunk_map, dir, old_cand, shwNext, exactHits)) any = true;
497 
498  if(any){
499  cand = flipped;
500  return true;
501  }
502 
503  if(shwLast) cand = bak;
504  } // end for shwLast
505  } // end for shwNext
506 
507  cand = bak;
508  return false;
509  }
bool TryExtendCandOne(Cand &cand, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, const Cand &prev_dir, bool shw, bool exactHits) const
TDirectory * dir
Definition: macro.C:5
bool dt::DiscreteTracker::TryExtendCandOne ( Cand cand,
const View v,
const dt::View::ChunkMap chunk_map,
Direction  dir,
const Cand prev_dir,
bool  shw,
bool  exactHits 
) const
protected

Definition at line 402 of file DiscreteTracker_module.cc.

References std::acos(), dt::Cand::AddChunk(), make_goodruns_defs::chunks, DotScore(), fAngle, fMinDotProduct, fSinglesAngle, dt::kDownstream, dt::kUpstream, dt::Segment::left, dt::Cand::NSegs(), PossibleNextSegs(), and dt::Cand::TryAddSeg().

Referenced by TryExtendCand().

408  {
409  // TODO: hoist this out of while loop below
410  std::vector<Chunk> chunks;
411  const std::vector<Segment> segs = PossibleNextSegs(v, chunk_map,
412  cand, dir, chunks, shw,
413  exactHits);
414  if(segs.empty()) return false;
415 
416  // Backup of cand, so we can restore it in case TryAddSeg succeeds
417  Cand bak = cand;
418 
419  std::vector<Cand> extended;
420 
421  const unsigned int segMax = segs.size();
422  for(unsigned int segIdx = 0; segIdx < segMax; ++segIdx){
423  const Segment& seg = segs[segIdx];
424 
425  if(cand.TryAddSeg(seg, dir)){
426 
427  // May need to add the chunk to match
428  if((dir == kUpstream && !seg.left) ||
429  (dir == kDownstream && seg.left)){
430  cand.AddChunk(chunks[segIdx], dir);
431  }
432 
433  const bool smallOverlap = (bak.NSegs() <= 2);
434  if(true){//smallOverlap){
435  // If there's only a small overlap, ensure dot product is sensible
436  const double dotscore = DotScore(prev_dir, cand);
437  fAngle->Fill(acos(dotscore));
438  if(smallOverlap) fSinglesAngle->Fill(acos(dotscore));
439  if(smallOverlap && dotscore < fMinDotProduct){
440  cand = bak;
441  continue;
442  }
443  }
444 
445  // Otherwise it's good
446  extended.push_back(cand);
447  cand = bak;
448  } // end if added seg
449  } // end for segIdx
450 
451  if(extended.empty()) return false;
452 
453  // TODO: Arbitrary for now
454  cand = extended.front();
455  if(extended.size() > 1){
456  // std::cout << "***WARNING: MORE THAN ONE EXTENSION CANDIDATE (" << extended.size() << "). SHOULD HAVE TRIED TO MAKE A SMARTER DECISION***" << std::endl;
457  }
458 
459  return true;
460  }
double DotScore(const Cand &c1, const Cand &c2, bool allowFlip=false) const
TH1 * fAngle
Opening angle of all cands considered for joining.
T acos(T number)
Definition: d0nt_math.hpp:54
double fMinDotProduct
Joined Cands that only have one chunk in common must have a dot-product better than this...
std::vector< Segment > PossibleNextSegs(const View &v, const dt::View::ChunkMap &chunk_map, const Cand &cand, Direction dir, std::vector< Chunk > &chunks_out, bool shw, bool exactHits) const
TH1 * fSinglesAngle
For all cands with only one overlapping chunk.
TDirectory * dir
Definition: macro.C:5
bool dt::DiscreteTracker::TryExtendChainOne ( Cand cur,
const Cand prev,
const View v,
const dt::View::ChunkMap chunk_map,
Direction  dir,
bool  exactHits 
) const
protected

Definition at line 512 of file DiscreteTracker_module.cc.

References dt::Cand::AllHits(), art::PtrVector< T >::empty(), dt::Cand::IsAllShower(), dt::Cand::NChunks(), dt::Cand::NSegs(), dt::Cand::ShortenOne(), and TryExtendCand().

Referenced by ExtendChain(), and SpliceChains().

518  {
519  while(true){
520  if(cur.NSegs() < 2) return false;
521  if(TryExtendCand(cur, prev, v, chunk_map, dir, exactHits)) return true;
522 
523  cur.ShortenOne(dir);
524 
525  // Too much freedom
526  if(cur.IsAllShower()) return false;
527  if(cur.NChunks() == 0) return false; // all dead
528  if(cur.AllHits().empty()) return false;
529  }
530  }
void prev()
Definition: show_event.C:91
TDirectory * dir
Definition: macro.C:5
bool TryExtendCand(Cand &cand, const Cand &old_cand, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
std::string art::EDProducer::workerType ( ) const
inherited

Member Data Documentation

TH1* dt::DiscreteTracker::fAngle
protected

Opening angle of all cands considered for joining.

Definition at line 143 of file DiscreteTracker_module.cc.

Referenced by beginJob(), and TryExtendCandOne().

std::string dt::DiscreteTracker::fClusterInput
protected

Module name to look for Clusters under.

Definition at line 130 of file DiscreteTracker_module.cc.

Referenced by produce(), and reconfigure().

bool dt::DiscreteTracker::fConsiderBadChans
protected

Passed to View.

Definition at line 131 of file DiscreteTracker_module.cc.

Referenced by CheckAndFixSparseChain(), RecoDiscrete(), reconfigure(), and SpliceChains().

FILE* dt::DiscreteTracker::fDebugFile
protected

File that track details are written to.

Definition at line 146 of file DiscreteTracker_module.cc.

Referenced by endJob(), produce(), and reconfigure().

double dt::DiscreteTracker::fMinDotProduct
protected

Joined Cands that only have one chunk in common must have a dot-product better than this.

Definition at line 136 of file DiscreteTracker_module.cc.

Referenced by reconfigure(), SpliceIsSane(), and TryExtendCandOne().

bool dt::DiscreteTracker::fObeyPreselection
protected

Check rb::IsFiltered?

Definition at line 132 of file DiscreteTracker_module.cc.

Referenced by produce(), and reconfigure().

bool dt::DiscreteTracker::fOnlyOneTrack
protected

Stop after the first track in each view.

Definition at line 141 of file DiscreteTracker_module.cc.

Referenced by RecoDiscrete(), and reconfigure().

TH1* dt::DiscreteTracker::fSinglesAngle
protected

For all cands with only one overlapping chunk.

Definition at line 144 of file DiscreteTracker_module.cc.

Referenced by beginJob(), and TryExtendCandOne().

int dt::DiscreteTracker::fVerticalTrackMinCells
protected

The number of adjacent cells in a single plane, otherwise not fit into any track, required to declare them a vertical track.

Definition at line 139 of file DiscreteTracker_module.cc.

Referenced by reconfigure().


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