Public Types | Public Member Functions | Static 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::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::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 ()
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

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
 
CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

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

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

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

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

Definition at line 166 of file DiscreteTracker_module.cc.

167  {
168  }

Member Function Documentation

void dt::DiscreteTracker::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 185 of file DiscreteTracker_module.cc.

References fAngle, fSinglesAngle, and art::TFileDirectory::make().

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

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

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

Referenced by RecoDiscrete().

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

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
double dt::DiscreteTracker::DotScore ( const Cand c1,
const Cand c2,
bool  allowFlip = false 
) const
protected

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

360  {
361  // Can't even determine a direction
362  if(c1.NSegs() < 2 || c2.NSegs() < 2) return -1;
363 
364  double bestDot = -2;
365 
366  for(bool flip1: {false, true}){
367  if(flip1 && !allowFlip) continue;
368 
369  bool ok = true;
370  const Cand c1f = flip1 ? c1.MaybeFlip(ok) : c1;
371  if(!ok) continue;
372 
373  for(bool flip2: {false, true}){
374  if(flip2 && !allowFlip) continue;
375 
376  bool ok = true;
377  const Cand c2f = flip2 ? c2.MaybeFlip(ok) : c2;
378  if(!ok) continue;
379 
380  const double m1lo = c1f.MinGradient();
381  const double m1hi = c1f.MaxGradient();
382  const double m2lo = c2f.MinGradient();
383  const double m2hi = c2f.MaxGradient();
384 
385  assert(m1lo <= m1hi);
386  assert(m2lo <= m2hi);
387 
388  // Overlapping ideas of direction, best is perfect
389  if((m1hi > m2lo && m1lo < m2hi) ||
390  (m1lo < m2hi && m1hi > m2lo)) return 1;
391 
392  bestDot = std::max(bestDot, DotScoreHelper(m1hi, m2lo));
393  bestDot = std::max(bestDot, DotScoreHelper(m1lo, m2hi));
394  } // end for flip2
395  } // end for flip1
396 
397  return bestDot;
398  }
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 352 of file DiscreteTracker_module.cc.

References stan::math::fabs(), and std::sqrt().

Referenced by DotScore().

