Trigger.cxx
Go to the documentation of this file.
2 
6 
7 #include <algorithm>
8 #include <bitset>
9 #include <chrono>
10 #include <iostream>
11 
12 using namespace std::chrono;
13 
14 
15 
17  max_cell_gap_cut_ (p.get<int>("max_cell_gap_cut")),
18  max_plane_gap_cut_ (p.get<int>("max_plane_gap_cut")),
19  sparsification_factor_ (p.get<unsigned>("sparsification_factor")),
20  track_max_beta_ (p.get<double>("track_max_beta")),
21  track_min_beta_ (p.get<double>("track_min_beta")),
22  trigger_time_buffer_ (p.get<unsigned>("trigger_time_buffer")),
23  trigger_timeout_ (p.get<double>("trigger_timeout")),
24  trigger_use_timeout_ (p.get<bool>("trigger_use_timeout")),
25  trigger_write_all_timeout_data_(p.get<bool>("trigger_write_all_timeout_data"))
26 {
27 }
28 
29 
30 
31 std::vector<novaddt::TriggerDecision>
33 {
34  return trigger_decisions_;
35 }
36 
37 
38 
41 {
42  trigger_decisions_.clear();
43  std::bitset<Extra_Info::BITSET_SIZE> extra_info;
44  extra_info.set(Extra_Info::CHECK_BIT, true);
45 
46  auto t_start = high_resolution_clock::now();
47  bool result = false;
48  novaddt::DAQHit trigger_entry_hit, trigger_exit_hit;
49 
50  std::map<std::string, novaddt::HitList> hitmap;
51  split_by_view(*hits, hitmap);
52 
53  if (hitmap["x"].size() < Constants::TRACK_MIN_HITS ||
54  hitmap["y"].size() < Constants::TRACK_MIN_HITS)
55  return false;
56 
57  hitmap["x_surface"];
58  hitmap["x_contained"];
59 
60  find_view_matched_hits(hitmap);
61 
62  if(hitmap["x_contained"].size() < Constants::TRACK_MIN_HITS)
63  return false;
64 
65  unsigned n_good_track = 0;
66  for (auto entry_hit = hitmap["x_surface"].cbegin();
67  entry_hit != hitmap["x_surface"].cend() && !result;
68  ++entry_hit)
69  {
70  // we add 1, so that we avoid combinations between the same hit
71  novaddt::HitList::const_reverse_iterator until(entry_hit + 1);
72 
73  auto contained_first = find(hitmap["x_contained"], entry_hit->TDC());
74 
75  std::find_if(hitmap["x_surface"].crbegin(), until,
76  [&](novaddt::DAQHit const& exit_hit)
77  {
78  auto contained_last =
79  find(contained_first, hitmap["x_contained"].cend(), exit_hit.TDC());
80 
81  if (std::distance(contained_first, contained_last) <
83  return false;
84 
85  Track track(*entry_hit, exit_hit);
86 
87  if (!good_track(track))
88  return false;
89 
90  if (++n_good_track % sparsification_factor_)
91  return false;
92 
93  auto time_spent = duration_cast<milliseconds>
94  (high_resolution_clock::now() - t_start).count();
95 
97  time_spent > trigger_timeout_)
98  {
99  result = true;
100  trigger_entry_hit = *entry_hit;
101 
103  trigger_exit_hit = hitmap["x_surface"].back();
104  else
105  trigger_exit_hit = *entry_hit;
106 
107  extra_info.set(Extra_Info::TIMEOUT_BIT, true);
108 
109  return true;
110  }
111 
112  std::set<int> plane_set = {entry_hit->Plane().val};
113  std::set<int> cell_set = {entry_hit->Cell().val};
114 
115  uint64_t t_max_plane_gap = track.start().TDC().val;
116  if (track.slope_plane_is_valid())
117  t_max_plane_gap += abs(track.slope_time_plane()) * max_plane_gap_cut_;
118 
119  uint64_t t_max_cell_gap = track.start().TDC().val;
120  if (track.slope_cell_is_valid())
121  t_max_cell_gap += abs(track.slope_time_cell()) * max_cell_gap_cut_;
122  uint64_t t_max_gap = std::min(t_max_plane_gap, t_max_cell_gap);
123 
124  //unsigned n_contained_hit = 0;
125  bool contains_gap_exceeding_cut = std::any_of
126  (contained_first, contained_last,
127  [&](novaddt::DAQHit const& hit)
128  {
129  if (hit.TDC().val > t_max_gap)
130  {
131  int64_t time_difference =
132  static_cast<int64_t>(hit.TDC().val) -
133  static_cast<int64_t>(track.start().TDC().val);
134 
135  if (hit.TDC().val > t_max_cell_gap)
136  {
137  double expected_cell =
138  static_cast<double>(track.start().Cell().val);
139  if (track.slope_time_is_valid())
140  expected_cell += track.slope_cell_time() * time_difference;
141 
142  if (gap_exceeds_cut(cell_set, expected_cell, max_cell_gap_cut_))
143  return true;
144  }
145 
146  if (hit.TDC().val > t_max_plane_gap)
147  {
148  double expected_plane =
149  static_cast<double>(track.start().Plane().val);
150  if (track.slope_time_is_valid())
151  expected_plane += track.slope_plane_time() * time_difference;
152 
153  if (gap_exceeds_cut
154  (plane_set, expected_plane, max_plane_gap_cut_))
155  return true;
156  }
157  }
158 
159  if (hit_is_in_time_with_road(hit, track))
160  {
161  cell_set.insert(hit.Cell().val);
162  plane_set.insert(hit.Plane().val);
163  }
164 
165  return false;
166  });
167 
168  if (contains_gap_exceeding_cut)
169  return false;
170 
171  plane_set.insert(exit_hit.Plane().val);
172  if (plane_set.size() < Constants::TRACK_MIN_UNIQUE_PLANES)
173  return false;
174 
175  if (find_max_gap(plane_set) > max_plane_gap_cut_)
176  return false;
177 
178  cell_set.insert(exit_hit.Cell().val);
179  if (find_max_gap(cell_set) > max_cell_gap_cut_)
180  return false;
181 
182  result = true;
183  trigger_entry_hit = *entry_hit;
184  trigger_exit_hit = exit_hit;
185 
186  return true;
187  });
188  }
189 
190  if (result)
191  {
192  auto t_begin = trigger_entry_hit.TDC().val - trigger_time_buffer_;
193  auto t_end = trigger_exit_hit.TDC().val + trigger_time_buffer_;
194 
196  emplace_back(t_begin, t_end - t_begin,
198  extra_info.to_ulong());
199  }
200 
201  return result;
202 }
203 
204 
205 
207 (novaddt::HitList hits, std::map<std::string, novaddt::HitList> & hitmap) const
208 {
209  auto bound = std::stable_partition
210  (hits.begin(), hits.end(),
211  [](novaddt::DAQHit const& h)
212  { return h.View().val == daqchannelmap::X_VIEW; });
213 
214  hitmap["x"] = novaddt::HitList(std::distance(hits.begin(), bound));
215  std::move(hits.begin(), bound, hitmap["x"].begin());
216 
217  hitmap["y"] = novaddt::HitList(std::distance(bound, hits.end()));
218  std::move(bound, hits.end(), hitmap["y"].begin());
219 }
220 
221 
222 
224 (std::map<std::string, novaddt::HitList> & hitmap) const
225 {
226  for (auto const& x_hit : hitmap["x"])
227  {
228  auto y_first =
229  find(hitmap["y"], x_hit.TDC().val - Constants::VIEW_MATCH_TDC_CUT);
230  auto y_last =
231  find(hitmap["y"], x_hit.TDC().val + Constants::VIEW_MATCH_TDC_CUT);
232 
233  auto y_match = std::find_if
234  (y_first, y_last, [&](novaddt::DAQHit const& y_hit)
235  { return hits_are_view_matched(x_hit, y_hit); });
236 
237  if (y_match != y_last)
238  {
239  // we found a matched hit
240  hitmap["x_contained"].push_back(x_hit);
241 
242  // is this hit on the surface?
243  if (hit_is_on_surface(x_hit))
244  {
245  hitmap["x_surface"].push_back(x_hit);
246  } else {
247  // now we need to look for additional matched y-hits and check whether
248  // they are on the surface
249  auto surface_y_match = std::find_if
250  (y_match, y_last, [&](novaddt::DAQHit const& y_hit)
251  {
252  if (!hits_are_view_matched(x_hit, y_hit))
253  return false;
254 
255  if (!hit_is_on_surface(y_hit))
256  return false;
257 
258  return true;
259  });
260 
261  if (surface_y_match != y_last)
262  hitmap["x_surface"].push_back(x_hit);
263  }
264  }
265  }
266 }
267 
268 
269 
271 (novaddt::DAQHit const& x_hit, novaddt::DAQHit const& y_hit) const
272 {
273  int plane_difference =
274  static_cast<int>(x_hit.Plane().val) - static_cast<int>(y_hit.Plane().val);
275 
276  if (std::abs(plane_difference) < Constants::VIEW_MATCH_PLANE_CUT)
277  return true;
278 
279  return false;
280 }
281 
282 
283 
285 (novaddt::DAQHit const& hit) const
286 {
288  return true;
289 
291  return true;
292 
294  return true;
295 
297  return true;
298 
299  return false;
300 }
301 
302 
303 
305 {
307  return false;
308 
309  if (track.beta() < track_min_beta_)
310  return false;
311 
312  if (track.beta() > track_max_beta_)
313  return false;
314 
315  if (abs(track.plane_difference()) <= max_plane_gap_cut_)
316  return false;
317 
318  return true;
319 }
320 
321 
322 
324 (novaddt::DAQHit const& hit, Track const& track) const
325 {
326  if (!hit_is_between_track_end_points(hit, track))
327  return false;
328 
329  // If the plane of the hit and track end points are all equal,
330  // the hit is on the road.
331  if (!track.slope_plane_is_valid())
332  return true;
333 
334  double road_center_cell = track.slope_cell_plane() *
335  (static_cast<double>(hit.Plane().val) -
336  static_cast<double>(track.start().Plane().val)) +
337  static_cast<double>(track.start().Cell().val);
338 
339  if (hit.Cell().val > road_center_cell + Constants::ROAD_HALF_WIDTH)
340  return false;
341 
342  if (hit.Cell().val < road_center_cell - Constants::ROAD_HALF_WIDTH)
343  return false;
344 
345  return true;
346 }
347 
348 
349 
351 (novaddt::DAQHit const& hit, Track const& track) const
352 {
353  if (!hit_is_on_road(hit, track))
354  return false;
355 
356  // If the plane of the hit and track end points are all equal,
357  // we cannot calculate the time_slope and make no statement on the time.
358  if (!track.slope_plane_is_valid())
359  return false;
360 
361 
362  int64_t relative_hit_time =
363  static_cast<int64_t>(hit.TDC().val) -
364  static_cast<int64_t>(track.start().TDC().val);
365 
366  double expected_time = track.slope_time_plane() *
367  (static_cast<double>(hit.Plane().val) -
368  static_cast<double>(track.start().Plane().val));
369 
370  if (relative_hit_time > expected_time + Constants::ROAD_TIME_VARIATION)
371  return false;
372 
373  if (relative_hit_time < expected_time - Constants::ROAD_TIME_VARIATION)
374  return false;
375 
376  return true;
377 }
378 
379 
380 
382 (std::set<int> collection, int const& object, int const& max_gap_cut) const
383 {
384  // Note that we copy the collection in here on purpose, so that we do not
385  // accidentally alter it.
386  collection.insert(object);
387 
388  // We need to encode the max_plane cut differently here, so as to keep
389  // this method as general as possible.
390  if (find_max_gap(collection) > max_gap_cut)
391  return true;
392 
393  return false;
394 }
395 
396 
397 
399 (std::set<int> const& collection) const
400 {
401  int result = -9999;
402 
403  std::adjacent_find
404  (collection.begin(), collection.end(), [&](int object, int next_object)
405  {
406  int gap = next_object - object;
407  result = std::max(gap, result);
408 
409  // This return is meaningless, but a boolean return is required by
410  // adjacent find. We choose false, so the loop goes until the end.
411  return false;
412  });
413 
414  return result;
415 }
416 
417 
418 
420 (novaddt::DAQHit const& hit, Track const& track) const
421 {
422  if (hit.Plane().val < track.min_plane())
423  return false;
424 
425  if (hit.Plane().val > track.max_plane())
426  return false;
427 
428  if (hit.Cell().val < track.min_cell())
429  return false;
430 
431  if (hit.Cell().val > track.max_cell())
432  return false;
433 
434  return true;
435 }
436 
437 
438 
439 novaddt::HitList::const_iterator novaddt::smt::Trigger::find
440 (novaddt::HitList::const_iterator begin,
441  novaddt::HitList::const_iterator end, novaddt::TDC tdc) const
442 {
443  return
444  std::lower_bound(begin, end, tdc, novaddt::CompareDAQHit<novaddt::TDC>());
445 }
446 
447 
448 
449 novaddt::HitList::const_iterator novaddt::smt::Trigger::find
450 (novaddt::HitList const& hits, novaddt::TDC tdc) const
451 {
452  return find(hits.cbegin(), hits.cend(), tdc);
453 }
T max(const caf::Proxy< T > &a, T b)
double track_max_beta_
Definition: Trigger.h:79
const unsigned ROAD_HALF_WIDTH
Definition: Constants.h:37
double slope_cell_plane() const
Definition: Track.cxx:127
bool hit_is_in_time_with_road(novaddt::DAQHit const &hit, Track const &track) const
Definition: Trigger.cxx:351
unsigned sparsification_factor_
Definition: Trigger.h:78
bool hit_is_on_road(novaddt::DAQHit const &hit, Track const &track) const
Definition: Trigger.cxx:324
value_type val
Definition: BaseProducts.h:34
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
bool hits_are_view_matched(novaddt::DAQHit const &x_hit, novaddt::DAQHit const &y_hit) const
Definition: Trigger.cxx:271
bool good_track(Track const &track) const
Definition: Trigger.cxx:304
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
std::vector< DAQHit > HitList
Definition: HitList.h:15
novaddt::DAQHit start() const
Definition: Track.cxx:202
const char * p
Definition: xmltok.h:285
value_type val
Definition: BaseProducts.h:109
unsigned min_cell() const
Definition: Track.cxx:237
Definition: event.h:19
bool hit_is_between_track_end_points(novaddt::DAQHit const &hit, Track const &track) const
Definition: Trigger.cxx:420
const unsigned SURFACE_MAX_CELL
Definition: Constants.h:24
void abs(TH1 *hist)
bool slope_cell_is_valid() const
Definition: Track.cxx:120
unsigned distance(const T &t1, const T &t2)
std::vector< novaddt::TriggerDecision > trigger_decisions_
Definition: Trigger.h:86
double distance() const
Definition: Track.cxx:78
float abs(float number)
Definition: d0nt_math.hpp:39
double slope_time_plane() const
Definition: Track.cxx:142
const unsigned SURFACE_MAX_PLANE
Definition: Constants.h:28
const double ROAD_TIME_VARIATION
Definition: Constants.h:38
double beta() const
Definition: Track.cxx:99
int find_max_gap(std::set< int > const &collection) const
Definition: Trigger.cxx:399
unsigned max_plane() const
Definition: Track.cxx:223
bool slope_plane_is_valid() const
Definition: Track.cxx:106
void find_view_matched_hits(std::map< std::string, novaddt::HitList > &hitmap) const
Definition: Trigger.cxx:224
void hits()
Definition: readHits.C:15
bool hit_is_on_surface(novaddt::DAQHit const &hit) const
Definition: Trigger.cxx:285
Identifier for the X measuring view of the detector (top)
const unsigned TRACK_MIN_UNIQUE_PLANES
Definition: Constants.h:34
void split_by_view(novaddt::HitList hits, std::map< std::string, novaddt::HitList > &hitmap) const
Definition: Trigger.cxx:207
const unsigned VIEW_MATCH_PLANE_CUT
Definition: Constants.h:31
double slope_plane_time() const
Definition: Track.cxx:172
std::vector< novaddt::TriggerDecision > trigger_decisions() const
Definition: Trigger.cxx:32
bool run_algorithm(art::Handle< novaddt::HitList > const &hits)
Definition: Trigger.cxx:40
double track_min_beta_
Definition: Trigger.h:80
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
bool trigger_write_all_timeout_data_
Definition: Trigger.h:84
const unsigned CHECK_BIT
Definition: Constants.h:45
double slope_cell_time() const
Definition: Track.cxx:157
bool slope_time_is_valid() const
Definition: Track.cxx:113
bool gap_exceeds_cut(std::set< int > collection, int const &object, int const &max_gap_cut) const
Definition: Trigger.cxx:382
bool trigger_use_timeout_
Definition: Trigger.h:83
unsigned min_plane() const
Definition: Track.cxx:216
int plane_difference() const
Definition: Track.cxx:230
Definition: structs.h:12
double trigger_timeout_
Definition: Trigger.h:82
double slope_time_cell() const
Definition: Track.cxx:187
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
novaddt::HitList::const_iterator find(novaddt::HitList::const_iterator begin, novaddt::HitList::const_iterator end, novaddt::TDC tdc) const
Definition: Trigger.cxx:440
const unsigned TIMEOUT_BIT
Definition: Constants.h:46
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
unsigned trigger_time_buffer_
Definition: Trigger.h:81
unsigned max_cell() const
Definition: Track.cxx:244
Float_t track
Definition: plot.C:35
const unsigned TRACK_MIN_HITS
Definition: Constants.h:33