Cluster.cxx
Go to the documentation of this file.
1 #include "DDTMonopole/Cluster.h"
2 
4 
5 #include <iostream>
6 
7 
8 
10  max_cell_gap_cut_ (p.get<int>("max_cell_gap_cut")),
11  max_plane_gap_cut_(p.get<int>("max_plane_gap_cut")),
12  track_min_beta_ (p.get<double>("track_min_beta")),
13  track_max_beta_ (p.get<double>("track_max_beta"))
14 {
15 }
16 
17 
18 
20 {
21  return monopole_hits_;
22 }
23 
24 
25 
27 {
29 
30  bool result = false;
31 
32  std::map<std::string, mono::HitList> hitmap;
33  split_by_view(hits, hitmap);
34 
35  std::cout << "MF: (Total Hits, X Hits, Y Hits) = ("
36  << hits.size()
37  << ", " << hitmap["x"].size()
38  << ", " << hitmap["y"].size()
39  << ")" << std::endl;
40 
41  if (hitmap["x"].size() < novaddt::smt::Constants::TRACK_MIN_HITS ||
43  return false;
44 
45  hitmap["x_surface"];
46  hitmap["x_contained"];
47 
48  find_view_matched_hits(hitmap);
49 
50  if(hitmap["x_contained"].size() < novaddt::smt::Constants::TRACK_MIN_HITS)
51  return false;
52 
53  for (auto entry_hit = hitmap["x_surface"].cbegin();
54  entry_hit != hitmap["x_surface"].cend() && !result;
55  ++entry_hit)
56  {
57  // we add 1, so that we avoid combinations between the same hit
58  mono::HitList::const_reverse_iterator until(entry_hit + 1);
59 
60  auto contained_first = find(hitmap["x_contained"], (*entry_hit)->TDC());
61 
62  std::find_if(hitmap["x_surface"].crbegin(), until,
63  [&](art::Ptr<rb::CellHit> const& exit_hit)
64  {
65  auto contained_last =
66  find(contained_first, hitmap["x_contained"].cend(), exit_hit->TDC());
67 
68  if (std::distance(contained_first, contained_last) <
70  return false;
71 
72  Track track(*entry_hit, exit_hit);
73 
74  if (!good_track(track))
75  return false;
76 
77  std::set<int> plane_set = { (*entry_hit)->Plane() };
78  std::set<int> cell_set = { (*entry_hit)->Cell() };
79 
80  int32_t t_max_plane_gap = track.start_rb()->TDC();
81  if (track.slope_plane_is_valid())
82  t_max_plane_gap += abs(track.slope_time_plane()) * max_plane_gap_cut_;
83 
84  int32_t t_max_cell_gap = track.start_rb()->TDC();
85  if (track.slope_cell_is_valid())
86  t_max_cell_gap += abs(track.slope_time_cell()) * max_cell_gap_cut_;
87  int32_t t_max_gap = std::min(t_max_plane_gap, t_max_cell_gap);
88 
89  unsigned n_contained_hit = 0;
90  bool contains_gap_exceeding_cut = std::any_of
91  (contained_first, contained_last,
92  [&](art::Ptr<rb::CellHit> const& hit)
93  {
94  if (hit->TDC() > t_max_gap)
95  {
96  int64_t time_difference =
97  static_cast<int64_t>(hit->TDC()) -
98  static_cast<int64_t>(track.start_rb()->TDC());
99 
100  if (hit->TDC() > t_max_cell_gap)
101  {
102  double expected_cell =
103  static_cast<double>(track.start_rb()->Cell());
104  if (track.slope_time_is_valid())
105  expected_cell += track.slope_cell_time() * time_difference;
106 
107  if (gap_exceeds_cut(cell_set, expected_cell, max_cell_gap_cut_))
108  return true;
109  }
110 
111  if (hit->TDC() > t_max_plane_gap)
112  {
113  double expected_plane =
114  static_cast<double>(track.start_rb()->Plane());
115  if (track.slope_time_is_valid())
116  expected_plane += track.slope_plane_time() * time_difference;
117 
118  if (gap_exceeds_cut
119  (plane_set, expected_plane, max_plane_gap_cut_))
120  return true;
121  }
122  }
123 
124  if (hit_is_in_time_with_road(hit, track))
125  {
126  cell_set.insert(hit->Cell());
127  plane_set.insert(hit->Plane());
128  }
129 
130  return false;
131  });
132 
133  if (contains_gap_exceeding_cut)
134  return false;
135 
136  plane_set.insert(exit_hit->Plane());
138  return false;
139 
140  if (find_max_gap(plane_set) > max_plane_gap_cut_)
141  return false;
142 
143  cell_set.insert(exit_hit->Cell());
144  if (find_max_gap(cell_set) > max_cell_gap_cut_)
145  return false;
146 
147  result = true;
148  monopole_hits_.Add(*entry_hit);
149  monopole_hits_.Add(exit_hit);
150  std::for_each
151  (contained_first, contained_last,
152  [&](art::Ptr<rb::CellHit> const& hit)
153  {
154  if (hit_is_in_time_with_road(hit, track))
155  monopole_hits_.Add(hit);
156  });
157 
158  return true;
159  });
160  }
161 
162  return result;
163 }
164 
165 
166 
168 (mono::HitList hits, std::map<std::string, mono::HitList> & hitmap) const
169 {
170  auto bound = std::stable_partition
171  (hits.begin(), hits.end(),
172  [](art::Ptr<rb::CellHit> const& h)
173  { return h->View() == geo::View_t::kX; });
174 
175  hitmap["x"] = mono::HitList(std::distance(hits.begin(), bound));
176  std::move(hits.begin(), bound, hitmap["x"].begin());
177 
178  hitmap["y"] = mono::HitList(std::distance(bound, hits.end()));
179  std::move(bound, hits.end(), hitmap["y"].begin());
180 }
181 
182 
183 
185 (std::map<std::string, mono::HitList> & hitmap) const
186 {
187  for (auto const& x_hit : hitmap["x"])
188  {
189  auto y_first =
190  find(hitmap["y"], x_hit->TDC() -
192  auto y_last =
193  find(hitmap["y"], x_hit->TDC() +
195 
196  auto y_match = std::find_if
197  (y_first, y_last, [&](art::Ptr<rb::CellHit> const& y_hit)
198  { return hits_are_view_matched(x_hit, y_hit); });
199 
200  if (y_match != y_last)
201  {
202  // we found a matched hit
203  hitmap["x_contained"].push_back(x_hit);
204 
205  // is this hit on the surface?
206  if (hit_is_on_surface(x_hit))
207  {
208  hitmap["x_surface"].push_back(x_hit);
209  } else {
210  // now we need to look for additional matched y-hits and check whether
211  // they are on the surface
212  auto surface_y_match = std::find_if
213  (y_match, y_last, [&](art::Ptr<rb::CellHit> const& y_hit)
214  {
215  if (!hits_are_view_matched(x_hit, y_hit))
216  return false;
217 
218  if (!hit_is_on_surface(y_hit))
219  return false;
220 
221  return true;
222  });
223 
224  if (surface_y_match != y_last)
225  hitmap["x_surface"].push_back(x_hit);
226  }
227  }
228  }
229 }
230 
231 
232 
234 (art::Ptr<rb::CellHit> const& x_hit, art::Ptr<rb::CellHit> const& y_hit) const
235 {
236  int plane_difference =
237  static_cast<int>(x_hit->Plane()) - static_cast<int>(y_hit->Plane());
238 
239  if (std::abs(plane_difference) <
241  return true;
242 
243  return false;
244 }
245 
246 
247 
250 {
252  return true;
253 
255  return true;
256 
258  return true;
259 
261  return true;
262 
263  return false;
264 }
265 
266 
267 
269 (art::Ptr<rb::CellHit> const& hit, Track const& track) const
270 {
271  if (hit->Plane() < track.min_plane())
272  return false;
273 
274  if (hit->Plane() > track.max_plane())
275  return false;
276 
277  if (hit->Cell() < track.min_cell())
278  return false;
279 
280  if (hit->Cell() > track.max_cell())
281  return false;
282 
283  return true;
284 }
285 
286 
287 
288 bool mono::Cluster::good_track(Track const& track) const
289 {
291  return false;
292 
293  if (track.beta() < track_min_beta_)
294  return false;
295 
296  if (track.beta() > track_max_beta_)
297  return false;
298 
299  if (abs(track.plane_difference()) <= max_plane_gap_cut_)
300  return false;
301 
302  return true;
303 }
304 
305 
306 
308 (art::Ptr<rb::CellHit> const& hit, Track const& track) const
309 {
310  if (!hit_is_between_track_end_points(hit, track))
311  return false;
312 
313  // If the plane of the hit and track end points are all equal,
314  // the hit is on the road.
315  if (!track.slope_plane_is_valid())
316  return true;
317 
318  double road_center_cell = track.slope_cell_plane() *
319  (static_cast<double>(hit->Plane()) -
320  static_cast<double>(track.start_rb()->Plane())) +
321  static_cast<double>(track.start_rb()->Cell());
322 
323  if (hit->Cell() > road_center_cell + novaddt::smt::Constants::ROAD_HALF_WIDTH)
324  return false;
325 
326  if (hit->Cell() < road_center_cell - novaddt::smt::Constants::ROAD_HALF_WIDTH)
327  return false;
328 
329  return true;
330 }
331 
332 
333 
335 (art::Ptr<rb::CellHit> const& hit, Track const& track) const
336 {
337  if (!hit_is_on_road(hit, track))
338  return false;
339 
340  // If the plane of the hit and track end points are all equal,
341  // we cannot calculate the time_slope and make no statement on the time.
342  if (!track.slope_plane_is_valid())
343  return false;
344 
345 
346  int32_t relative_hit_time = hit->TDC() - track.start_rb()->TDC();
347 
348  double expected_time = track.slope_time_plane() *
349  (static_cast<double>(hit->Plane()) -
350  static_cast<double>(track.start_rb()->Plane()));
351 
352  if (relative_hit_time >
354  return false;
355 
356  if (relative_hit_time <
358  return false;
359 
360  return true;
361 }
362 
363 
364 
366 (std::set<int> collection, int const& object, int const& max_gap_cut) const
367 {
368  // Note that we copy the collection in here on purpose, so that we do not
369  // accidentally alter it.
370  collection.insert(object);
371 
372  // We need to encode the max_plane cut differently here, so as to keep
373  // this method as general as possible.
374  if (find_max_gap(collection) > max_gap_cut)
375  return true;
376 
377  return false;
378 }
379 
380 
381 
383 (std::set<int> const& collection) const
384 {
385  int result = -9999;
386 
387  std::adjacent_find
388  (collection.begin(), collection.end(), [&](int object, int next_object)
389  {
390  int gap = next_object - object;
391  result = std::max(gap, result);
392 
393  // This return is meaningless, but a boolean return is required by
394  // adjacent find. We choose false, so the loop goes until the end.
395  return false;
396  });
397 
398  return result;
399 }
400 
401 
402 
403 mono::HitList::const_iterator mono::Cluster::find
404 (mono::HitList::const_iterator begin,
405  mono::HitList::const_iterator end, int32_t tdc) const
406 {
407  return std::lower_bound
408  (begin, end, tdc,
409  [](art::Ptr<rb::CellHit> const& h, int32_t const& tdc_val)
410  { return h->TDC() < tdc_val; });
411 }
412 
413 
414 
415 mono::HitList::const_iterator mono::Cluster::find
416 (mono::HitList const& hits, int32_t tdc) const
417 {
418  return find(hits.cbegin(), hits.cend(), tdc);
419 }
double beta() const
Definition: Track.cxx:129
T max(const caf::Proxy< T > &a, T b)
int max_plane_gap_cut_
Definition: Cluster.h:53
const unsigned ROAD_HALF_WIDTH
Definition: Constants.h:37
int plane_difference() const
Definition: Track.cxx:232
bool slope_plane_is_valid() const
Definition: Track.cxx:136
unsigned max_plane() const
Definition: Track.cxx:245
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
bool gap_exceeds_cut(std::set< int > collection, int const &object, int const &max_gap_cut) const
Definition: Cluster.cxx:366
unsigned short Plane() const
Definition: CellHit.h:39
const char * p
Definition: xmltok.h:285
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
const unsigned SURFACE_MAX_CELL
Definition: Constants.h:24
void abs(TH1 *hist)
std::vector< art::Ptr< rb::CellHit > > HitList
Definition: Cluster.h:15
unsigned distance(const T &t1, const T &t2)
void find_view_matched_hits(std::map< std::string, mono::HitList > &hitmap) const
Definition: Cluster.cxx:185
bool slope_time_is_valid() const
Definition: Track.cxx:143
float abs(float number)
Definition: d0nt_math.hpp:39
const unsigned SURFACE_MAX_PLANE
Definition: Constants.h:28
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
const double ROAD_TIME_VARIATION
Definition: Constants.h:38
int max_cell_gap_cut_
Definition: Cluster.h:53
unsigned short Cell() const
Definition: CellHit.h:40
double slope_plane_time() const
Definition: Track.cxx:202
unsigned max_cell() const
Definition: Track.cxx:266
void hits()
Definition: readHits.C:15
mono::HitList::const_iterator find(mono::HitList const &hits, int32_t tdc) const
Definition: Cluster.cxx:416
bool hit_is_on_surface(art::Ptr< rb::CellHit > const &hit) const
Definition: Cluster.cxx:249
void split_by_view(mono::HitList hits, std::map< std::string, mono::HitList > &hitmap) const
Definition: Cluster.cxx:168
const unsigned TRACK_MIN_UNIQUE_PLANES
Definition: Constants.h:34
bool hit_is_in_time_with_road(art::Ptr< rb::CellHit > const &hit, Track const &track) const
Definition: Cluster.cxx:335
const unsigned VIEW_MATCH_PLANE_CUT
Definition: Constants.h:31
bool hit_is_between_track_end_points(art::Ptr< rb::CellHit > const &hit, Track const &track) const
Definition: Cluster.cxx:269
rb::Cluster get() const
Definition: Cluster.cxx:19
art::Ptr< rb::CellHit > start_rb() const
Definition: Track.cxx:344
const unsigned SURFACE_MIN_PLANE
Definition: Constants.h:27
bool make(std::vector< art::Ptr< rb::CellHit > > const &hits)
Definition: Cluster.cxx:26
const double TRACK_MIN_DISTANCE
Definition: Constants.h:35
bool good_track(Track const &track) const
Definition: Cluster.cxx:288
OStream cout
Definition: OStream.cxx:6
bool slope_cell_is_valid() const
Definition: Track.cxx:150
double slope_cell_plane() const
Definition: Track.cxx:157
double slope_time_plane() const
Definition: Track.cxx:172
Definition: structs.h:12
rb::Cluster monopole_hits_
Definition: Cluster.h:55
double slope_time_cell() const
Definition: Track.cxx:217
Definition: geo.h:1
Cluster(fhicl::ParameterSet const &p)
Definition: Cluster.cxx:9
unsigned min_cell() const
Definition: Track.cxx:259
bool hit_is_on_road(art::Ptr< rb::CellHit > const &hit, Track const &track) const
Definition: Cluster.cxx:308
int find_max_gap(std::set< int > const &collection) const
Definition: Cluster.cxx:383
T min(const caf::Proxy< T > &a, T b)
const unsigned SURFACE_MIN_CELL
Definition: Constants.h:23
const unsigned VIEW_MATCH_TDC_CUT
Definition: Constants.h:30
double track_min_beta_
Definition: Cluster.h:54
double distance() const
Definition: Track.cxx:108
double slope_cell_time() const
Definition: Track.cxx:187
Float_t track
Definition: plot.C:35
bool hits_are_view_matched(art::Ptr< rb::CellHit > const &x_hit, art::Ptr< rb::CellHit > const &y_hit) const
Definition: Cluster.cxx:234
virtual void Clear()
Forget about all owned cell hits.
Definition: Cluster.cxx:279
const unsigned TRACK_MIN_HITS
Definition: Constants.h:33
unsigned min_plane() const
Definition: Track.cxx:238
double track_max_beta_
Definition: Cluster.h:54