HoughTrackMaker.cxx
Go to the documentation of this file.
2 #include <TVector2.h>
3 #include <utility>
4 
5 #include <iostream>
6 #include <set>
7 
9  : _tdcSigma (p.get<uint64_t>("max_tdc_separation"))
10  , _planeSigma (p.get<int>("max_plane_separation"))
11  , _cellSigma (p.get<int>("max_cell_separation"))
12  , _vCut (p.get<float>("velocity_cut"))
13  , _minPointsToTag (p.get<int>("minimum_points_to_tag"))
14  , _minHitsPerTrack (p.get<int>("minimum_hits_per_track"))
15  , _maxHitsInSlice (p.get<int>("maximum_hits_per_slice"))
16 {
17  tracks_ =
18  std::unique_ptr<std::vector<Track> >(new std::vector<Track>);
19  track_hits_ =
20  std::unique_ptr<std::vector<HitList> >(new std::vector<HitList>);
21  track_slice_.clear();
22 }
23 
24 
25 
28  : _tdcSigma(other._tdcSigma)
29  , _planeSigma(other._planeSigma)
30  , _cellSigma(other._cellSigma)
31  , _vCut(other._vCut)
35 {
36  tracks_ =
37  std::unique_ptr<std::vector<Track> >(new std::vector<Track>);
38  track_hits_ =
39  std::unique_ptr<std::vector<HitList> >(new std::vector<HitList>);
40  track_slice_.clear();
41 }
42 
43 
44 
46 {
47  for (int c_i = 0; c_i != fcGrid_; ++c_i)
48  {
49  for (int d_i = 0; d_i != fdGrid_; ++d_i)
50  {
51  for (auto pair : pair_pool_[0][c_i][d_i])
52  delete pair;
53 
54  for (auto pair : pair_pool_[1][c_i][d_i])
55  delete pair;
56  }
57  }
58 }
59 
60 
61 
63 {
64  for (auto const& slice : slices)
65  do_like_zukai(slice);
66 }
67 
68 
69 
71 {
72  for (auto const& slice : slices)
73  do_like_zukai(slice);
74 }
75 
76 
77 
79 {
80  tracks_->insert(tracks_->end(),
81  other.tracks_->begin(), other.tracks_->end());
82 
83  track_hits_->insert(track_hits_->end(),
84  other.track_hits_->begin(), other.track_hits_->end());
85 
86  track_slice_.insert(track_slice_.end(),
87  other.track_slice_.begin(), other.track_slice_.end());
88 }
89 
90 
91 
94  art::Event & e,
97 {
98  if (tracks_->size() != track_hits_->size() ||
99  tracks_->size() != track_slice_.size())
100  {
101  std::cerr << "\n\tThe number of tracks, associated track hit lists, and "
102  << "slice indexes is not the same!"
103  << "\n\t\ttracks: " << tracks_->size()
104  << "\n\t\thit lists: " << track_hits_->size()
105  << "\n\t\tslice indexes: " << track_slice_.size()
106  << std::endl;
107  assert(false);
108  }
109 
110  for (unsigned index = 0; index < tracks_->size(); ++index)
111  {
112  // get the Ptr for one of the tracks
113  art::ProductID t_id = filter.getProductID<std::vector<Track> >();
114  art::Ptr<Track> t_ptr(t_id, index, e.productGetter(t_id));
115 
116  // get the Ptr for one of the track_hits
117  art::ProductID th_id = filter.getProductID<std::vector<HitList> >();
118  art::Ptr<HitList> th_ptr(th_id, index, e.productGetter(th_id));
119 
120  // Association I
121  assn.addSingle(t_ptr, th_ptr);
122 
123  // Association II
124  assn_b.addSingle(t_ptr, track_slice_[index]);
125  }
126 }
127 
128 
129 
131 {
132  e.put(std::move(tracks_));
133  e.put(std::move(track_hits_));
134 }
135 
136 
137 
139 {
140 
141  // STAGE 0: Sanity Checks
142  //======================================================================================================================
143  // Bail if the number of hits in the slice are to high
144  if(slice->size() > _maxHitsInSlice){
145  return false;
146  }
147 
148  //PART 1: Preparations
149  //======================================================================================================================
150  for (int c_i = 0; c_i != fcGrid_; ++c_i)
151  {
152  for (int d_i = 0; d_i != fdGrid_; ++d_i)
153  {
154  flag_[0][c_i][d_i] = false;
155  flag_[1][c_i][d_i] = false;
156 
157  pool_location_[0][c_i][d_i] = 0;
158  pool_location_[1][c_i][d_i] = 0;
159  }
160  }
161 
162  int n_tracks = 0;
163 
164  const float cos_max = 1;
165  const float rho_max = 828.6;
166  const float cos_wid = 2.0*cos_max/(fcGrid_ - 0.);
167  const float rho_wid = 2.0*rho_max/(fdGrid_ - 0.);
168 
169  int view_index = 0;
170  int c = 0;
171  int d = 0;
172 
173  std::vector<int> c_tag;
174  std::vector<int> d_tag;
175  std::vector<int> view_tag;
176 
177  //Before voting, the coordinate transform is to be done:
178  std::vector<TVector2> Coordinates;
179  for (auto hit = slice->begin(); hit != slice->end() ; ++hit) {
180  TVector2 coord(1.6795*(hit->Plane().val-480),
181  hit->Cell().val-191.5);
182 
183  Coordinates.emplace_back(coord);
184  }
185 
186  //======================================================================================================================
187  //PART 2: Voting inside this slice
188  //======================================================================================================================
189  unsigned idx_i=0;
190  unsigned idx_j=0;
191 
192  for (auto hit1 = slice->begin(); hit1 != slice->end() - 1; ++hit1)
193  {
194  // uint8_t view = hit1->View().val;
195  view_index = 0;
196  if (hit1->View().val != daqchannelmap::X_VIEW)
197  view_index = 1;
198 
199  idx_j=idx_i+1;
200  for (auto hit2 = hit1 + 1; hit2 != slice->end(); ++hit2)
201  {
202 
203  if (PairCutter(hit1,hit2)) {
204  ++idx_j;
205  continue;
206  }
207 
208  // Notice: this requires later delete the memory:
209  HoughPoint* hp_ptr(new HoughPoint(*hit1,Coordinates.at(idx_i),*hit2,Coordinates.at(idx_j)));
210  ++idx_j;
211  float const& vpro = fabs(hp_ptr->getVpro());
212  if (vpro > _vCut)
213  {
214  delete hp_ptr;
215  continue;
216  }
217 
218  float const& cosv = hp_ptr->getCosv();
219  float const& doca = hp_ptr->getDoca();
220 
221  //Figure out which grid this pair should be push back:
222  c = int((cosv + cos_max) / cos_wid);
223  d = int((doca + rho_max) / rho_wid);
224 
225  // Only discovered once:
226  // but this kind of labeling is having the danger of boundary problems:
227  if (c == fcGrid_)
228  c = fcGrid_ - 1;
229  if (d == fdGrid_)
230  d = fdGrid_ - 1;
231 
232  unsigned location = pool_location_[view_index][c][d];
233  unsigned size = pair_pool_[view_index][c][d].size();
234 
235  if (location < size)
236  *pair_pool_[view_index][c][d][location] = *hp_ptr;
237  else
238  pair_pool_[view_index][c][d].emplace_back(new HoughPoint(*hp_ptr));
239 
240  pool_location_[view_index][c][d]++;
241 
242  delete hp_ptr;
243 
244  if (flag_[view_index][c][d])
245  continue;
246  else if (location < _minPointsToTag)
247  continue;
248 
249  flag_[view_index][c][d]=true;
250 
251  //tag the index
252  c_tag.emplace_back(c);
253  d_tag.emplace_back(d);
254  view_tag.emplace_back(view_index);
255 
256  } // inner hit loop end
257  ++idx_i;
258  } // Voting Loop finished.
259 
260  unsigned ntag = c_tag.size();
261 
262  if (ntag == 0)
263  return false;
264  // else
265  // result = true;
266 
267  //======================================================================================================================
268  //Voting finished
269  //PART 3: Clustering in hough space
270  //======================================================================================================================
271  std::vector<std::vector<int> > tag_lists;
272  bool visit_flag[ntag];
273 
274  for (unsigned i = 0; i!= ntag; ++i) visit_flag[i] = 0;
275 
276  for (unsigned i = 0; i!= ntag; ++i) {
277  std::vector<int> tags;
278 
279  //whether this tag has been visited?
280  if (visit_flag[i]) continue;
281  visit_flag[i] = true;
282  tags.emplace_back(i);
283 
284  //Search for its neighbors:
285  for (unsigned j = i+1; j!= ntag; ++j) {
286  if (visit_flag[j]) continue;
287  //judge whether a neighbor:
288  if (view_tag[i]!=view_tag[j]) continue;
289 
290  //This configuration is for vertical tracks:
291  if (c_tag[i]<3 && c_tag[j]<3) {
292  if (abs(d_tag[i]+d_tag[j]-100)<3 || abs(c_tag[i]-c_tag[j])<3){
293  visit_flag[j] = true;
294  tags.emplace_back(j);
295  }
296  }
297 
298  else if (abs(c_tag[i]-c_tag[j])<3 && abs(d_tag[i]-d_tag[j])<3) {
299  visit_flag[j] = true;
300  tags.emplace_back(j);
301 
302  //2nd order combination:
303  for (unsigned k = j-1; k!=ntag;++k) {
304  if (visit_flag[k]) continue;
305  //judge whether a neighbor:
306  if (view_tag[k]!=view_tag[j]) continue;
307  if (abs(c_tag[k]-c_tag[j])<3 && abs(d_tag[k]-d_tag[j])<3) {
308  visit_flag[k] = true;
309  tags.emplace_back(k);
310  }
311  }
312  }
313  }
314  tag_lists.emplace_back(tags);
315  } // Clustering Loop finished
316 
317  //======================================================================================================================
318  //PART 4: Producing Tracks
319  //======================================================================================================================
320 
321  unsigned n_track_candidates = tag_lists.size();
322 
323  for (unsigned i=0; i != n_track_candidates; ++i)
324  {
325  std::set<const DAQHit*, hitcomp> track_candidate;
326  std::vector<const DAQHit*> test_vec;
327  int current_view = view_tag[tag_lists[i][0]];
328  for (unsigned j = 0; j != tag_lists[i].size(); ++j)
329  {
330  for (unsigned p = 0; p != pool_location_[current_view]
331  [c_tag[tag_lists[i][j]]]
332  [d_tag[tag_lists[i][j]]]; ++p)
333  {
334  HoughPoint* pair = pair_pool_[current_view]
335  [c_tag[tag_lists[i][j]]]
336  [d_tag[tag_lists[i][j]]][p];
337 
338  track_candidate.insert(&pair->getCh1());
339  track_candidate.insert(&pair->getCh2());
340  }
341  }
342 
343  if (track_candidate.size() < _minHitsPerTrack)
344  continue;
345 
346  // By now, the track candidate hit list is sorted by TDC!
347  // Prepare the hitlist:
348 
349  HitList this_hit_list;
350 
351  //Initialize these quantities:
352  TDC startT = (*track_candidate.begin())->TDC();
353  TDC endT = (*track_candidate.rbegin())->TDC();
354  float startZ = ((*track_candidate.begin())->Plane()).val;
355  float endZ = startZ;
356  float startV = ((*track_candidate.begin())->Cell()).val;
357  float endV = startV;
358 
359  //Determine start and end by Z:
360  if (abs(c_tag[tag_lists[i][0]] - 49) < 26)
361  {
362  for (auto const& hitptr : track_candidate)
363  {
364  this_hit_list.emplace_back(*hitptr);
365  if (hitptr->Plane().val < startZ)
366  {
367  startZ = hitptr->Plane().val;
368  startV = hitptr->Cell().val;
369  startT = hitptr->TDC();
370  }
371 
372  if (hitptr->Plane().val > endZ)
373  {
374  endZ = hitptr->Plane().val;
375  endV = hitptr->Cell().val;
376  endT = hitptr->TDC();
377  }
378  }
379  } // end of start determination by z
380 
381  //Determine start and end by V:
382  else {
383  for (auto const& hitptr : track_candidate )
384  {
385  this_hit_list.emplace_back(*hitptr);
386  if (hitptr->Cell().val < startV)
387  {
388  startZ = hitptr->Plane().val;
389  startV = hitptr->Cell().val;
390  startT = hitptr->TDC();
391  }
392 
393  if (hitptr->Cell().val > endV)
394  {
395  endZ = hitptr->Plane().val;
396  endV = hitptr->Cell().val;
397  endT = hitptr->TDC();
398  }
399  }
400  } // end of start determination by v
401 
402  uint8_t the_view = daqchannelmap::X_VIEW;
403  if (current_view == 1)
404  the_view = daqchannelmap::Y_VIEW;
405 
406 
407  novaddt::Track track2D (the_view,
408  0,
409  startV,
410  startZ,
411  startT,
412  endV,
413  endZ,
414  endT);
415 
416  tracks_->emplace_back(track2D);
417  track_hits_->emplace_back(this_hit_list);
418  track_slice_.emplace_back(slice);
419 
420  ++n_tracks;
421  } // end of track production
422  return true;
423 }
424 
425 
426 
428 (HitList::const_iterator hit_ptr1, HitList::const_iterator hit_ptr2) const
429 {
430  if (hit_ptr1->View().val != hit_ptr2->View().val) return true;
431  // possible for hits of same position but at different times
432  if ( (hit_ptr1->Plane() == hit_ptr2->Plane()) &&
433  (hit_ptr1->Cell() == hit_ptr2->Cell()) ) return true;
434  else if ((std::max(hit_ptr1->Plane().val,hit_ptr2->Plane().val) -
435  std::min(hit_ptr1->Plane().val,hit_ptr2->Plane().val)) > _planeSigma)
436  return true;
437 
438  else if ((std::max(hit_ptr1->Cell().val,hit_ptr2->Cell().val) -
439  std::min(hit_ptr1->Cell().val,hit_ptr2->Cell().val)) > _cellSigma)
440  return true;
441 
442  // Just this one line: partitioning by TDC!
443  else if ((std::max(hit_ptr1->TDC().val,hit_ptr2->TDC().val) -
444  std::min(hit_ptr1->TDC().val,hit_ptr2->TDC().val)) > _tdcSigma)
445  return true;
446 
447  else return false;
448 }
EDProductGetter const * productGetter(ProductID const) const
T max(const caf::Proxy< T > &a, T b)
void split(double tt, double *fr)
bool flag_[2][fcGrid_][fdGrid_]
float const & getCosv() const
Definition: HoughPoint.h:70
std::vector< art::Ptr< HitList > > track_slice_
#define location
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::unique_ptr< std::vector< Track > > tracks_
bool do_like_zukai(art::Ptr< HitList > slice)
std::vector< DAQHit > HitList
Definition: HitList.h:15
ProductID getProductID(std::string const &instanceName={}) const
Definition: EDFilter.h:131
const char * p
Definition: xmltok.h:285
float const & getDoca() const
Definition: HoughPoint.h:71
void create_associations(art::EDFilter const &filter, art::Event &e, art::Assns< novaddt::Track, novaddt::HitList > &assn, art::Assns< novaddt::Track, novaddt::HitList > &assn_b) const
DAQHit const & getCh2() const
Definition: HoughPoint.h:77
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
std::vector< HoughPoint * > pair_pool_[2][fcGrid_][fdGrid_]
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.
Identifier for the Y measuring view of the detector (side)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
TDC()=default
void join(HoughTrackMaker const &other)
Identifier for the X measuring view of the detector (top)
bool PairCutter(HitList::const_iterator hit_ptr1, HitList::const_iterator hit_ptr2) const
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
void test_vec()
void operator()(br_t const &slices)
void do_serial(art::PtrVector< HitList > const &slices)
unsigned pool_location_[2][fcGrid_][fdGrid_]
tbb::blocked_range< art::PtrVector< HitList >::const_iterator > br_t
static const int fdGrid_
Definition: structs.h:12
void addSingle(Ptr< left_t > const &left, Ptr< right_t > const &right, data_t const &data)
Definition: Assns.h:489
assert(nhit_max >=nhit_nbins)
void move_collections_into_event(art::Event &e)
std::unique_ptr< std::vector< HitList > > track_hits_
DAQHit const & getCh1() const
Definition: HoughPoint.h:76
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
Definition: fwd.h:28
float const & getVpro() const
Definition: HoughPoint.h:74
static const int fcGrid_
Maximum Number of hits in a slice that will be processes.