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

Match 2D tracks into 3D tracks. More...

Inheritance diagram for dt::ViewMerger:
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

 ViewMerger (const fhicl::ParameterSet &pset)
 
 ~ViewMerger ()
 
virtual void reconfigure (const fhicl::ParameterSet &pset)
 
virtual void beginJob ()
 
virtual void endJob ()
 
virtual void produce (art::Event &evt)
 
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

std::vector< rb::TrackMatchTracks (const rb::HitMap &sliceMap, std::list< rb::Track > &xtracks, std::list< rb::Track > &ytracks) const
 Do most of the work. Called from produce. More...
 
std::vector< rb::TrackMatchVerticals (std::list< rb::Track > &xtracks, std::list< rb::Track > &ytracks) const
 Find vertical tracks seperated by one plane. Return results of merge. More...
 
rb::Track ZipVerticalTracks (const rb::Track &xtrk, const rb::Track &ytrk) const
 Tilt tracks arbitrarily and merge them into a 3D track. More...
 
void TotalChargePerView (const rb::Track &zip, double &totx, double &toty) const
 
double MissingChargeByExtension (const rb::HitMap &sliceMap, const rb::Track &aTrk, const rb::Track &bTrk) const
 
double ScoreForExtension (const rb::HitMap &sliceMap, const rb::Track &xTrk, const rb::Track &yTrk) const
 
double ScoreForJoinPlusExtension (const rb::HitMap &sliceMap, const rb::Track &single, const rb::Track &parta, const rb::Track &partb) const
 
rb::Track JoinTracks (const rb::Track &parta, const rb::Track &partb) const
 Append two tracks in the same view (crudely). More...
 
int CountUnexplainedOnLine (const rb::HitMap &hmap, TVector3 start, TVector3 end, geo::View_t view) const
 hmap should contain cells that are hit and so OK to pass through More...
 
double ChargeDifferenceBetweenViews (const rb::Track &zip) const
 
void MaybeExtendToVertex (const rb::Vertex *vtx, std::vector< rb::Track > &tracks) const
 
double MatchScore (double unexp, const rb::Track &aTrk, const rb::Track &bTrk) const
 Find number of overlapping planes, return unexp/overlap. More...
 
rb::Track ZipTracks (const rb::Track &aTrk, const rb::Track &bTrk) const
 Make a 3D track based on two 2D tracks. More...
 
void Extremes (const rb::Track &trk, std::map< int, int > &mins, std::map< int, int > &maxs) const
 
bool Adjacent (const rb::Track &a, const rb::Track &b) const
 
void MergeAdjacentTracks (std::list< rb::Track > &trks) const
 
rb::VertexFindVertex (std::vector< rb::Track > trks) const
 Toy vertex finder. More...
 
rb::Track FakeThirdDimension (const rb::Track &trk, const std::vector< rb::Track > &merged) const
 Turn a 2D track into a 3D track. More...
 
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 fSliceInput
 Module name to look for slices under. More...
 
std::string fTrackInput
 Module name to look for 2D tracks under. More...
 
bool fObeyPreselection
 Check rb::IsFiltered? More...
 
FILE * fDebugFile
 

Detailed Description

Match 2D tracks into 3D tracks.

Definition at line 38 of file ViewMerger_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::ViewMerger::ViewMerger ( const fhicl::ParameterSet pset)
explicit

Definition at line 151 of file ViewMerger_module.cc.

References reconfigure().

152  : fDebugFile(0)
153  {
154  produces<std::vector<rb::Track> >();
155  produces<std::vector<rb::Vertex> >();
156  // Tracks are associated with the slice they came from
157  produces<art::Assns<rb::Track, rb::Cluster> >();
158 
159  reconfigure(pset);
160  }
virtual void reconfigure(const fhicl::ParameterSet &pset)
dt::ViewMerger::~ViewMerger ( )

Definition at line 163 of file ViewMerger_module.cc.

164  {
165  }

Member Function Documentation

bool dt::ViewMerger::Adjacent ( const rb::Track a,
const rb::Track b 
) const
protected

Helper for MergeAdjacentTracks. Are all hits on one track adjacent to the other one?

Definition at line 816 of file ViewMerger_module.cc.

References ana::assert(), rb::Cluster::ExtentPlane(), Extremes(), rb::Cluster::Is2D(), it, and rb::Cluster::View().

Referenced by MergeAdjacentTracks().

