MonopoleTrack_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MonopoleTrack
3 // Module Type: filter
4 // File: MonopoleTrack_module.cc
5 //
6 // Generated at Wed Jul 8 13:21:15 2015 by Martin Frank using artmod
7 // from cetpkgsupport v1_08_06.
8 ////////////////////////////////////////////////////////////////////////
9 
11 // #include <canvas/Persistency/Common/FindManyP.h>
12 #include <art/Framework/Core/FindManyP.h>
19 // #include <canvas/Utilities/InputTag.h>
20 #include <art/Utilities/InputTag.h>
21 #include <fhiclcpp/ParameterSet.h>
23 
24 #include "Geometry/Geometry.h"
25 #include "Monopole/Track2D.h"
26 #include "Monopole/Track3D.h"
27 #include "RecoBase/CellHit.h"
28 #include "RecoBase/Cluster.h"
29 #include "RecoBase/HoughResult.h"
30 #include "RecoBase/Track.h"
31 #include "Utilities/func/MathUtil.h"
32 
33 #include <iostream>
34 #include <memory>
35 #include <vector>
36 
37 namespace mono
38 {
39  struct Track2D;
40  struct Track_Match;
41  class MonopoleTrack;
42 }
43 
45 {
46  Track_Match(std::vector<Track2D>::const_iterator x,
47  std::vector<Track2D>::const_iterator y, double s) :
48  x_track(x), y_track(y), score(s) {}
49 
50  std::vector<Track2D>::const_iterator x_track, y_track;
51  double score;
52 };
53 
55 {
56 public:
57  explicit MonopoleTrack(fhicl::ParameterSet const & p);
58 
59  MonopoleTrack(MonopoleTrack const &) = delete;
60  MonopoleTrack(MonopoleTrack &&) = delete;
61  MonopoleTrack & operator = (MonopoleTrack const &) = delete;
62  MonopoleTrack & operator = (MonopoleTrack &&) = delete;
63 
64  bool filter(art::Event & e) override;
65 
66 private:
67  bool hit_is_on_track
69  double const& slope, double const& intercept) const;
70 
75  unsigned track_min_hits_;
76 
78 };
79 
80 
81 
83  hough_label_ (p.get<std::string>("hough_label")),
84  hough_peak_cut_ (p.get<double>("hough_peak_cut")),
85  slice_label_ (p.get<std::string>("slice_label")),
86  track_matching_min_score_(p.get<double>("track_matching_min_score")),
87  track_min_hits_ (p.get<unsigned>("track_min_hits"))
88 {
89  produces<std::vector<rb::Cluster> >();
90  produces<std::vector<mono::Track3D> >();
91 }
92 
93 
94 
96 {
97  //
98  // Extract Objects from the Event
99  //
101  e.getByLabel(slice_label_, mono_slices);
102 
104  find_hough_results(mono_slices, e, hough_label_);
105 
106 
107 
108  //
109  // Track the Monopoles
110  //
111  std::map<geo::View_t, std::vector<Track2D> > mono_tracks;
112  for (size_t n_slice = 0; n_slice != mono_slices->size(); ++n_slice)
113  {
114  if (!mono_slices->at(n_slice).IsNoise())
115  {
116  std::vector<art::Ptr<rb::HoughResult> > hough_results;
117  find_hough_results.get(n_slice, hough_results);
118  for (auto const& hough_result : hough_results)
119  {
120  unsigned n_peak = 0;
121  for (auto const& peak : hough_result->fPeak)
122  {
123  if (peak.fH > hough_peak_cut_)
124  {
125  geo::View_t view = hough_result->fView;
126  double slope = -9999;
127  double intercept = -9999;
128  hough_result->SlopeIcept(n_peak++, &slope, &intercept);
129 
130  Track2D mono_track(view, slope, intercept);
131  for (auto const& hit : mono_slices->at(n_slice).AllCells())
132  if (hit_is_on_track(hit, view, slope, intercept))
133  mono_track.cluster.Add(hit);
134  mono_tracks[view].push_back(mono_track);
135  }
136  }
137  }
138  }
139  }
140 
141 
142 
143  //
144  // Merge Tracks
145  //
146  std::vector<Track_Match> matches;
147  for (auto x_track = mono_tracks[geo::kX].begin();
148  x_track != mono_tracks[geo::kX].end(); ++x_track)
149  {
150  for (auto y_track = mono_tracks[geo::kY].begin();
151  y_track != mono_tracks[geo::kY].end(); ++y_track)
152  {
153  // for each x-hit, we need to see whether there is a matching y-hit
154  // let's count one point per match
155  // this can be encapsulated into a nice separate function
156  double score = 0.0;
157  for (auto const& x_hit : x_track->cluster.AllCells())
158  {
159  for (auto const& y_hit : y_track->cluster.AllCells())
160  {
161  if (abs(int(x_hit->Plane()) - int(y_hit->Plane())) <= 1)
162  score += 1;
163  }
164  }
165 
166  if (score >= track_matching_min_score_)
167  matches.emplace_back(x_track, y_track, score);
168  }
169  }
170 
171  std::sort(matches.begin(), matches.end(),
172  [](Track_Match const& lhs, Track_Match const& rhs)
173  { return lhs.score > rhs.score; });
174 
175  std::set<std::vector<Track2D>::const_iterator> already_matched_tracks;
176  std::vector<Track_Match> good_track_matches;
177  for (auto const& match : matches)
178  {
179  if (already_matched_tracks.find(match.x_track) ==
180  already_matched_tracks.end() &&
181  already_matched_tracks.find(match.y_track) ==
182  already_matched_tracks.end())
183  {
184  good_track_matches.push_back(match);
185  already_matched_tracks.insert(match.x_track);
186  already_matched_tracks.insert(match.y_track);
187  }
188  }
189 
190 
191 
192  //
193  // Add Objects to the Event
194  //
195  std::unique_ptr<std::vector<Track3D> >
196  tracks_3D(new std::vector<Track3D>);
197  for (auto const& good_match : good_track_matches)
198  {
199  // perhaps we can pass the iterators through to make things more efficient
200  Track3D track3D(*(good_match.x_track), *(good_match.y_track));
201 
202  if (track3D.cluster().NCell() >= track_min_hits_)
203  tracks_3D->push_back(track3D);
204  }
205 
206  std::sort(tracks_3D->begin(), tracks_3D->end(),
207  [](Track3D const& lhs, Track3D const& rhs)
208  { return lhs.beta() < rhs.beta(); });
209 
210  std::unique_ptr<std::vector<rb::Cluster> >
211  clusters_3D(new std::vector<rb::Cluster>);
212  for (auto const& track : *tracks_3D)
213  clusters_3D->push_back(track.cluster());
214 
215  e.put(std::move(clusters_3D));
216  e.put(std::move(tracks_3D));
217 
218  return true;
219 }
220 
221 
222 
225  double const& slope, double const& intercept) const
226 {
227  if (hit->View() != view)
228  return false;
229 
230  double xyz[3];
231  geometry_->CellInfo(hit->Plane(), hit->Cell(), &view, xyz);
232  TVector3 hit_vector(xyz);
233 
234  double hit_position = -9999;
235  if (view == geo::View_t::kX)
236  hit_position = hit_vector.x();
237  else if (view == geo::View_t::kY)
238  hit_position = hit_vector.y();
239 
240  double track_position = slope * hit_vector.z() + intercept;
241  double distance_to_track = std::fabs(track_position - hit_position);
242 
243  if (distance_to_track > 20)
244  return false;
245 
246  return true;
247 }
248 
249 
250 
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< Track2D >::const_iterator y_track
Definition: Cluster.h:13
double beta() const
Definition: Track3D.cxx:178
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
std::vector< Track2D >::const_iterator x_track
void abs(TH1 *hist)
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
Track_Match(std::vector< Track2D >::const_iterator x, std::vector< Track2D >::const_iterator y, double s)
rb::Cluster cluster() const
Definition: Track3D.cxx:185
bool filter(art::Event &e) override
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Data resulting from a Hough transform on the cell hit positions.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
art::ServiceHandle< geo::Geometry > geometry_
rb::Cluster cluster
Definition: Track2D.h:26
Definition: geo.h:1
bool hit_is_on_track(art::Ptr< rb::CellHit > const &hit, geo::View_t view, double const &slope, double const &intercept) const
size_type get(size_type i, reference item, data_reference data) const
Definition: FindManyP.h:469
Float_t e
Definition: plot.C:35
Encapsulate the geometry of one entire detector (near, far, ndos)
MonopoleTrack(fhicl::ParameterSet const &p)