SlowMonopoleTrigger.cxx
Go to the documentation of this file.
2 #include "DDTTools/DDTHelpers.h"
4 #include "DDTMonopole/Track.h"
5 
8 
9 #include <boost/format.hpp>
10 
11 #include <algorithm>
12 #include <chrono>
13 #include <iostream>
14 
15 using namespace std::chrono;
16 
17 
18 
20 (novaddt::HitList const& h, unsigned long long const& event_time,
21  Configuration const& c) :
22  hits_(h),
23  event_time_(event_time),
24  max_plane_gap_cut_(c.max_plane_gap_cut),
25  sparsification_factor_(c.sparsification_factor)
26 {
27  max_cell_gap_cut_ = 20;
28  track_max_beta_ = 5e-3;
29  track_min_beta_ = 4e-5;
30  trigger_time_buffer_ = 256;
31 
33  algorithm();
35 
36  parameter_["algorithm_time"] = duration_cast<milliseconds>(t1-t0).count();
37 
38  std::cout << "MF: Algorithm Time: "
39  << duration_cast<milliseconds>(t1-t0).count() << " ms"
40  << std::endl;
41 }
42 
43 
44 
45 std::vector<novaddt::TriggerDecision>
47 {
48  return trigger_decisions_;
49 }
50 
51 
52 
54 {
55  return trigger_time_buffer_;
56 }
57 
58 
59 
60 
61 std::map<std::string, novaddt::HitList>
63 {
64  return essential_hit_lists_;
65 }
66 
67 
68 
70 {
71  if (parameter_.find(name) == parameter_.end())
72  return false;
73 
74  return true;
75 }
76 
77 
78 
80 {
81  if (parameter_.find(name) == parameter_.end())
82  {
83  std::cerr << "\n\tParameter " << name << " is not available.\n"
84  << std::endl;
85  assert(false);
86  }
87 
88  return parameter_[name];
89 }
90 
91 
92 
94 {
95  // Timing Cheat Sheet
96  //
97  // std::map<std::string, double> dur =
98  // { {"setup", 0}, {"find", 0}, {"track", 0}, {"sets", 0} };
99  // auto t_start = high_resolution_clock::now();
100  // dur["setup"] +=
101  // duration_cast<nanoseconds>(high_resolution_clock::now() - t_start).count();
102 
103  bool result = false;
104  novaddt::DAQHit trigger_entry_hit, trigger_exit_hit;
105 
106  std::map<std::string, novaddt::HitList> hitmap;
107  split_by_view(hits_, hitmap);
108 
109  // std::cout << "MF: (x hits, y hits, min_hits) = ("
110  // << hitmap["x"].size()
111  // << ", " << hitmap["y"].size()
112  // << ", " << min_hits_
113  // << ")" << std::endl;
114 
115  if (hitmap["x"].size() < novaddt::smt::Constants::TRACK_MIN_HITS ||
117  return false;
118 
119  hitmap["x_surface"];
120  hitmap["x_contained"];
121 
122  find_view_matched_hits(hitmap);
123 
124  // std::cout << "MF: (x surface, x contained, min_hits) = ("
125  // << hitmap["x_surface"].size()
126  // << ", " << hitmap["x_contained"].size()
127  // << ", " << min_hits_
128  // << ")" << std::endl;
129 
130  // for (auto const& hit : hitmap["x_surface"])
131  // ddthelpers::print(hit, event_time_);
132 
133  // std::cout << "\n\n\nContained Hits:" << std::endl;
134  // for (auto const& hit : hitmap["x_contained"])
135  // ddthelpers::print(hit, event_time_);
136 
137  if(hitmap["x_contained"].size() < novaddt::smt::Constants::TRACK_MIN_HITS)
138  return false;
139 
140  unsigned n_good_track = 0;
141  for (auto entry_hit = hitmap["x_surface"].cbegin();
142  entry_hit != hitmap["x_surface"].cend() && !result;
143  ++entry_hit)
144  {
145  // we add 1, so that we avoid combinations between the same hit
146  novaddt::HitList::const_reverse_iterator until(entry_hit + 1);
147 
148  auto contained_first = find(hitmap["x_contained"], entry_hit->TDC());
149 
150  std::find_if(hitmap["x_surface"].crbegin(), until,
151  [&](novaddt::DAQHit const& exit_hit)
152  {
153  auto contained_last =
154  find(contained_first, hitmap["x_contained"].cend(), exit_hit.TDC());
155 
156  if (std::distance(contained_first, contained_last) <
158  return false;
159 
160  Track track(*entry_hit, exit_hit);
161 
162  // std::cout << "\nMF (Track): (distance, beta, good) = ("
163  // << track.distance_in_cm() << ", "
164  // << track.beta() << ", "
165  // << std::boolalpha << good_track(track) << ")" << std::endl;
166  // ddthelpers::print(*entry_hit, event_time_);
167  // ddthelpers::print(exit_hit, event_time_);
168 
169  if (!good_track(track))
170  return false;
171 
172  if (++n_good_track % sparsification_factor_)
173  return false;
174 
175  std::set<int> plane_set = {entry_hit->Plane().val};
176  std::set<int> cell_set = {entry_hit->Cell().val};
177 
178  unsigned long long t_max_plane_gap = track.start_ddt().TDC().val;
179  if (track.slope_plane_is_valid())
180  t_max_plane_gap += abs(track.slope_time_plane()) * max_plane_gap_cut_;
181 
182  unsigned long long t_max_cell_gap = track.start_ddt().TDC().val;
183  if (track.slope_cell_is_valid())
184  t_max_cell_gap += abs(track.slope_time_cell()) * max_cell_gap_cut_;
185  unsigned long long t_max_gap = std::min(t_max_plane_gap, t_max_cell_gap);
186 
187  unsigned n_contained_hit = 0;
188  bool contains_gap_exceeding_cut = std::any_of
189  (contained_first, contained_last,
190  [&](novaddt::DAQHit const& hit)
191  {
192  if (hit.TDC().val > t_max_gap)
193  {
194  int64_t time_difference =
195  static_cast<int64_t>(hit.TDC().val) -
196  static_cast<int64_t>(track.start_ddt().TDC().val);
197 
198  if (hit.TDC().val > t_max_cell_gap)
199  {
200  double expected_cell =
201  static_cast<double>(track.start_ddt().Cell().val);
202  if (track.slope_time_is_valid())
203  expected_cell += track.slope_cell_time() * time_difference;
204 
205  if (gap_exceeds_cut(cell_set, expected_cell, max_cell_gap_cut_))
206  return true;
207  }
208 
209  if (hit.TDC().val > t_max_plane_gap)
210  {
211  double expected_plane =
212  static_cast<double>(track.start_ddt().Plane().val);
213  if (track.slope_time_is_valid())
214  expected_plane += track.slope_plane_time() * time_difference;
215 
216  if (gap_exceeds_cut
217  (plane_set, expected_plane, max_plane_gap_cut_))
218  return true;
219  }
220  }
221 
222  if (hit_is_in_time_with_road(hit, track))
223  {
224  cell_set.insert(hit.Cell().val);
225  plane_set.insert(hit.Plane().val);
226  }
227 
228  return false;
229  });
230 
231  // std::cout << "MF: contains_gap_exceeding_cut_ = "
232  // << std::boolalpha << contains_gap_exceeding_cut
233  // << std::endl;
234 
235  if (contains_gap_exceeding_cut)
236  return false;
237 
238  // std::cout << "MF: contains_gap_exceeding_cut_ = "
239  // << std::boolalpha << contains_gap_exceeding_cut
240  // << std::endl;
241 
242  plane_set.insert(exit_hit.Plane().val);
243  // std::cout << "MF: (plane_set_size, plane_max_gap) = ("
244  // << plane_set.size()
245  // << ", " << find_max_gap(plane_set)
246  // << ")" << std::endl;
247 
249  return false;
250 
251  if (find_max_gap(plane_set) > max_plane_gap_cut_)
252  return false;
253 
254  cell_set.insert(exit_hit.Cell().val);
255  // std::cout << "MF: (cell_set_size, cell_max_gap) = ("
256  // << cell_set.size()
257  // << ", " << find_max_gap(cell_set)
258  // << ")" << std::endl;
259 
260  if (find_max_gap(cell_set) > max_cell_gap_cut_)
261  return false;
262 
263  result = true;
264  trigger_entry_hit = *entry_hit;
265  trigger_exit_hit = exit_hit;
266 
267  std::cout << "MF: Trigger Decision Hits:" << std::endl;
268  ddthelpers::print(trigger_entry_hit, event_time_);
269  ddthelpers::print(trigger_exit_hit, event_time_);
270  std::cout << "MF: Trigger Decision Track:"
271  << "\nMF: (distance, velocity) = ("
272  << track.distance() << ", "
273  << track.beta() << ")" << std::endl;
274 
275  return true;
276  });
277  }
278 
279  // std::cout << "MF: n_good_track = " << n_good_track << std::endl;
280 
281 
282  if (result)
283  {
284  std::cout << "MF: Trigger Decision!" << std::endl;
285 
286  auto t0 = trigger_entry_hit.TDC().val;
287  auto dt = trigger_exit_hit.TDC().val - t0;
288 
289  trigger_decisions_.
290  emplace_back(t0 - 2 * trigger_time_buffer_,
291  dt + 4 * trigger_time_buffer_,
293 
294  essential_hit_lists_["mf_entry_exit"].push_back(trigger_entry_hit);
295  essential_hit_lists_["mf_entry_exit"].push_back(trigger_exit_hit);
296  }
297 
298  return result;
299 }
300 
301 
302 
305  std::map<std::string, novaddt::HitList> & hitmap) const
306 {
307  auto bound = std::stable_partition
308  (hits.begin(), hits.end(),
309  [](novaddt::DAQHit const& h)
310  { return h.View().val == daqchannelmap::X_VIEW; });
311 
312  hitmap["x"] = novaddt::HitList(std::distance(hits.begin(), bound));
313  std::move(hits.begin(), bound, hitmap["x"].begin());
314 
315  hitmap["y"] = novaddt::HitList(std::distance(bound, hits.end()));
316  std::move(bound, hits.end(), hitmap["y"].begin());
317 }
318 
319 
320 
322 (std::map<std::string, novaddt::HitList> & hitmap) const
323 {
324  for (auto const& x_hit : hitmap["x"])
325  {
326  auto y_first =
327  find(hitmap["y"],
329  auto y_last =
330  find(hitmap["y"],
332 
333  auto y_match = std::find_if
334  (y_first, y_last, [&](novaddt::DAQHit const& y_hit)
335  {
336  return hits_are_view_matched(x_hit, y_hit);
337  });
338 
339  if (y_match != y_last)
340  {
341  // we found a matched hit
342  hitmap["x_contained"].push_back(x_hit);
343 
344  // is this hit on the surface?
345  if (hit_is_on_surface(x_hit))
346  {
347  hitmap["x_surface"].push_back(x_hit);
348  } else {
349  // now we need to look for additional matched y-hits and check whether
350  // they are on the surface
351  auto surface_y_match = std::find_if
352  (y_match, y_last, [&](novaddt::DAQHit const& y_hit)
353  {
354  if (!hits_are_view_matched(x_hit, y_hit))
355  return false;
356 
357  if (!hit_is_on_surface(y_hit))
358  return false;
359 
360  return true;
361  });
362 
363  if (surface_y_match != y_last)
364  hitmap["x_surface"].push_back(x_hit);
365  }
366  }
367  }
368 }
369 
370 
371 
373 (novaddt::DAQHit const& x_hit, novaddt::DAQHit const& y_hit) const
374 {
375  int plane_difference =
376  static_cast<int>(x_hit.Plane().val) - static_cast<int>(y_hit.Plane().val);
377 
378  if (std::abs(plane_difference) <
380  return true;
381 
382  return false;
383 }
384 
385 
386 
388 (novaddt::DAQHit const& hit) const
389 {
391  return true;
392 
394  return true;
395 
397  return true;
398 
400  return true;
401 
402  return false;
403 }
404 
405 
406 
408 {
410  return false;
411 
412  if (track.velocity() < track_min_beta_)
413  return false;
414 
415  if (track.velocity() > track_max_beta_)
416  return false;
417 
418  // For simplicity, let's punt on all events that cross fewer planes
419  // than the maximum gap cut:
420  if (abs(track.plane_difference()) <= max_plane_gap_cut_)
421  return false;
422 
423  return true;
424 }
425 
426 
427 
429 (novaddt::DAQHit const& hit, Track const& track) const
430 {
431  if (!hit_is_between_track_end_points(hit, track))
432  return false;
433 
434  // If the plane of the hit and track end points are all equal,
435  // the hit is on the road.
436  if (!track.slope_plane_is_valid())
437  return true;
438 
439  double road_center_cell = track.slope_cell_plane() *
440  (static_cast<double>(hit.Plane().val) -
441  static_cast<double>(track.start_ddt().Plane().val)) +
442  static_cast<double>(track.start_ddt().Cell().val);
443 
444  if (hit.Cell().val >
445  road_center_cell + novaddt::smt::Constants::ROAD_HALF_WIDTH)
446  return false;
447 
448  if (hit.Cell().val <
449  road_center_cell - novaddt::smt::Constants::ROAD_HALF_WIDTH)
450  return false;
451 
452  return true;
453 }
454 
455 
456 
458 (novaddt::DAQHit const& hit, Track const& track) const
459 {
460  if (!hit_is_on_road(hit, track))
461  return false;
462 
463  // If the plane of the hit and track end points are all equal,
464  // we cannot calculate the time_slope and make no statement on the time.
465  if (!track.slope_plane_is_valid())
466  return false;
467 
468 
469  int64_t relative_hit_time =
470  static_cast<int64_t>(hit.TDC().val) -
471  static_cast<int64_t>(track.start_ddt().TDC().val);
472 
473  double expected_time = track.slope_time_plane() *
474  (static_cast<double>(hit.Plane().val) -
475  static_cast<double>(track.start_ddt().Plane().val));
476 
477  if (relative_hit_time >
479  return false;
480 
481  if (relative_hit_time <
483  return false;
484 
485  return true;
486 }
487 
488 
489 
491 (std::set<int> collection, int const& object, int const& max_gap_cut) const
492 {
493  // Note that we copy the collection in here on purpose, so that we do not
494  // accidentally alter it.
495  collection.insert(object);
496 
497  // std::cout << "MF: Collection = ";
498  // for (auto const& element : collection)
499  // std::cout << element << " ";
500  // std::cout << std::endl;
501  // std::cout << "MF: (max_gap, max_gap_cut) = ("
502  // << find_max_gap(collection)
503  // << ", " << max_gap_cut
504  // << ")" << std::endl;
505 
506  // We need to encode the max_plane cut differently here, so as to keep
507  // this method as general as possible.
508  if (find_max_gap(collection) > max_gap_cut)
509  return true;
510 
511  return false;
512 }
513 
514 
515 
517 (std::set<int> const& collection) const
518 {
519  int result = -9999;
520 
521  std::adjacent_find
522  (collection.begin(), collection.end(), [&](int object, int next_object)
523  {
524  int gap = next_object - object;
525  result = std::max(gap, result);
526 
527  // This return is meaningless, but a boolean return is required by
528  // adjacent find. We choose false, so the loop goes until the end.
529  return false;
530  });
531 
532  return result;
533 }
534 
535 
536 
538 (novaddt::DAQHit const& hit, Track const& track) const
539 {
540  if (hit.Plane().val < track.min_plane())
541  return false;
542 
543  if (hit.Plane().val > track.max_plane())
544  return false;
545 
546  if (hit.Cell().val < track.min_cell())
547  return false;
548 
549  if (hit.Cell().val > track.max_cell())
550  return false;
551 
552  return true;
553 }
554 
555 
556 
557 
558 novaddt::HitList::const_iterator mono::SlowMonopoleTrigger::find
559 (novaddt::HitList::const_iterator begin,
560  novaddt::HitList::const_iterator end, novaddt::TDC tdc) const
561 {
562  return std::lower_bound(begin, end, tdc,
564 }
565 
566 
567 
568 novaddt::HitList::const_iterator mono::SlowMonopoleTrigger::find
569 (novaddt::HitList const& hits, novaddt::TDC tdc) const
570 {
571  return find(hits.cbegin(), hits.cend(), tdc);
572 }
double beta() const
Definition: Track.cxx:129
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
const unsigned ROAD_HALF_WIDTH
Definition: Constants.h:37
bool hit_is_in_time_with_road(novaddt::DAQHit const &hit, Track const &track) const
bool hits_are_view_matched(novaddt::DAQHit const &x_hit, novaddt::DAQHit const &y_hit) const
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
novaddt::DAQHit start_ddt() const
Definition: Track.cxx:316
std::vector< novaddt::TriggerDecision > trigger_decisions() const
value_type val
Definition: BaseProducts.h:34
SlowMonopoleTrigger(novaddt::HitList const &h, unsigned long long const &event_time, Configuration const &c)
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
std::vector< DAQHit > HitList
Definition: HitList.h:15
value_type val
Definition: BaseProducts.h:109
bool hit_is_on_surface(novaddt::DAQHit const &hit) const
void find_view_matched_hits(std::map< std::string, novaddt::HitList > &hitmap) const
Definition: event.h:19
OStream cerr
Definition: OStream.cxx:7
const unsigned SURFACE_MAX_CELL
Definition: Constants.h:24
bool parameter_exists(std::string const name) const
unsigned distance(const T &t1, const T &t2)
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
std::map< std::string, novaddt::HitList > essential_hit_lists() const
const double ROAD_TIME_VARIATION
Definition: Constants.h:38
novaddt::HitList::const_iterator find(novaddt::HitList::const_iterator begin, novaddt::HitList::const_iterator end, novaddt::TDC tdc) const
double slope_plane_time() const
Definition: Track.cxx:202
int find_max_gap(std::set< int > const &collection) const
unsigned max_cell() const
Definition: Track.cxx:266
bool good_track(Track const &track) const
Definition: Cand.cxx:23
void hits()
Definition: readHits.C:15
void print(novaddt::DAQHit const &h, unsigned long long const &event_time=0)
Definition: DDTHelpers.cxx:84
unsigned trigger_decision_time_buffer() const
Identifier for the X measuring view of the detector (top)
const unsigned TRACK_MIN_UNIQUE_PLANES
Definition: Constants.h:34
const unsigned VIEW_MATCH_PLANE_CUT
Definition: Constants.h:31
const unsigned SURFACE_MIN_PLANE
Definition: Constants.h:27
value_type val
Definition: BaseProducts.h:84
const double TRACK_MIN_DISTANCE
Definition: Constants.h:35
OStream cout
Definition: OStream.cxx:6
bool hit_is_on_road(novaddt::DAQHit const &hit, Track const &track) const
bool slope_cell_is_valid() const
Definition: Track.cxx:150
void split_by_view(novaddt::HitList hits, std::map< std::string, novaddt::HitList > &hitmap) const
double parameter(std::string const name)
double slope_cell_plane() const
Definition: Track.cxx:157
double slope_time_plane() const
Definition: Track.cxx:172
double velocity() const
Definition: Track.cxx:122
Definition: structs.h:12
bool gap_exceeds_cut(std::set< int > collection, int const &object, int const &max_gap_cut) const
double slope_time_cell() const
Definition: Track.cxx:217
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
assert(nhit_max >=nhit_nbins)
bool hit_is_between_track_end_points(novaddt::DAQHit const &hit, Track const &track) const
unsigned min_cell() const
Definition: Track.cxx:259
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
Float_t e
Definition: plot.C:35
double distance() const
Definition: Track.cxx:108
double slope_cell_time() const
Definition: Track.cxx:187
Float_t track
Definition: plot.C:35
const unsigned TRACK_MIN_HITS
Definition: Constants.h:33
unsigned min_plane() const
Definition: Track.cxx:238
enum BeamMode string