817  {
818  assert(a.View() == b.View());
819  assert(a.Is2D());
820 
821  if(a.ExtentPlane() < b.ExtentPlane()) return Adjacent(b, a);
822  // Guaranteed b is the shorter one
823 
824  // Don't attach verticals
825  if(b.ExtentPlane() == 1) return false;
826 
827  std::map<int, int> amins, amaxs, bmins, bmaxs;
828  Extremes(a, amins, amaxs);
829  Extremes(b, bmins, bmaxs);
830 
831  // Maybe a is all below b
832  bool all = true;
833  for(std::map<int, int>::iterator it = bmins.begin(); it != bmins.end(); ++it){
834  if(amaxs[it->first] != it->second-1){
835  all = false;
836  break;
837  }
838  }
839  if(all) return true;
840 
841  // Or maybe a is all above b
842  for(std::map<int, int>::iterator it = bmaxs.begin(); it != bmaxs.end(); ++it){
843  if(amins[it->first] != it->second+1) return false;
844  }
845  return true;
846  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
set< int >::iterator it
bool Is2D() const
Definition: Cluster.h:96
void Extremes(const rb::Track &trk, std::map< int, int > &mins, std::map< int, int > &maxs) const
bool Adjacent(const rb::Track &a, const rb::Track &b) const
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
assert(nhit_max >=nhit_nbins)
void dt::ViewMerger::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 180 of file ViewMerger_module.cc.

181  {
182  }
double dt::ViewMerger::ChargeDifferenceBetweenViews ( const rb::Track zip) const
protected

Definition at line 681 of file ViewMerger_module.cc.

References stan::math::fabs(), and TotalChargePerView().

Referenced by ScoreForJoinPlusExtension().

682  {
683  double totx, toty;
684  TotalChargePerView(zip, totx, toty);
685 
686  return fabs(totx-toty);
687  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void TotalChargePerView(const rb::Track &zip, double &totx, double &toty) const
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
int dt::ViewMerger::CountUnexplainedOnLine ( const rb::HitMap hmap,
TVector3  start,
TVector3  end,
geo::View_t  view 
) const
protected

hmap should contain cells that are hit and so OK to pass through

Definition at line 642 of file ViewMerger_module.cc.

References getBrightness::cell, rb::HitMap::CellExists(), geo::GeometryBase::CountCellsOnLineFast(), geom(), chaninfo::BadChanList::IsBad(), chaninfo::BadChanList::IsLowEfficiency(), geo::kX, geo::kY, getGoodRuns4SAM::n, NDAPDHVSetting::plane, and runNovaSAM::ret.

Referenced by MaybeExtendToVertex(), MissingChargeByExtension(), and ScoreForJoinPlusExtension().

645  {
646  // TODO: we need to penalize the case where the track now runs outside the
647  // detector for some of its length. This is actually leading to
648  // mismerges...
651 
652  std::vector<geo::OfflineChan> xcells, ycells;
653  geom->CountCellsOnLineFast(start.X(), start.Y(), start.Z(),
654  end.X(), end.Y(), end.Z(),
655  xcells, ycells);
656 
657  int ret = 0;
658  if(view != geo::kY){ // X or 3D
659  for(unsigned int n = 0; n < xcells.size(); ++n){
660  const int plane = xcells[n].Plane();
661  const int cell = xcells[n].Cell();
662  if(!hmap.CellExists(plane, cell) &&
663  !bcl->IsBad(plane, cell) &&
664  !bcl->IsLowEfficiency(plane, cell)) ++ret;
665  }
666  }
667  if(view != geo::kX){ // Y or 3D
668  for(unsigned int n = 0; n < ycells.size(); ++n){
669  const int plane = ycells[n].Plane();
670  const int cell = ycells[n].Cell();
671  if(!hmap.CellExists(plane, cell) &&
672  !bcl->IsBad(plane, cell) &&
673  !bcl->IsLowEfficiency(plane, cell)) ++ret;
674  }
675  }
676 
677  return ret;
678  }
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CountCellsOnLineFast(double x1, double y1, double z1, double x2, double y2, double z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
Make list of cells in each view traversed by a line segment.
void geom(int which=0)
Definition: geom.C:163
bool CellExists(unsigned int planeIdx, unsigned int cellIdx) const
Does the map contain any cell at this position?
Definition: HitMap.cxx:81
bool IsBad(int plane, int cell)
bool IsLowEfficiency(int plane, int cell)
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
void dt::ViewMerger::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 185 of file ViewMerger_module.cc.

References fclose(), and fDebugFile.

186  {
188  }
fclose(fg1)
void dt::ViewMerger::Extremes ( const rb::Track trk,
std::map< int, int > &  mins,
std::map< int, int > &  maxs 
) const
protected

Helper for Adjacent. Fills mins and maxs with minimum and maximum cell number in each plane of trk

Definition at line 801 of file ViewMerger_module.cc.

References plot_validation_datamc::c, rb::CellHit::Cell(), rb::Cluster::Cell(), getGoodRuns4SAM::n, rb::Cluster::NCell(), and rb::CellHit::Plane().

Referenced by Adjacent().

804  {
805  for(unsigned int n = 0; n < trk.NCell(); ++n){
806  const int p = trk.Cell(n)->Plane();
807  const int c = trk.Cell(n)->Cell();
808 
809  if(c > maxs[p]) maxs[p] = c;
810  // Careful with zero default value
811  if(mins.find(p) == mins.end() || c < mins[p]) mins[p] = c;
812  }
813  }
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
const char * p
Definition: xmltok.h:285
unsigned short Cell() const
Definition: CellHit.h:40
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
rb::Track dt::ViewMerger::FakeThirdDimension ( const rb::Track trk,
const std::vector< rb::Track > &  merged 
) const
protected

Turn a 2D track into a 3D track.

Uses merged to guess a W position

Definition at line 942 of file ViewMerger_module.cc.

References rb::Cluster::AllCells(), ana::assert(), DEFINE_ART_MODULE(), rb::Prong::Dir(), rb::Cluster::ExtentPlane(), rb::Track::InterpolateXY(), rb::Prong::Is3D(), geo::kX, getGoodRuns4SAM::n, rb::Track::NTrajectoryPoints(), gen_flatrecord::pt, runNovaSAM::ret, febshutoff_auto::start, rb::Prong::Start(), rb::Track::TrajectoryPoint(), POTSpillRate::view, rb::Cluster::View(), w, submit_syst::x, submit_syst::y, and test::z.

944  {
945  assert(!trk.Is3D());
946 
947  const double z = trk.Start().Z();
948 
949  const rb::Track* best = 0;
950 
951  // The best source for W information is the longest overlapping track
952  for(unsigned int n = 0; n < merged.size(); ++n){
953  if(merged[n].Start().Z()-5 < z &&
954  merged[n].Stop().Z()+5 > z){
955  if(!best || merged[n].ExtentPlane() > best->ExtentPlane())
956  best = &merged[n];
957  }
958  }
959 
960  // Failing that, just the longest track
961  if(!best){
962  for(unsigned int n = 0; n < merged.size(); ++n){
963  if(!best || merged[n].ExtentPlane() > best->ExtentPlane())
964  best = &merged[n];
965  }
966  }
967 
968  const geo::View_t view = trk.View();
969 
970  double w = 0;
971  if(best){
972  double x, y;
973  best->InterpolateXY(z, x, y);
974  w = (view == geo::kX) ? y : x;
975  }
976 
977  // Fill in W information everywhere
978 
979  TVector3 start = trk.Start();
980  start[1-view] = w;
981 
982  rb::Track ret(trk.AllCells(), start, trk.Dir());
983 
984  for(unsigned int n = 1; n < trk.NTrajectoryPoints(); ++n){
985  TVector3 pt = trk.TrajectoryPoint(n);
986  pt[1-view] = w;
987  ret.AppendTrajectoryPoint(pt);
988  }
989 
990  return ret;
991  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
size_t NTrajectoryPoints() const
Definition: Track.h:83
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
virtual TVector3 Start() const
Definition: Prong.h:73
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
z
Definition: test.py:28
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
assert(nhit_max >=nhit_nbins)
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
Float_t w
Definition: plot.C:20
virtual bool Is3D() const
Definition: Prong.h:71
rb::Vertex * dt::ViewMerger::FindVertex ( std::vector< rb::Track trks) const
protected

Toy vertex finder.

Uses the intersection of the 3D tracks whose intersection is proportionally nearest their ends.

Definition at line 876 of file ViewMerger_module.cc.

References dt::CompareByLengthVM(), geo::kXorY, geo::LineIntersection(), and util::pythag().

Referenced by produce().

877  {
878  if(trks.empty()) return 0;
879 
880  // Longest tracks first. Only used for fallback at the end
881  std::sort(trks.rbegin(), trks.rend(), CompareByLengthVM);
882 
883  double bestVx, bestVy, bestVz;
884  // I don't believe there's a path that doesn't set these, but need this to
885  // shut the compiler up in optimized builds.
886  bestVx = bestVy = bestVz = 0;
887 
888  double bestScore = 1e10;
889  for(unsigned int idx1 = 0; idx1 < trks.size()-1; ++idx1){
890  if(trks[idx1].View() != geo::kXorY) continue;
891  if(trks[idx1].ExtentPlane() < 3) continue;
892  for(bool end1: {false, true}){
893  for(unsigned int idx2 = idx1+1; idx2 < trks.size(); ++idx2){
894  if(trks[idx2].View() != geo::kXorY) continue;
895  if(trks[idx2].ExtentPlane() < 3) continue;
896  for(bool end2: {false, true}){
897 
898  // Start and direction defining point of these two lines
899  const TVector3 s1 = end1 ? trks[idx1].Stop() : trks[idx1].Start();
900  const TVector3 s2 = end2 ? trks[idx2].Stop() : trks[idx2].Start();
901  const TVector3 d1 = trks[idx1].TrajectoryPoint(end1 ? trks[idx1].NTrajectoryPoints()-2 : 1);
902  const TVector3 d2 = trks[idx2].TrajectoryPoint(end2 ? trks[idx2].NTrajectoryPoints()-2 : 1);
903 
904  double vx, vy, vz1, vz2;
905  bool ok = geo::LineIntersection(s1.X(), s1.Z(), d1.X(), d1.Z(),
906  s2.X(), s2.Z(), d2.X(), d2.Z(),
907  vx, vz1);
908  if(!ok) continue;
909  ok = geo::LineIntersection(s1.Y(), s1.Z(), d1.Y(), d1.Z(),
910  s2.Y(), s2.Z(), d2.Y(), d2.Z(),
911  vy, vz2);
912  if(!ok) continue;
913  const double vz = (vz1+vz2)/2;
914 
915  const double dist1 = util::pythag(s1.X()-vx, s1.Y()-vy, s1.Z()-vz);
916  const double dist2 = util::pythag(s2.X()-vx, s2.Y()-vy, s2.Z()-vz);
917 
918  // Favour both long tracks and hits close to the ends
919  const double score = dist1/trks[idx1].TotalLength()+dist2/trks[idx2].TotalLength();
920  if(score < bestScore){
921  bestScore = score;
922  bestVx = vx;
923  bestVy = vy;
924  bestVz = vz;
925  } // end if better
926  } // end for end2
927  } // end for idx2
928  } // end for end1
929  } // end for idx1
930  // Too far from the ends, fallback
931  if(bestScore > 6){
932  bestVx = trks[0].Start().X();
933  bestVy = trks[0].Start().Y();
934  bestVz = trks[0].Start().Z();
935  }
936 
937  // Time is an afterthought
938  return new rb::Vertex(bestVx, bestVy, bestVz, trks[0].MinTNS());
939  }
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
X or Y views.
Definition: PlaneGeo.h:30
Definition: View.py:1
bool CompareByLengthVM(const rb::Track &a, const rb::Track &b)
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  }
rb::Track dt::ViewMerger::JoinTracks ( const rb::Track parta,
const rb::Track partb 
) const
protected

Append two tracks in the same view (crudely).

Definition at line 619 of file ViewMerger_module.cc.

References rb::Cluster::AllCells(), ana::assert(), rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), getGoodRuns4SAM::n, rb::Track::NTrajectoryPoints(), rb::Track::TrajectoryPoint(), and rb::Cluster::View().

Referenced by MatchTracks(), and ScoreForJoinPlusExtension().

621  {
622  assert(parta.View() == partb.View());
623 
624  // Must be in this order, must not overlap
625  assert(partb.MinPlane() >= parta.MaxPlane());
626 
627  // This is pretty crude. Could do a neater job with tighter interface into
628  // the 2D tracker.
629 
630  rb::Track join(parta);
631  // Right now, no two tracks can share hits, so this is safe. Might not be
632  // so in future, would have to remove duplicates.
633  join.Add(partb.AllCells());
634  for(unsigned int n = 0; n < partb.NTrajectoryPoints(); ++n){
635  join.AppendTrajectoryPoint(partb.TrajectoryPoint(n));
636  }
637 
638  return join;
639  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
size_t NTrajectoryPoints() const
Definition: Track.h:83
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
double dt::ViewMerger::MatchScore ( double  unexp,
const rb::Track aTrk,
const rb::Track bTrk 
) const
protected

Find number of overlapping planes, return unexp/overlap.

Definition at line 721 of file ViewMerger_module.cc.

References abs(), ana::assert(), e, rb::Cluster::ExtentPlane(), rb::Cluster::MaxPlane(), and rb::Cluster::MinPlane().

Referenced by ScoreForExtension(), and ScoreForJoinPlusExtension().

724  {
725  // These casts to int are absolutely necessary :-@
726  const int dStart = abs(int(aTrk.MinPlane())-int(bTrk.MinPlane()));
727  assert(dStart > 0);
728  const int dEnd = abs(int(aTrk.MaxPlane())-int(bTrk.MaxPlane()));
729  assert(dEnd > 0);
730  const int totL = aTrk.ExtentPlane()+bTrk.ExtentPlane();
731  const int overlap = totL-dStart-dEnd;
732  if(overlap <= 0) return 1e12;
733 
734  // Add a small amount to miss so we can distinguish different zeroes
735  return (miss+1e-3)/overlap;
736  }
void abs(TH1 *hist)
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
Float_t e
Definition: plot.C:35
std::vector< rb::Track > dt::ViewMerger::MatchTracks ( const rb::HitMap sliceMap,
std::list< rb::Track > &  xtracks,
std::list< rb::Track > &  ytracks 
) const
protected

Do most of the work. Called from produce.

Parameters
sliceMapAll the hits in this slice
xtracksThe 2D tracks from the x-view to merge. Those used will be removed from the list. Any left unmerged will remain in the list.
ytracksDitto
Returns
The 3D tracks produced

Definition at line 296 of file ViewMerger_module.cc.

References ana::assert(), allTimeWatchdog::endl, JoinTracks(), geo::kX, geo::kY, LOG_DEBUG, MatchVerticals(), submit_nova_art::mode, ScoreForExtension(), ScoreForJoinPlusExtension(), and ZipTracks().

Referenced by produce().

299  {
300  // This is what we return
301  std::vector<rb::Track> matchedTracks;
302 
303  typedef std::list<rb::Track>::iterator it_t;
304 
305  // Keep going until we stop being able to match anything
306  while(true){
307  double bestScore = 1e10;
308  it_t bestX = xtracks.end();
309  it_t bestY = ytracks.end();
310  it_t bestJoin = xtracks.end(); // Picking this end() is just convention
311 
312  for(it_t xTrk = xtracks.begin(); xTrk != xtracks.end(); ++xTrk){
313  assert(xTrk->View() == geo::kX);
314  // No way to sensibly extend the trajectory for these right now
315  if(xTrk->ExtentPlane() == 1) continue;
316 
317  for(it_t yTrk = ytracks.begin(); yTrk != ytracks.end(); ++yTrk){
318  assert(yTrk->View() == geo::kY);
319  if(yTrk->ExtentPlane() == 1) continue;
320 
321  // Simple merge, join in x, join in y
322  for(int mode = 0; mode < 3; ++mode){
323  // Don't try to join tracks within a view. DiscreteTracker got
324  // better at this, and right now we're just making a mess when we
325  // try.
326  if(mode == 1 || mode == 2) continue;
327 
328  it_t joinTrk;
329  if(mode == 0) joinTrk = xtracks.end();
330  if(mode == 1){joinTrk = xTrk; ++joinTrk;}
331  if(mode == 2){joinTrk = yTrk; ++joinTrk;}
332 
333  // For simple mode, go once with joinTrk invalid. For other modes,
334  // loop through all joinable tracks in the right view.
335  for(; mode == 0 ||
336  (mode == 1 && joinTrk != xtracks.end()) ||
337  (mode == 2 && joinTrk != ytracks.end()); ++joinTrk){
338  // This could be done smarter if we interfaced with
339  // DiscreteTracker properly. It has a lot more information about
340  // these tracks that we could use.
341  double score;
342  if(mode == 0) score = ScoreForExtension(sliceMap, *xTrk, *yTrk);
343  if(mode == 1) score = ScoreForJoinPlusExtension(sliceMap, *yTrk,
344  *xTrk, *joinTrk);
345  if(mode == 2) score = ScoreForJoinPlusExtension(sliceMap, *xTrk,
346  *yTrk, *joinTrk);
347 
348  if(score < bestScore){
349  bestScore = score;
350  bestX = xTrk;
351  bestY = yTrk;
352  bestJoin = joinTrk;
353  }
354 
355  if(mode == 0) break;
356  } // end for joinTrk
357  } // end for mode
358  } // end for yTrk
359  } // end for xTrk
360 
361  // No matches left
362  if(bestX == xtracks.end()) break;
363 
364  LOG_DEBUG("ViewMerger") << "Best score: " << bestScore << std::endl;
365 
366  if(bestJoin != xtracks.end()){
367  if(bestJoin->View() == geo::kX)
368  matchedTracks.push_back(ZipTracks(JoinTracks(*bestX, *bestJoin),
369  *bestY));
370  else
371  matchedTracks.push_back(ZipTracks(*bestX,
372  JoinTracks(*bestY, *bestJoin)));
373  }
374  else{
375  matchedTracks.push_back(ZipTracks(*bestX, *bestY));
376  }
377 
378  // Can't use either 2D track again
379  xtracks.erase(bestX);
380  ytracks.erase(bestY);
381  if(bestJoin != xtracks.end()){
382  if(bestJoin->View() == geo::kX)
383  xtracks.erase(bestJoin);
384  else
385  ytracks.erase(bestJoin);
386  }
387  } // end while
388 
389  return matchedTracks;
390  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
rb::Track ZipTracks(const rb::Track &aTrk, const rb::Track &bTrk) const
Make a 3D track based on two 2D tracks.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double ScoreForExtension(const rb::HitMap &sliceMap, const rb::Track &xTrk, const rb::Track &yTrk) const
double ScoreForJoinPlusExtension(const rb::HitMap &sliceMap, const rb::Track &single, const rb::Track &parta, const rb::Track &partb) const
assert(nhit_max >=nhit_nbins)
rb::Track JoinTracks(const rb::Track &parta, const rb::Track &partb) const
Append two tracks in the same view (crudely).
std::vector< rb::Track > dt::ViewMerger::MatchVerticals ( std::list< rb::Track > &  xtracks,
std::list< rb::Track > &  ytracks 
) const
protected

Find vertical tracks seperated by one plane. Return results of merge.

Definition at line 394 of file ViewMerger_module.cc.

References std::abs(), ana::assert(), geo::kX, geo::kY, file_size_ana::progress, runNovaSAM::ret, and ZipVerticalTracks().

Referenced by MatchTracks(), and produce().

396  {
397  std::vector<rb::Track> ret;
398 
399  typedef std::list<rb::Track>::iterator it_t;
400 
401  bool progress = true;
402  while(progress){
403  progress = false;
404 
405  for(it_t xTrk = xtracks.begin(); xTrk != xtracks.end(); ++xTrk){
406  assert(xTrk->View() == geo::kX);
407  if(xTrk->ExtentPlane() != 1) continue;
408 
409  for(it_t yTrk = ytracks.begin(); yTrk != ytracks.end(); ++yTrk){
410  assert(yTrk->View() == geo::kY);
411  if(yTrk->ExtentPlane() != 1) continue;
412 
413  if(std::abs(1.*xTrk->MinPlane() - 1.*yTrk->MinPlane()) != 1.) continue;
414 
415  ret.push_back(ZipVerticalTracks(*xTrk, *yTrk));
416  xtracks.erase(xTrk);
417  ytracks.erase(yTrk);
418  progress = true;
419  break;
420  } // end for yTrk
421  if(progress) break;
422  } // end for xTrk
423  } // end while
424 
425  return ret;
426  }
rb::Track ZipVerticalTracks(const rb::Track &xtrk, const rb::Track &ytrk) const
Tilt tracks arbitrarily and merge them into a 3D track.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
bool progress
Insert implicit nodes #####.
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
assert(nhit_max >=nhit_nbins)
void dt::ViewMerger::MaybeExtendToVertex ( const rb::Vertex vtx,
std::vector< rb::Track > &  tracks 
) const
protected

Extend each track to the vertex point, only if it's completely uncontroversial.

Definition at line 690 of file ViewMerger_module.cc.

References rb::HitMap::Add(), CountUnexplainedOnLine(), Dot(), rb::Vertex::GetXYZ(), getGoodRuns4SAM::n, and Unit().

Referenced by produce().

692  {
693  rb::HitMap hmap;
694  for(unsigned int n = 0; n < tracks.size(); ++n) hmap.Add(&tracks[n]);
695 
696  for(unsigned int n = 0; n < tracks.size(); ++n){
697  // Only if we don't pass through a single non-dead cell
698  if(CountUnexplainedOnLine(hmap, tracks[n].Start(), vtx->GetXYZ(),
699  tracks[n].View()) == 0){
700 
701  TVector3 xyz = vtx->GetXYZ();
702  if(!tracks[n].Is3D()) xyz[1-tracks[n].View()] = 0;
703 
704  // And only if the dot-product is reasonable
705  if((tracks[n].Start()-xyz).Unit().Dot(tracks[n].Dir()) < .94) continue;
706 
707  tracks[n].PrependTrajectoryPoint(xyz);
708  if(tracks[n].Is3D()){
709  // Overwrites the point we just prepended
710  tracks[n].SetStart(vtx->GetXYZ());
711  }
712  else{
713  // Overwrites the point we just prepended
714  tracks[n].SetStart(xyz[tracks[n].View()], xyz.Z());
715  }
716  }
717  } // end for n
718  }
float Dot(const Proxy &v) const
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
TVector3 GetXYZ() const
Definition: Vertex.cxx:45
TVector3 Unit() const
Definition: View.py:1
int CountUnexplainedOnLine(const rb::HitMap &hmap, TVector3 start, TVector3 end, geo::View_t view) const
hmap should contain cells that are hit and so OK to pass through
void Add(const art::Ptr< rb::CellHit > &chit)
Definition: HitMap.cxx:37
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::ViewMerger::MergeAdjacentTracks ( std::list< rb::Track > &  trks) const
protected

Any tracks that are adjacent along the whole length of the shorter one will be merged

Definition at line 849 of file ViewMerger_module.cc.

References Adjacent(), allTimeWatchdog::endl, and LOG_DEBUG.

Referenced by produce().

850  {
851  for(std::list<rb::Track>::iterator iti = trks.begin(); iti != trks.end(); ++iti){
852  for(std::list<rb::Track>::iterator itj = trks.begin(); itj != trks.end(); ++itj){
853  if(iti == itj) continue;
854 
855  if(iti->ExtentPlane() < itj->ExtentPlane()) continue;
856 
857  if(Adjacent(*iti, *itj)){
858  LOG_DEBUG("ViewMerger") << "Adjacent: " << *iti << " and " << *itj << std::endl;
859 
860  iti->Add(itj->AllCells());
861  trks.erase(itj);
862  MergeAdjacentTracks(trks); // Start again from the top
863  return;
864  }
865  } // end for itj
866  } // end for iti
867  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void MergeAdjacentTracks(std::list< rb::Track > &trks) const
bool Adjacent(const rb::Track &a, const rb::Track &b) const
double dt::ViewMerger::MissingChargeByExtension ( const rb::HitMap sliceMap,
const rb::Track aTrk,
const rb::Track bTrk 
) const
protected

If we extend tracks to match in length, how many cells do they pass over but don't hit? Use the mean energy deposit in these tracks to estimate how much charge has gone missing, in GeV.

Definition at line 514 of file ViewMerger_module.cc.

References ana::assert(), CountUnexplainedOnLine(), geo::kX, geo::kY, rb::Cluster::NCell(), rb::Prong::Start(), rb::Track::Stop(), TotalChargePerView(), rb::Cluster::View(), and ZipTracks().

Referenced by ScoreForExtension(), and ScoreForJoinPlusExtension().

517  {
518  // UnexplainedByJoinPlusExtension doesn't know any better. Easier to fix it
519  // here.
520  if(xTrk.View() == geo::kY){
521  assert(yTrk.View() == geo::kX);
522  return MissingChargeByExtension(slice, yTrk, xTrk);
523  }
524 
525  const rb::Track& zip = ZipTracks(xTrk, yTrk);
526 
527  const int unexpX =
528  CountUnexplainedOnLine(slice, zip.Start(), xTrk.Start(), geo::kX)+
529  CountUnexplainedOnLine(slice, zip.Stop(), xTrk.Stop(), geo::kX);
530 
531  const int unexpY =
532  CountUnexplainedOnLine(slice, zip.Start(), yTrk.Start(), geo::kY)+
533  CountUnexplainedOnLine(slice, zip.Stop(), yTrk.Stop(), geo::kY);
534 
535  // End up passing through too many live cells that were unhit
536  if(unexpX+unexpY > 4) return 1e12;
537 
538  double xdens, ydens;
539  TotalChargePerView(zip, xdens, ydens);
540  xdens /= xTrk.NCell();
541  ydens /= yTrk.NCell();
542 
543  return unexpX*xdens+unexpY*ydens;
544  }
rb::Track ZipTracks(const rb::Track &aTrk, const rb::Track &bTrk) const
Make a 3D track based on two 2D tracks.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
virtual TVector3 Start() const
Definition: Prong.h:73
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void TotalChargePerView(const rb::Track &zip, double &totx, double &toty) const
assert(nhit_max >=nhit_nbins)
int CountUnexplainedOnLine(const rb::HitMap &hmap, TVector3 start, TVector3 end, geo::View_t view) const
hmap should contain cells that are hit and so OK to pass through
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
double MissingChargeByExtension(const rb::HitMap &sliceMap, const rb::Track &aTrk, const rb::Track &bTrk) const
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
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void dt::ViewMerger::produce ( art::Event evt)
virtual

Implements art::EDProducer.

Definition at line 191 of file ViewMerger_module.cc.

References util::CreateAssn(), art::Event::event(), fDebugFile, FindVertex(), PandAna.Demos.pi0_spectra::fmt, fObeyPreselection, fSliceInput, fTrackInput, art::DataViewImpl::getByLabel(), MECModelEnuComparisons::i, rb::IsFiltered(), it, geo::kX, MatchTracks(), MatchVerticals(), MaybeExtendToVertex(), MergeAdjacentTracks(), getGoodRuns4SAM::n, art::PtrVector< T >::push_back(), art::Event::put(), and art::PtrVector< T >::size().

192  {
193  std::unique_ptr<std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
194  std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
195  // Associate tracks with the slice they came from
196  std::unique_ptr<art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
197 
199  evt.getByLabel(fSliceInput, hSlices);
200 
202  for(unsigned int i = 0; i < hSlices->size(); ++i)
203  slices.push_back(art::Ptr<rb::Cluster>(hSlices, i));
204 
205  const unsigned int sliceMax = slices.size();
207  for(unsigned int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
208  if(slices[sliceIdx]->IsNoise()) continue;
209 
210  if(fObeyPreselection && rb::IsFiltered(evt, slices[sliceIdx])) continue;
211 
212  // Get all the tracks associated with this slice
213  std::vector< art::Ptr<rb::Track> > pTracks = fmt.at(sliceIdx);
214 
215  if(pTracks.empty()) continue;
216 
217  const rb::HitMap sliceMap(slices[sliceIdx]->AllCells());
218 
219  // In a nicer form
220  std::list<rb::Track> xtracks, ytracks;
221  for(unsigned int trkIdx = 0; trkIdx < pTracks.size(); ++trkIdx){
222  // No sensible way to merge these, and they make a mess
223  if(pTracks[trkIdx]->NCell() == 1) continue;
224  if(pTracks[trkIdx]->View() == geo::kX)
225  xtracks.push_back(*pTracks[trkIdx]);
226  else
227  ytracks.push_back(*pTracks[trkIdx]);
228  } // end for trkIdx
229 
230  MergeAdjacentTracks(xtracks);
231  MergeAdjacentTracks(ytracks);
232 
233  // Actually do all the work
234  std::vector<rb::Track> matchedTracks = MatchTracks(sliceMap,
235  xtracks, ytracks);
236  std::vector<rb::Track> verts = MatchVerticals(xtracks, ytracks);
237  matchedTracks.insert(matchedTracks.end(), verts.begin(), verts.end());
238 
239  // Write summary to debug file
240  fprintf(fDebugFile, "%d: ", evt.event());
241 
242  if(!matchedTracks.empty()){
243  fprintf(fDebugFile, "3d: ");
244  for(unsigned int n = 0; n < matchedTracks.size(); ++n){
245  fprintf(fDebugFile, "%u ", matchedTracks[n].NCell());
246  }
247  }
248 
249  if(!xtracks.empty() || !ytracks.empty()){
250  fprintf(fDebugFile, "2d: ");
251  typedef std::list<rb::Track>::iterator it_t;
252  for(it_t it = xtracks.begin(); it != xtracks.end(); ++it)
253  fprintf(fDebugFile, "%u(x) ", it->NCell());
254  for(it_t it = ytracks.begin(); it != ytracks.end(); ++it)
255  fprintf(fDebugFile, "%u(y) ", it->NCell());
256  }
257 
258  fprintf(fDebugFile, "\n");
259 
260  const rb::Vertex* vtx = FindVertex(matchedTracks);
261  if(vtx){
262  MaybeExtendToVertex(vtx, matchedTracks);
263  vtxcol->push_back(*vtx);
264  }
265  delete vtx;
266 
267  // This can look pretty silly and I'm not sure it adds anything.
268  /*
269  // Invent a third dimension for the unmatched 2D tracks and include them
270  // in the list.
271  std::vector<rb::Track> fakes;
272 
273  for(std::list<rb::Track>::iterator it = xtracks.begin(); it != xtracks.end(); ++it){
274  fakes.push_back(FakeThirdDimension(*it, matchedTracks));
275  }
276  for(std::list<rb::Track>::iterator it = ytracks.begin(); it != ytracks.end(); ++it){
277  fakes.push_back(FakeThirdDimension(*it, matchedTracks));
278  }
279 
280  matchedTracks.insert(matchedTracks.end(), fakes.begin(), fakes.end());
281  */
282 
283  for(unsigned int n = 0; n < matchedTracks.size(); ++n){
284  trackcol->push_back(matchedTracks[n]);
285  util::CreateAssn(*this, evt, *trackcol, slices[sliceIdx], *assns);
286  }
287  } // end for sliceIdx
288 
289  evt.put(std::move(trackcol));
290  evt.put(std::move(vtxcol));
291  evt.put(std::move(assns));
292  }
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.
void MergeAdjacentTracks(std::list< rb::Track > &trks) const
set< int >::iterator it
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::string fTrackInput
Module name to look for 2D tracks under.
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
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::string fSliceInput
Module name to look for slices under.
Definition: View.py:1
size_type size() const
Definition: PtrVector.h:308
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool fObeyPreselection
Check rb::IsFiltered?
std::vector< rb::Track > MatchVerticals(std::list< rb::Track > &xtracks, std::list< rb::Track > &ytracks) const
Find vertical tracks seperated by one plane. Return results of merge.
void MaybeExtendToVertex(const rb::Vertex *vtx, std::vector< rb::Track > &tracks) const
std::vector< rb::Track > MatchTracks(const rb::HitMap &sliceMap, std::list< rb::Track > &xtracks, std::list< rb::Track > &ytracks) const
Do most of the work. Called from produce.
rb::Vertex * FindVertex(std::vector< rb::Track > trks) const
Toy vertex finder.
void dt::ViewMerger::reconfigure ( const fhicl::ParameterSet pset)
virtual

Definition at line 168 of file ViewMerger_module.cc.

References fclose(), fDebugFile, fObeyPreselection, fSliceInput, fTrackInput, fhicl::ParameterSet::get(), and string.

Referenced by ViewMerger().

169  {
170  fSliceInput = pset.get<std::string>("SliceInput");
171  fTrackInput = pset.get<std::string>("TrackInput");
172 
174  fDebugFile = fopen(pset.get<std::string>("DebugFile").c_str(), "w");
175 
176  fObeyPreselection = pset.get<bool>("ObeyPreselection");
177  }
std::string fTrackInput
Module name to look for 2D tracks under.
fclose(fg1)
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fSliceInput
Module name to look for slices under.
bool fObeyPreselection
Check rb::IsFiltered?
enum BeamMode string
double dt::ViewMerger::ScoreForExtension ( const rb::HitMap sliceMap,
const rb::Track xTrk,
const rb::Track yTrk 
) const
protected

Definition at line 547 of file ViewMerger_module.cc.

References LOG_DEBUG, MatchScore(), rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), and MissingChargeByExtension().

Referenced by MatchTracks().

550  {
551  // No overlap: MatchScore requires that anyway
552  if(xTrk.MinPlane() > yTrk.MaxPlane()) return 1e12;
553  if(yTrk.MinPlane() > xTrk.MaxPlane()) return 1e12;
554 
555  const double miss = MissingChargeByExtension(sliceMap, xTrk, yTrk);
556 
557  // End up passing through too many live cells that were unhit
558  if(miss > 5e11) return 1e12;
559 
560  // Disable charge difference component for now.
561  // const double diff = ChargeDifferenceBetweenViews(ZipTracks(xTrk, yTrk));
562 
563  const double score = MatchScore(miss/*+diff*/, xTrk, yTrk);
564 
565  LOG_DEBUG("ViewMerger") << "Missing: " << miss
566  << ", score = " << score;
567  // << ", charge difference " << diff << std::endl;
568 
569  return score;
570  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
double MatchScore(double unexp, const rb::Track &aTrk, const rb::Track &bTrk) const
Find number of overlapping planes, return unexp/overlap.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
double MissingChargeByExtension(const rb::HitMap &sliceMap, const rb::Track &aTrk, const rb::Track &bTrk) const
double dt::ViewMerger::ScoreForJoinPlusExtension ( const rb::HitMap sliceMap,
const rb::Track single,
const rb::Track parta,
const rb::Track partb 
) const
protected

Definition at line 573 of file ViewMerger_module.cc.

References ChargeDifferenceBetweenViews(), CountUnexplainedOnLine(), release_diff::diff, allTimeWatchdog::endl, JoinTracks(), LOG_DEBUG, MatchScore(), rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), MissingChargeByExtension(), rb::Cluster::NCell(), rb::Prong::Start(), rb::Track::Stop(), TotalChargePerView(), rb::Cluster::View(), and ZipTracks().

Referenced by MatchTracks().

577  {
578  // Must be in this order, must not overlap
579  if(partb.MinPlane() < parta.MaxPlane()) return 1e12;
580 
581  // No overlap: MatchScore requires that anyway
582  if(parta.MinPlane() > single.MaxPlane()) return 1e12;
583  if(single.MinPlane() > partb.MaxPlane()) return 1e12;
584 
585  const rb::Track& join = JoinTracks(parta, partb);
586 
587  const rb::Track& zip = ZipTracks(single, join);
588  double totx, toty;
589  TotalChargePerView(zip, totx, toty);
590  const double avgDens = (totx+toty)/zip.NCell();
591 
592  // Don't allow join to use existing hits, only extension
593  const int unexpJoin = CountUnexplainedOnLine(rb::HitMap(),
594  parta.Stop(), partb.Start(),
595  join.View());
596  // Join shouldn't go through any unhit cells or others' cells. This
597  // provides some protection against crazy merges.
598  if(unexpJoin > 0) return 1e12;
599 
600  // Charge missing at the ends
601  const double miss = MissingChargeByExtension(sliceMap, single, join);
602 
603  // End up passing through too many live cells that were unhit
604  if(miss > 4*avgDens) return 1e12;
605 
606  const double diff = ChargeDifferenceBetweenViews(ZipTracks(single, join));
607 
608  // Small penalty for having a join
609  const double score = MatchScore(miss+avgDens+diff, single, join);
610 
611  LOG_DEBUG("ViewMerger") << "Join: missing charge " << miss
612  << ", score = " << score
613  << ", charge difference: " << diff << std::endl;
614 
615  return score;
616  }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
rb::Track ZipTracks(const rb::Track &aTrk, const rb::Track &bTrk) const
Make a 3D track based on two 2D tracks.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
virtual TVector3 Start() const
Definition: Prong.h:73
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
void TotalChargePerView(const rb::Track &zip, double &totx, double &toty) const
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
double MatchScore(double unexp, const rb::Track &aTrk, const rb::Track &bTrk) const
Find number of overlapping planes, return unexp/overlap.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
rb::Track JoinTracks(const rb::Track &parta, const rb::Track &partb) const
Append two tracks in the same view (crudely).
int CountUnexplainedOnLine(const rb::HitMap &hmap, TVector3 start, TVector3 end, geo::View_t view) const
hmap should contain cells that are hit and so OK to pass through
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
double ChargeDifferenceBetweenViews(const rb::Track &zip) const
double MissingChargeByExtension(const rb::HitMap &sliceMap, const rb::Track &aTrk, const rb::Track &bTrk) const
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

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

void dt::ViewMerger::TotalChargePerView ( const rb::Track zip,
double &  totx,
double &  toty 
) const
protected

Corrected for direction cosines, so if you divide this by extent in Z you'd get dE/dx

Definition at line 486 of file ViewMerger_module.cc.

References rb::Cluster::Cell(), dz, e, rb::RecoHit::GeV(), rb::Track::InterpolateXY(), rb::RecoHit::IsCalibrated(), geo::kX, getGoodRuns4SAM::n, rb::Cluster::NCell(), util::pythag(), rb::Cluster::RecoHit(), rb::CellHit::View(), x1, y1, and rb::RecoHit::Z().

Referenced by ChargeDifferenceBetweenViews(), MissingChargeByExtension(), and ScoreForJoinPlusExtension().

488  {
489  // TODO: accounting for uncalibrated hits and dead channels
490 
491  totx = 0;
492  toty = 0;
493 
494  for(unsigned int n = 0; n < zip.NCell(); ++n){
495  const art::Ptr<rb::CellHit> chit = zip.Cell(n);
496  const rb::RecoHit rhit = zip.RecoHit(chit);
497  if(rhit.IsCalibrated()){
498  // Estimate track direction here
499  const double eps = 1e-3;
500  double x0, y0, x1, y1;
501  zip.InterpolateXY(rhit.Z(), x0, y0);
502  zip.InterpolateXY(rhit.Z()+eps, x1, y1);
503  const double dz = eps/util::pythag(eps, x1-x0, y1-y0);
504 
505  if(chit->View() == geo::kX)
506  totx += rhit.GeV()*dz;
507  else
508  toty += rhit.GeV()*dz;
509  }
510  }
511  }
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
double dz[NP][NC]
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float GeV() const
Definition: RecoHit.cxx:69
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
Float_t e
Definition: plot.C:35
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited
rb::Track dt::ViewMerger::ZipTracks ( const rb::Track aTrk,
const rb::Track bTrk 
) const
protected

Make a 3D track based on two 2D tracks.

Definition at line 739 of file ViewMerger_module.cc.

References rb::Track::ZipWith().

Referenced by MatchTracks(), MissingChargeByExtension(), and ScoreForJoinPlusExtension().

741  {
742  return xTrk.ZipWith(yTrk);
743 
744  // TODO: restore protection against going up to a plane too far. Can just
745  // truncate that returned track.
746 
747  /*
748  if(xTrk.View() == geo::kY){
749  assert(yTrk.View() == geo::kX);
750  return ZipTracks(yTrk, xTrk);
751  }
752 
753  art::ServiceHandle<geo::Geometry> geom;
754 
755  std::vector<TVector3> trajPts;
756 
757  double minZ = std::min(xTrk.Start().Z(), yTrk.Start().Z());
758  double maxZ = std::max(xTrk.Stop().Z(), yTrk.Stop().Z());
759 
760  const int minPlane = std::min(xTrk.MinPlane(), yTrk.MinPlane());
761  const int maxPlane = std::max(xTrk.MaxPlane(), yTrk.MaxPlane());
762 
763  // 2D tracks often stick halfway into the next plane. Don't allow the final
764  // track to go further than we have any evidence of hits. In that case put
765  // it to within 1mm of the end of the last plane with hits.
766  double xyz[3], dxyz[3];
767  geom->CellInfo(minPlane, 0, 0, xyz, dxyz);
768  minZ = std::max(minZ, xyz[2]-dxyz[2]+.1);
769 
770  geom->CellInfo(maxPlane, 0, 0, xyz, dxyz);
771  maxZ = std::min(maxZ, xyz[2]+dxyz[2]-.1);
772 
773  const double kZStep = 5; // cm
774  const int stepMax = int((maxZ-minZ)/kZStep);
775  for(int step = 0; step <= stepMax; ++step){
776  const double z = minZ+(double(step)/stepMax)*(maxZ-minZ);
777 
778  double xyz[3];
779  xyz[2] = z;
780 
781  double junk;
782  xTrk.InterpolateXY(z, xyz[0], junk);
783  yTrk.InterpolateXY(z, junk, xyz[1]);
784 
785  trajPts.push_back(xyz);
786  } // end for step
787 
788  rb::Track ret(art::PtrVector<rb::CellHit>(),
789  trajPts[0], trajPts[1]-trajPts[0]);
790  ret.Add(xTrk.AllCells());
791  ret.Add(yTrk.AllCells());
792 
793  for(unsigned int n = 1; n < trajPts.size(); ++n)
794  ret.AppendTrajectoryPoint(trajPts[n]);
795 
796  return ret;
797  */
798  }
rb::Track dt::ViewMerger::ZipVerticalTracks ( const rb::Track xtrk,
const rb::Track ytrk 
) const
protected

Tilt tracks arbitrarily and merge them into a 3D track.

Definition at line 429 of file ViewMerger_module.cc.

References rb::Cluster::AllCells(), geo::GeometryBase::CellInfo(), geom(), geo::kX, geo::kY, rb::Cluster::MaxCell(), rb::Cluster::MinCell(), rb::Cluster::MinPlane(), runNovaSAM::ret, febshutoff_auto::start, x1, and y1.

Referenced by MatchVerticals().

431  {
433 
434  double xyz[3], dxyz[3];
435  geom->CellInfo(xtrk.MinPlane(), xtrk.MinCell(geo::kX), 0, xyz, dxyz);
436  const double x0 = xyz[0];
437  const double xz0 = xyz[2]-dxyz[2];
438  geom->CellInfo(xtrk.MinPlane(), xtrk.MaxCell(geo::kX), 0, xyz, dxyz);
439  const double x1 = xyz[0];
440  const double xz1 = xyz[2]+dxyz[2];
441 
442  geom->CellInfo(ytrk.MinPlane(), ytrk.MinCell(geo::kY), 0, xyz, dxyz);
443  const double y0 = xyz[1];
444  const double yz0 = xyz[2]-dxyz[2];
445  geom->CellInfo(ytrk.MinPlane(), ytrk.MaxCell(geo::kY), 0, xyz, dxyz);
446  const double y1 = xyz[1];
447  const double yz1 = xyz[2]+dxyz[2];
448 
449  // If we had any idea where the vertex was we could potentially swap x0, x1
450  // or y0, y1.
451 
452  const double dxdz = (x1-x0)/(xz1-xz0);
453  const double dydz = (y1-y0)/(yz1-yz0);
454 
455  TVector3 start, stop;
456 
457  if(xz0+xz1 < yz0+yz1){
458  // x is upstream of y
459  start.SetX(x0);
460  start.SetY(y0+dydz*(xz0-yz0));
461  start.SetZ(xz0);
462 
463  stop.SetX(x0+dxdz*(yz1-xz0));
464  stop.SetY(y1);
465  stop.SetZ(yz1);
466  }
467  else{
468  // y is upstream of x
469  start.SetX(x0+dxdz*(yz0-xz0));
470  start.SetY(y0);
471  start.SetZ(yz0);
472 
473  stop.SetX(x1);
474  stop.SetY(y1+dydz*(xz1-yz1));
475  stop.SetZ(xz1);
476  }
477 
478  rb::Track ret(xtrk.AllCells(), start, stop-start);
479  ret.Add(ytrk.AllCells());
480  ret.AppendTrajectoryPoint(stop);
481 
482  return ret;
483  }
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
void geom(int which=0)
Definition: geom.C:163

Member Data Documentation

FILE* dt::ViewMerger::fDebugFile
protected

Definition at line 147 of file ViewMerger_module.cc.

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

bool dt::ViewMerger::fObeyPreselection
protected

Check rb::IsFiltered?

Definition at line 145 of file ViewMerger_module.cc.

Referenced by produce(), and reconfigure().

std::string dt::ViewMerger::fSliceInput
protected

Module name to look for slices under.

Definition at line 142 of file ViewMerger_module.cc.

Referenced by produce(), and reconfigure().

std::string dt::ViewMerger::fTrackInput
protected

Module name to look for 2D tracks under.

Definition at line 143 of file ViewMerger_module.cc.

Referenced by produce(), and reconfigure().


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