353  {
354  return fabs(1+m1*m2)/sqrt((1+m1*m1)*(1+m2*m2));
355  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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 194 of file DiscreteTracker_module.cc.

References fclose(), and fDebugFile.

195  {
197  }
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 532 of file DiscreteTracker_module.cc.

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

Referenced by RecoOneTrack().

535  {
536  Cand cur = chain.ExtremalCand(dir);
537 
538  while(TryExtendChainOne(cur, chain.ExtremalCand(dir), v, chunk_map, dir, false)){
539  chain.Add(cur, dir);
540  }
541  }
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 924 of file DiscreteTracker_module.cc.

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

Referenced by produce().

926  {
927  if(xTrks.empty() || yTrks.empty()) return 0;
928 
929  // Longest tracks first. Only used for fallback at the end
930  std::sort(xTrks.rbegin(), xTrks.rend(), CompareByLength);
931  std::sort(yTrks.rbegin(), yTrks.rend(), CompareByLength);
932 
933  double bestVx, bestVy, bestVz1, bestVz2;
934  // I don't believe there's a path that doesn't set these, but need this to
935  // shut the compiler up in optimized builds.
936  bestVx = bestVy = bestVz1 = bestVz2 = 0;
937 
938  for(geo::View_t view: {geo::kX, geo::kY}){
939  double bestScore = 1e10;
940  const std::vector<rb::Track>& trks = (view == geo::kX) ? xTrks : yTrks;
941  for(unsigned int idx1 = 0; idx1 < trks.size()-1; ++idx1){
942  if(trks[idx1].ExtentPlane() < 3) continue;
943  for(bool end1: {false, true}){
944  for(unsigned int idx2 = idx1+1; idx2 < trks.size(); ++idx2){
945  if(trks[idx2].ExtentPlane() < 3) continue;
946  for(bool end2: {false, true}){
947 
948  // Start and direction defining point of these two lines
949  const TVector3 s1 = end1 ? trks[idx1].Stop() : trks[idx1].Start();
950  const TVector3 s2 = end2 ? trks[idx2].Stop() : trks[idx2].Start();
951  const TVector3 d1 = trks[idx1].TrajectoryPoint(end1 ? trks[idx1].NTrajectoryPoints()-2 : 1);
952  const TVector3 d2 = trks[idx2].TrajectoryPoint(end2 ? trks[idx2].NTrajectoryPoints()-2 : 1);
953 
954  double vxy, vz;
955  const bool ok = geo::LineIntersection(s1[view], s1.Z(), d1[view], d1.Z(),
956  s2[view], s2.Z(), d2[view], d2.Z(),
957  vxy, vz);
958  if(!ok) continue;
959 
960  const double dist1 = util::pythag(s1[view]-vxy, s1.Z()-vz);
961  const double dist2 = util::pythag(s2[view]-vxy, s2.Z()-vz);
962 
963  // Favour both long tracks and hits close to the ends
964  const double score = dist1/trks[idx1].TotalLength()+dist2/trks[idx2].TotalLength();
965  if(score < bestScore){
966  bestScore = score;
967  if(view == geo::kX){
968  bestVx = vxy;
969  bestVz1 = vz;
970  }
971  else{
972  bestVy = vxy;
973  bestVz2 = vz;
974  }
975  } // end if better
976  } // end for end2
977  } // end for idx2
978  } // end for end1
979  } // end for idx1
980  // Too far from the ends, fallback
981  if(bestScore > 6){
982  if(view == geo::kX){
983  bestVx = trks[0].Start().X();
984  bestVz1 = trks[0].Start().Z();
985  }
986  else{
987  bestVy = trks[0].Start().Y();
988  bestVz2 = trks[0].Start().Z();
989  }
990  }
991  } // end for view
992 
993  // Time is an afterthought
994  return new rb::Vertex(bestVx, bestVy, (bestVz1+bestVz2)/2, xTrks[0].MinTNS());
995  }
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
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
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
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 266 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().

267  {
270 
271  // get the total number of non masked off channels by looping over all
272  // planes/cells in the detector
273  int numGoodChannels = 0;
274 
275  for(unsigned int iPlane = 0; iPlane < geom->NPlanes(); ++iPlane){
276  const geo::PlaneGeo* p = geom->Plane(iPlane);
277  for(unsigned int iCell = 0; iCell < p->Ncells(); ++iCell){
278  if(!bcl->IsBad(iPlane, iCell)) ++numGoodChannels;
279  }
280  }
281 
282  return clust->NCell() > 0.15*numGoodChannels;
283  }
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 = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
void dt::DiscreteTracker::MergeParasiteTracks ( std::vector< rb::Track > &  tracks) const
protected

Definition at line 998 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, getGoodRuns4SAM::n, rb::Cluster::NCell(), and rb::CellHit::Plane().

Referenced by RecoDiscrete().

999  {
1000  // Worried this could be too aggressive around the vertex
1001 
1002  const unsigned int N = tracks.size();
1003  if(N < 2) return;
1004 
1005  // Longest first
1006  std::sort(tracks.begin(), tracks.end(),
1007  [](const rb::Track& a, const rb::Track& b){
1008  return a.ExtentPlane() > b.ExtentPlane();
1009  });
1010 
1011  for(unsigned int hostIdx = 0; hostIdx < N-1; ++hostIdx){
1012  const rb::HitMap hmap(&tracks[hostIdx]);
1013  for(unsigned int paraIdx = hostIdx+1; paraIdx < N; ++paraIdx){
1014  const rb::Track& para = tracks[paraIdx];
1015  bool match = true;
1016  for(unsigned int i = 0; i < para.NCell(); ++i){
1017  const art::Ptr<rb::CellHit> chit = para.Cell(i);
1018  if(!hmap.CellExists(chit->Plane(), chit->Cell()-1) &&
1019  !hmap.CellExists(chit->Plane(), chit->Cell()+1)){
1020  match = false;
1021  break;
1022  }
1023  } // end for i
1024  if(match){
1025  // Make sure there are no duplicates in the final host track
1026  rb::HitMap hmap(tracks[hostIdx].AllCells());
1027  for(unsigned int n = 0; n < para.NCell(); ++n){
1028  const art::Ptr<rb::CellHit> chit = para.Cell(n);
1029  if(!hmap.CellExists(chit->Plane(), chit->Cell()))
1030  tracks[hostIdx].Add(chit);
1031  }
1032  tracks.erase(tracks.begin()+paraIdx);
1033  MergeParasiteTracks(tracks);
1034  return;
1035  }
1036  } // end for paraIdx
1037  } // end for hostIdx
1038  }
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
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
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
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 287 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().

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

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

201  {
202  std::unique_ptr<std::vector<rb::Track> > debug_seeds_col(new std::vector<rb::Track>);
203  std::unique_ptr<std::vector<rb::Track> > debug_unspliced_col(new std::vector<rb::Track>);
204  std::unique_ptr<std::vector<rb::Track> > debug_chains_col(new std::vector<rb::Track>);
205  std::unique_ptr<std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
206  std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
207  // Associate tracks with the slice they came from
208  std::unique_ptr<art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
209 
211  evt.getByLabel(fClusterInput, clusts);
212 
213  const int clustMax = clusts->size();
214  for(int clustIdx = 0; clustIdx < clustMax; ++clustIdx){
215  art::Ptr<rb::Cluster> clust(clusts, clustIdx);
216  if(fObeyPreselection && rb::IsFiltered(evt, clust)) continue;
217 
218  if(clust->IsNoise()) continue;
219 
220  if(HighActivity(clust)){
221  mf::LogWarning("DiscreteTracker") << "Skipping event due to very high activity";
222  continue;
223  }
224 
225  std::vector<rb::Track> tracks[2];
226  for(geo::View_t view: {geo::kX, geo::kY}){
227  const char viewC = (view == geo::kX) ? 'x' : 'y';
228 
230  const int Nx = clust->NCell(view);
231  for(int n = 0; n < Nx; ++n) cells.push_back(clust->Cell(view, n));
232 
233  fprintf(fDebugFile, "%d: %c: ", evt.event(), viewC);
234 
235  tracks[view] = RecoDiscrete(cells,
236  debug_seeds_col.get(),
237  debug_unspliced_col.get(),
238  debug_chains_col.get());
239 
240  for(unsigned int n = 0; n < tracks[view].size(); ++n)
241  fprintf(fDebugFile, "%u ", tracks[view][n].NCell());
242  fprintf(fDebugFile, "\n");
243 
244  for(std::vector<rb::Track>::const_iterator it = tracks[view].begin();
245  it != tracks[view].end(); ++it){
246  trackcol->push_back(*it);
247  // Associate most recently added track with the slice
248  util::CreateAssn(*this, evt, *trackcol, clust, *assns);
249  } // end for it
250  } // end for view
251 
252  const rb::Vertex* vtx = FindVertex(tracks[geo::kX], tracks[geo::kY]);
253  if(vtx) vtxcol->push_back(*vtx);
254  delete vtx;
255  } // end for clustIdx
256 
257  evt.put(std::move(debug_seeds_col), "seeds");
258  evt.put(std::move(debug_unspliced_col), "unspliced");
259  evt.put(std::move(debug_chains_col), "chains");
260  evt.put(std::move(trackcol));
261  evt.put(std::move(vtxcol));
262  evt.put(std::move(assns));
263  }
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.
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
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
bool HighActivity(const art::Ptr< rb::Cluster > &clust) const
Cut to prevent choking on large air shower events.
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
EventNumber_t event() const
Definition: Event.h:67
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
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 562 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(), getGoodRuns4SAM::n, dt::Cand::NSegs(), RecoOneTrack(), runNovaSAM::ret, art::PtrVector< T >::size(), SpliceChains(), dt::Chain::ToTrack(), dt::Chain::TrimEnds(), and registry_explorer::v.

Referenced by produce().

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

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

Referenced by DiscreteTracker().

172  {
173  fClusterInput = pset.get<std::string>("ClusterInput");
174  fConsiderBadChans = pset.get<bool>("ConsiderBadChans");
175  fMinDotProduct = pset.get<double>("MinDotProduct");
176  fVerticalTrackMinCells = pset.get<int>("VerticalTrackMinCells");
177  fOnlyOneTrack = pset.get<bool>("OnlyOneTrack");
178  fObeyPreselection = pset.get<bool>("ObeyPreselection");
179 
181  fDebugFile = fopen(pset.get<std::string>("DebugFile").c_str(), "w");
182  }
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 544 of file DiscreteTracker_module.cc.

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

Referenced by RecoDiscrete().

548  {
549  const Cand& seed = BestSeed(v, ok);
550  if(!ok) return Chain();
551  seeds_out->push_back(seed.ToTrack());
552 
553  Chain chain(seed);
554  ExtendChain(v, chunk_map, chain, kUpstream);
555  ExtendChain(v, chunk_map, chain, kDownstream);
556  ok = true;
557  return chain;
558  }
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::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void dt::DiscreteTracker::SpliceChains ( std::vector< Chain > &  chains) const
protected

Definition at line 671 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(), getGoodRuns4SAM::n, dt::Chain::NHitPlanes(), dt::OtherDir(), file_size_ana::parent, dt::SpliceDesc::parentIdx, sd(), dt::SpliceDesc::spliced, SpliceIsSane(), TryExtendChainOne(), and registry_explorer::v.

Referenced by RecoDiscrete().

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

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

Referenced by SpliceChains().

634  {
635  // Must not overlap. Can't jump more than half the whole total length
636 
637  if(a.LastPlane() < b.FirstPlane()){
638  if(DotScore(a.LastCand(), b.FirstCand(), true) < fMinDotProduct) return false;
639 
640  return (2*(b.FirstPlane()-a.LastPlane()) < b.LastPlane()-a.FirstPlane());
641  }
642  if(b.LastPlane() < a.FirstPlane()){
643  if(DotScore(b.LastCand(), a.FirstCand(), true) < fMinDotProduct) return false;
644 
645  return (2*(a.FirstPlane()-b.LastPlane()) < a.LastPlane()-b.FirstPlane());
646  }
647 
648  return false;
649  }
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 462 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().

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

407  {
408  // TODO: hoist this out of while loop below
409  std::vector<Chunk> chunks;
410  const std::vector<Segment> segs = PossibleNextSegs(v, chunk_map,
411  cand, dir, chunks, shw,
412  exactHits);
413  if(segs.empty()) return false;
414 
415  // Backup of cand, so we can restore it in case TryAddSeg succeeds
416  Cand bak = cand;
417 
418  std::vector<Cand> extended;
419 
420  const unsigned int segMax = segs.size();
421  for(unsigned int segIdx = 0; segIdx < segMax; ++segIdx){
422  const Segment& seg = segs[segIdx];
423 
424  if(cand.TryAddSeg(seg, dir)){
425 
426  // May need to add the chunk to match
427  if((dir == kUpstream && !seg.left) ||
428  (dir == kDownstream && seg.left)){
429  cand.AddChunk(chunks[segIdx], dir);
430  }
431 
432  const bool smallOverlap = (bak.NSegs() <= 2);
433  if(true){//smallOverlap){
434  // If there's only a small overlap, ensure dot product is sensible
435  const double dotscore = DotScore(prev_dir, cand);
436  fAngle->Fill(acos(dotscore));
437  if(smallOverlap) fSinglesAngle->Fill(acos(dotscore));
438  if(smallOverlap && dotscore < fMinDotProduct){
439  cand = bak;
440  continue;
441  }
442  }
443 
444  // Otherwise it's good
445  extended.push_back(cand);
446  cand = bak;
447  } // end if added seg
448  } // end for segIdx
449 
450  if(extended.empty()) return false;
451 
452  // TODO: Arbitrary for now
453  cand = extended.front();
454  if(extended.size() > 1){
455  // std::cout << "***WARNING: MORE THAN ONE EXTENSION CANDIDATE (" << extended.size() << "). SHOULD HAVE TRIED TO MAKE A SMARTER DECISION***" << std::endl;
456  }
457 
458  return true;
459  }
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 511 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().

517  {
518  while(true){
519  if(cur.NSegs() < 2) return false;
520  if(TryExtendCand(cur, prev, v, chunk_map, dir, exactHits)) return true;
521 
522  cur.ShortenOne(dir);
523 
524  // Too much freedom
525  if(cur.IsAllShower()) return false;
526  if(cur.NChunks() == 0) return false; // all dead
527  if(cur.AllHits().empty()) return false;
528  }
529  }
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
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

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: