Merge2DTracks_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: Merge2DTracks
3 // Module Type: filter
4 // File: Merge2DTracks_module.cc
5 //
6 // Generated at Tue Jun 25 06:01:01 2013 by Matthew Tamsett using artmod
7 // from cetpkgsupport v1_02_00.
8 ////////////////////////////////////////////////////////////////////////
9 
18 #include "fhiclcpp/ParameterSet.h"
20 
25 
26 #include "DAQChannelMap/DAQChannelMap.h"
28 
29 #include <unordered_set>
30 #include <memory>
31 #include <vector>
32 //------------------------------------------------------------------------
33 namespace novaddt {
34  class Merge2DTracks;
35 
36  const float DELTA = 1e-5;
37  const float DELTA_TDC = 1e9;
38 
39  inline bool same_track(const Track *lhs, const Track *rhs) {
40  return lhs == rhs;
41  /*(fabs(lhs.StartT()-rhs.StartT())<DELTA_TDC &&
42  fabs(lhs.EndT()-rhs.EndT()) <DELTA_TDC &&
43  fabs(lhs.StartV()-rhs.StartV())<DELTA &&
44  fabs(lhs.StartZ()-rhs.StartZ())<DELTA &&
45  fabs(lhs.EndV()-rhs.EndV()) <DELTA &&
46  fabs(lhs.EndZ()-rhs.EndZ()) <DELTA &&
47  lhs.View() == rhs.View());*/
48  }
49 }
50 //------------------------------------------------------------------------
52  public:
53  explicit Merge2DTracks(fhicl::ParameterSet const & p);
54  virtual ~Merge2DTracks();
55  bool filter(art::Event & e) override;
56  void endJob() override;
57  void fillDummyTracks(std::vector< art::Ptr<novaddt::Track> > tracks, std::unique_ptr< std::vector<novaddt::Track3D> > & merged_tracks);
58  void simpleMerge(std::vector< art::Ptr<novaddt::Track> > tracks,
59  art::Handle< std::vector<novaddt::Track> > all_tracks,
60  std::unique_ptr< std::vector<novaddt::Track3D> > & merged_tracks,
61  std::unique_ptr< std::vector<novaddt::HitList> > & hitListcol,
63  void twoViewMerge(std::vector< art::Ptr<novaddt::Track> > tracks,
64  art::Handle< std::vector<novaddt::Track> > all_tracks,
65  std::unique_ptr< std::vector<novaddt::Track3D> > & merged_tracks,
66  std::unique_ptr< std::vector<novaddt::HitList> > & hitListcol,
69  double nHitsOfTrackPair(art::FindOneP<novaddt::HitList> fohl, unsigned int i, unsigned int j);
70  private:
71  std::string _trackModuleLabel; ///< label of module making the Tracks
72  std::string _trackInstanceLabel; ///< instance label making the Tracks
73  std::string _sliceModuleLabel; ///< label of module making the HitList
74  std::string _sliceInstanceLabel; ///< instance label making the HitList
75  std::string _trackToSliceInstanceLabel;///< instance of track to slice association produced by track module
76  bool _simpleTwoViewMode; ///< A boolean switch to a mode where we expect exactly two tracks and we combine these without checking metrics
77  bool _alwaysMergeMin; ///< A boolean switch to force the merge of the track pair with the minimum metric
78  std::string _matchingMetric; ///< the method used to match tracks
79  ///< currently configured options are:
80  ///< dZ - sqrt(dZ_start^2 + dZ_end^2)
81  std::string _ambiguityMetric; ///< the method used to break matching degeneracies
82  ///< currently configured options are:
83  ///< length - the longest pair are selected
84  ///< nhits - the pair with the most hits are selected
85  double _metricThreshold; ///< the maximum value of metric which can be considered matched
86  unsigned int _minhits; ///< the maximum value of metric which can be considered matched
87  int _nSlices = 0; ///< number of slices seen
88  int _nTracks = 0; ///< number of tracks seen
89  int _nMergedTracks_3D = 0; ///< number of merged tracks produced that are 3D
90  int _nMergedTracks_2D = 0; ///< number of merged tracks produced that are 2D
91  int _nMergedTracks = 0; ///< number of merged tracks produced
92 
93  std::unordered_set<size_t> unseenTracks;
94 };
95 //------------------------------------------------------------------------
97 : _trackModuleLabel (p.get<std::string >("TrackModuleLabel" ))
98  , _trackInstanceLabel (p.get<std::string >("TrackInstanceLabel", ""))
99  , _sliceModuleLabel (p.get<std::string >("SliceModuleLabel" ))
100  , _sliceInstanceLabel (p.get<std::string >("SliceInstanceLabel", ""))
101  , _trackToSliceInstanceLabel (p.get<std::string >("TrackToSliceInstanceLabel", ""))
102  , _simpleTwoViewMode (p.get<bool> ("SimpleTwoViewMode" ))
103  , _alwaysMergeMin (p.get<bool> ("AlwaysMergeMin" ))
104  , _matchingMetric (p.get<std::string >("MatchingMetric", ""))
105  , _ambiguityMetric (p.get<std::string >("AmbiguityMetric", ""))
106  , _metricThreshold (p.get<double> ("MetricThreshold" ))
107  , _minhits (p.get<int> ("MinHits" ))
108 {
109  std::cout << "=== novaddt::Merge2DTracks instantiate" << std::endl;
110  std::cout << "\t TrackModuleLabel: " << _trackModuleLabel << std::endl;
111  std::cout << "\t TrackInstanceLabel: " << _trackInstanceLabel << std::endl;
112  std::cout << "\t SliceModuleLabel: " << _sliceModuleLabel << std::endl;
113  std::cout << "\t SliceInstanceLabel: " << _sliceInstanceLabel << std::endl;
114  std::cout << "\t TrackToSliceInstanceLabel: " << _trackToSliceInstanceLabel << std::endl;
115  std::cout << "\t SimpleTwoViewMode: " << _simpleTwoViewMode << std::endl;
116  std::cout << "\t AlwaysMergeMin: " << _alwaysMergeMin << std::endl;
117  std::cout << "\t MatchingMetric: " << _matchingMetric << std::endl;
118  std::cout << "\t AmbiguityMetric: " << _ambiguityMetric << std::endl;
119  std::cout << "\t MetricThreshold: " << _metricThreshold << std::endl;
120 
121  // check that the requested matching metric is configured
122  assert( (_matchingMetric == "dZ") );
123  // check that the requested ambiguity metric is configured
124  assert( (_ambiguityMetric == "length") ||
125  (_ambiguityMetric == "nhits")
126  );
127 
128  produces< std::vector<novaddt::Track3D> >();
129  produces< std::vector<novaddt::HitList> >();
130 
131  produces< art::Assns<novaddt::Track3D, novaddt::HitList> >();
132  produces< art::Assns<novaddt::Track3D, novaddt::HitList> >(_trackToSliceInstanceLabel);
133 }
134 //------------------------------------------------------------------------
136 {
137  // Clean up dynamic memory and other resources here.
138 }
139 //------------------------------------------------------------------------
141 {
142  // Implementation of required member function here.
143  //LOG_DEBUG("Merge2DTracks") << "=== novaddt::Merge2DTracks filter. Event: "
144  // << event.id().event()
145  // << std::endl;
146 
147  // make the final data products
148  std::unique_ptr< std::vector<novaddt::Track3D> > merged_tracks (new std::vector<novaddt::Track3D>);
149  std::unique_ptr< std::vector<novaddt::HitList> > merged_hit_lists (new std::vector<novaddt::HitList>);
150  std::unique_ptr< art::Assns<novaddt::Track3D, novaddt::HitList> > assn(new art::Assns<novaddt::Track3D, novaddt::HitList>);
151  std::unique_ptr< art::Assns<novaddt::Track3D, novaddt::HitList> > assn_b(new art::Assns<novaddt::Track3D, novaddt::HitList>);
152 
153  // get the slices
155  event.getByLabel(_sliceModuleLabel, _sliceInstanceLabel, slices);
156  //LOG_DEBUG("Merge2DTracks") << "\tgot " << slices->size() << " slices" << std::endl;
157  _nSlices+=slices->size();
158 
159  //Now make it an art::PtrVector so we can use associations later
161  for(unsigned int i = 0; i < slices->size(); ++i){
162  art::Ptr<novaddt::HitList>slice(slices,i);
163  slicelist.push_back(slice);
164  }
165 
166  // ---- start of bug fix ----
167  // Due to a bug in ART v1_08_09 in order to use the one to many association later
168  // we first have to access the one to one.
169  // This finder can be used later so it's not all a waste of time
170  // in fact this might speed up the algorithm!
172  event.getByLabel(_trackModuleLabel, _trackInstanceLabel, all_tracks);
173  //// get the hit list for each track
174  art::FindOneP<novaddt::HitList> fohl(all_tracks, event, _trackModuleLabel);
175  //assert(fohl.isValid());
176  // ---- end of bug fix ----
177 
178 
179 
180  // populate unseenTracks set:
181  unseenTracks.clear();
182  for (size_t i=0; i<all_tracks->size(); ++i) unseenTracks.insert(i);
183 
184 
186  // get the tracks for each hit list
187  //art::FindManyP<novaddt::Track> find_tracks(slices, event, _trackModuleLabel);
188  art::FindManyP<novaddt::Track> find_tracks(slices, event, the_slices);
189  assert(find_tracks.isValid());
190 
191  // loop over slices
192  for(size_t i_slice = 0; i_slice < slices->size(); ++i_slice){
193  // get the hit list
194  novaddt::HitList hit_list = slices->at(i_slice);
195  if(hit_list.size() < _minhits) continue;
196  //LOG_DEBUG("Merge2DTracks") << "\tslice[" << i_slice << "]: " << hit_list.size() << " hits" << std::endl;
197  // loop over hits, useful for debugging only at the moment
198  //for(size_t i_hit = 0; i_hit < hit_list.size(); ++i_hit){
199  //novaddt::DAQHit hit(hit_list[i_hit]);
200  //std::string view = "y";
201  //if( hit.View().val == daqchannelmap::X_VIEW ) view = "x";
202  //LOG_DEBUG("Merge2DTracks") << "\t\thit[" << i_hit
203  //<< "]: TDC: " << hit.TDC().val
204  //<< ", ADC: " << hit.ADC().val
205  //<< ", plane: " << hit.Plane().val
206  //<< ", " << view
207  //<< "-cell: " << hit.Cell().val
208  //<< std::endl;
209  //} // end of loop on hits
210  // find the tracks for this hit list
211  std::vector< art::Ptr<novaddt::Track> > tracks = find_tracks.at(i_slice);
212 
213 
214  //LOG_DEBUG("Merge2DTracks") << "\t - got " << tracks.size() << " tracks" << std::endl;
215  // if there are no tracks we continue
216  if (!tracks.size()) continue;
217  _nTracks+=tracks.size();
218 
219  // if simple two view mode then we must get exactly two tracks
220  assert(!_simpleTwoViewMode || (tracks.size()==2));
221 
222  // get the hit list for each track
223  // removed in order to improve speed, see Merge functions
224  //art::FindOneP<novaddt::HitList> fohl(tracks, event, _trackModuleLabel);
225  //assert(fohl.isValid());
226  //assert(fohl.size() == tracks.size());
227 
228  //for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
229  //art::Ptr<novaddt::Track> track = tracks[i_track];
230  //std::string view = "y";
231  //if( track->View() == daqchannelmap::X_VIEW ) view = "x";
232  //// get the hit list for this track
233  //art::Ptr<novaddt::HitList> this_hit_list = fohl.at(i_track);
234 
235  ////LOG_DEBUG("Merge2DTracks") << "\t\ttrack[" << i_track
236  ////<< "]\t view: " << view
237  ////<< ", hits: " << this_hit_list->size()
238  //////<< ", time: " << track->StartT() << " - " << track->EndT()
239  ////<< ", V: " << track->StartV() << " - " << track->EndV()
240  ////<< ", Z: " << track->StartZ() << " - " << track->EndZ()
241  ////<< std::endl;
242 
243  ////// loop over hits, useful for debugging only at the moment
244  ////for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
245  ////std::string h_view = "y";
246  ////if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW) h_view = "x";
247  ////LOG_DEBUG("Merge2DTracks") << "\t\t\thit[" << i_hit
248  ////<< "]: \tplane: " << this_hit_list->at(i_hit).Plane().val
249  ////<< ",\tcell: " << this_hit_list->at(i_hit).Cell().val
250  ////<< ",\t view: " << h_view
251  ////<< ",\tTDC: " << this_hit_list->at(i_hit).TDC().val
252  ////<< ",\tADC: " << this_hit_list->at(i_hit).ADC().val
253  ////<< std::endl;
254  ////} // end of loop on hits
255  //} // end of loop on tracks
256 
257  // merge tracks
258  //fillDummyTracks(tracks, merged_tracks);
259  unsigned int current_size = merged_tracks->size();
260  if (_simpleTwoViewMode){
261  twoViewMerge(tracks, all_tracks, merged_tracks, merged_hit_lists, fohl);
262  } else {
263  simpleMerge(tracks, all_tracks, merged_tracks, merged_hit_lists, fohl);
264  }
265 
266  //LOG_DEBUG("Merge2DTracks") << "\t - produced " << merged_tracks->size() << " tracks" << std::endl;
267  //LOG_DEBUG("Merge2DTracks") << "\t - produced " << merged_hit_lists->size() << " hit lists" << std::endl;
268 
269  // the number of merged tracks should ALWAYS be equal to the number of hit lists
270  assert(merged_tracks->size()==merged_hit_lists->size());
271 
272  //LOG_DEBUG("Merge2DTracks") << "\t - saving associations" << std::endl;
273  // loop over new tracks and hit lists, save final products and associations
274  for (unsigned int i = current_size; i < merged_tracks->size(); ++i){
275  //LOG_DEBUG("Merge2DTracks") << "\t\tassociating track[" << i << "] to slice[" << i_slice << "]" << std::endl;
276  // association between this slice and each track
277  util::CreateAssn(*this, event, *merged_tracks, slicelist[i_slice], *assn_b, i);
278  // association between each track and it's hits
279  //LOG_DEBUG("Merge2DTracks") << "\t\tassociating track[" << i << "] to hit list[" << i << " - " << i+1 << "]" << std::endl;
280  util::CreateAssn(*this, event, *merged_tracks, *merged_hit_lists, *assn, i, i+1, i);
281  }
282 
283  } // end of loop on slices
284 
285  // store the merged tracks
286  _nMergedTracks += merged_tracks->size();
287  event.put(std::move(merged_tracks));
288  event.put(std::move(merged_hit_lists));
289  event.put(std::move(assn));
290  event.put(std::move(assn_b), _trackToSliceInstanceLabel);
291  return true;
292 }
293 //------------------------------------------------------------------------
295  std::unique_ptr< std::vector<novaddt::Track3D> > & merged_tracks)
296 {
297  //LOG_DEBUG("Merge2DTracks") << "\tFilling dummy tracks" << std::endl;
298  for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
299  art::Ptr<novaddt::Track> track = tracks[i_track];
300  int view = 0;
301  TVector3 start;
302  TVector3 end;
303 
304  if( track->View() == daqchannelmap::X_VIEW ){
305  view = 1;
306  start.SetXYZ(track->StartV(), 0., track->StartZ());
307  end.SetXYZ( track->EndV(), 0., track->EndZ());
308  } else {
309  start.SetXYZ(0., track->StartV(), track->StartZ());
310  end.SetXYZ( 0., track->EndV(), track->EndZ());
311  }
312 
313  novaddt::Track3D merged_track = Track3D(false, view, start, end);
314 
315  //LOG_DEBUG("Merge2DTracks") << "\t\tmerged track[" << i_track
316  //<< "]\t Is3D: " << merged_track.Is3D()
317  //<< ", view: " << merged_track.View()
318  //<< ", X: " << merged_track.Start().X() << " - " << merged_track.End().X()
319  //<< ", Y: " << merged_track.Start().Y() << " - " << merged_track.End().Y()
320  //<< ", Z: " << merged_track.Start().Z() << " - " << merged_track.End().Z()
321  //<< std::endl;
322 
323  merged_tracks->push_back(merged_track);
324  } // end of loop on tracks
325 }
326 //------------------------------------------------------------------------
328  art::Handle< std::vector<novaddt::Track> > all_tracks,
329  std::unique_ptr< std::vector<novaddt::Track3D> > & merged_tracks,
330  std::unique_ptr< std::vector<novaddt::HitList> > & merged_hit_lists,
332 {
333  //LOG_DEBUG("Merge2DTracks") << "\tMerging tracks" << std::endl;
334  int nXTracks = 0;
335  int nYTracks = 0;
336  for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
337  if (tracks[i_track]->View() == daqchannelmap::X_VIEW){
338  nXTracks++;
339  } else {
340  nYTracks++;
341  }
342  }
343 
344  //LOG_DEBUG("Merge2DTracks") << "\t - " << nXTracks << " x-tracks" << std::endl;
345  //LOG_DEBUG("Merge2DTracks") << "\t - " << nYTracks << " y-tracks" << std::endl;
346 
347  if ( (nXTracks == 0) && (nYTracks == 0) ) return;
348 
349  //LOG_DEBUG("Merge2DTracks") << "\t - make metric matrix" << std::endl;
350  // construct metric matrix
351  double dummy_metric = 1e9;
352  double metric_matrix[tracks.size()][tracks.size()];
353  // fill with large dummy values
354  for (unsigned int i = 0; i < tracks.size(); ++i){
355  for (unsigned int j = 0; j < tracks.size(); ++j){
356  metric_matrix[i][j] = dummy_metric;
357  }
358  }
359 
360  // for each track in slice, find the index into fohl that we must use
361  std::vector<size_t> fohl_index;
362 
363  for (size_t i_track = 0; i_track < tracks.size(); ++i_track) {
364  for (auto j_it = unseenTracks.begin(); j_it != unseenTracks.end(); ++j_it)
365  {
366  // this shouldn't happen, but here for safety
367  if (all_tracks->size() <= *j_it) continue;
368 
369  if (same_track(&(*tracks.at(i_track)), &all_tracks->at(*j_it))) {
370  fohl_index.push_back(*j_it);//assign(i_track, j_track);
371  // erase track at position j_track
372  unseenTracks.erase(j_it);
373  break;
374  }
375  }
376  }
377 
378  if (fohl_index.size() != tracks.size()) return; // can't match all tracks, bail
379 
380 
381  // double loop over tracks to fill metric matrix
382  //LOG_DEBUG("Merge2DTracks") << "\t - fill metric matrix" << std::endl;
383  for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
384  //LOG_DEBUG("Merge2DTracks") << "\t - x-track[" << i_track << "]" << std::endl;
385  if (tracks[i_track]->View() != daqchannelmap::X_VIEW) continue;
386  art::Ptr<novaddt::Track> x_track = tracks[i_track];
387 
388  for(size_t j_track = 0; j_track < tracks.size(); ++j_track){
389  //LOG_DEBUG("Merge2DTracks") << "\t\t - y-track[" << j_track << "]" << std::endl;
390  if (tracks[j_track]->View() == daqchannelmap::X_VIEW) continue;
391  art::Ptr<novaddt::Track> y_track = tracks[j_track];
392 
393  // define the distance measure
394  // this could be anything we want
395  // start with sqrt(dZ_low^2 + dZ_high^2)
396  double x_low_z = std::min(x_track->StartZ(),x_track->EndZ());
397  double x_high_z = std::max(x_track->StartZ(),x_track->EndZ());
398  double y_low_z = std::min(y_track->StartZ(),y_track->EndZ());
399  double y_high_z = std::max(y_track->StartZ(),y_track->EndZ());
400 
401  double dZ_low = x_low_z - y_low_z;
402  double dZ_high = x_high_z - y_high_z;
403 
404  if (_matchingMetric=="dZ") metric_matrix[i_track][j_track] = std::sqrt( (dZ_low*dZ_low) + (dZ_high*dZ_high) );
405 
406  //LOG_DEBUG("Merge2DTracks") << "\t\tmetric: " << _matchingMetric
407  //<< " between x-track[" << i_track
408  //<< "] and y-track[" << j_track
409  //<< "]: " << metric_matrix[i_track][j_track]
410  //<< std::endl;
411  } // end of loop on y
412  } // end of loop on x
413 
414  // now create merged tracks for those for which the metric is below threshold
415  int iteration_counter = 0;
416  int n_merged_tracks = 0;
417  std::vector<unsigned int> used_tracks;
418 
419  while (true){
420  iteration_counter++;
421  //LOG_DEBUG("Merge2DTracks") << "iteration: " << iteration_counter << std::endl;
422  double min_metric = 9e9;
423  int min_x = 1e9;
424  size_t min_x_index = 1e9;
425  int min_y = 1e9;
426  size_t min_y_index = 1e9;
427  double ambiguity_metric = 1e9;
428  bool any_min = false;
429 
430  // loop over x
431  for (unsigned int i = 0; i < tracks.size(); ++i){
432  if (tracks[i]->View() != daqchannelmap::X_VIEW) continue;
433  // if used continue
434  if(std::find(used_tracks.begin(), used_tracks.end(), i)!=used_tracks.end()){
435  //LOG_DEBUG("Merge2DTracks") << "\tx-track[" << i << "]: used" << std::endl;
436  continue;
437  }
438 
439  size_t i_index = fohl_index.at(i);
440  // loop over y
441  for (unsigned int j = 0; j < tracks.size(); ++j){
442  if (tracks[j]->View() == daqchannelmap::X_VIEW) continue;
443  // if used continue
444  if(std::find(used_tracks.begin(), used_tracks.end(), j)!=used_tracks.end()){
445  //LOG_DEBUG("Merge2DTracks") << "\ty-track[" << j << "]: used" << std::endl;
446  continue;
447  }
448  // get metric
449  //LOG_DEBUG("Merge2DTracks") << "\t - metric between x-track[" << i
450  //<< "] and y-track[" << j
451  //<< "]: " << metric_matrix[i][j] <<std::endl;
452  // if the value wasn't filled continue
453  if (metric_matrix[i][j] == dummy_metric) continue;
454 
455  size_t j_index = fohl_index.at(j);
456 
457  // check if new minimum
458  bool new_min = false;
459  if (metric_matrix[i][j] < min_metric) new_min = true;
460  // if identical pick the pair with the lowest ambiguity metric
461  if (metric_matrix[i][j] == min_metric){
462  double this_ambiguity_metric = 0;
463  if (_ambiguityMetric=="length") this_ambiguity_metric = lengthOfTrackPair(tracks[i],tracks[j]);
464  if (_ambiguityMetric=="nhits") this_ambiguity_metric = nHitsOfTrackPair(fohl,i_index,j_index);
465  //LOG_DEBUG("Merge2DTracks") << "\t\t - indentical metrics, picking based on " << _ambiguityMetric << ": " << this_ambiguity_metric << std::endl;
466  if (this_ambiguity_metric > ambiguity_metric) new_min = true;
467  } // end of ambiguity solver
468 
469  if (new_min){
470  // save new min
471  min_metric = metric_matrix[i][j];
472  min_x = i;
473  min_x_index = i_index;
474  min_y = j;
475  min_y_index = j_index;
476  any_min = true;
477  if (_ambiguityMetric=="length") ambiguity_metric = lengthOfTrackPair(tracks[i],tracks[j]);
478  if (_ambiguityMetric=="nhits") ambiguity_metric = nHitsOfTrackPair(fohl,i_index,j_index);
479 
480  //LOG_DEBUG("Merge2DTracks") << "\t\t - new min, " << _ambiguityMetric
481  //<< ": " << ambiguity_metric << std::endl;
482  }
483  } // end of loop on y
484  } // end of loop on x
485 
486  if ((min_metric<_metricThreshold) || \
487  (_alwaysMergeMin && (n_merged_tracks == 0) && any_min)){
488  n_merged_tracks++;
489  //LOG_DEBUG("Merge2DTracks") << "Saving minimum: " << min_metric
490  //<< " from [" << min_x
491  //<< "," << min_y
492  //<< "]" << std::endl;
493  used_tracks.push_back(min_x);
494  used_tracks.push_back(min_y);
495 
496  // save merged track
497  TVector3 start;
498  TVector3 end;
499  // assume north going, that is the min Z is the start point
500  // find the minimum x-z
501  // assume that track is already north going
502  double start_x = tracks[min_x]->StartV();
503  double end_x = tracks[min_x]->EndV();
504  double min_xz_z = tracks[min_x]->StartZ();
505  double max_xz_z = tracks[min_x]->EndZ();
506  // but check for south going
507  if (tracks[min_x]->StartZ() > tracks[min_x]->EndZ() ){
508  start_x = tracks[min_x]->EndV();
509  end_x = tracks[min_x]->StartV();
510  min_xz_z = tracks[min_x]->EndZ();
511  max_xz_z = tracks[min_x]->StartZ();
512  }
513  // repeat for y-x
514  double start_y = tracks[min_y]->StartV();
515  double end_y = tracks[min_y]->EndV();
516  double min_yz_z = tracks[min_y]->StartZ();
517  double max_yz_z = tracks[min_y]->EndZ();
518  // but check for south going
519  if (tracks[min_y]->StartZ() > tracks[min_y]->EndZ() ){
520  start_y = tracks[min_y]->EndV();
521  end_y = tracks[min_y]->StartV();
522  min_yz_z = tracks[min_y]->EndZ();
523  max_yz_z = tracks[min_y]->StartZ();
524  }
525 
526  double start_z = std::min(min_xz_z, min_yz_z);
527  double end_z = std::max(max_xz_z, max_yz_z);
528 
529  start.SetXYZ(start_x,
530  start_y,
531  start_z);
532  end.SetXYZ( end_x,
533  end_y,
534  end_z);
535 
536  novaddt::Track3D merged_track = Track3D(true, 2, start, end);
537 
538  //LOG_DEBUG("Merge2DTracks") << "\t\tmerged track:"
539  //<< ", Is3D: " << merged_track.Is3D()
540  //<< ", view: " << merged_track.View()
541  //<< ", X: " << merged_track.Start().X() << " - " << merged_track.End().X()
542  //<< ", Y: " << merged_track.Start().Y() << " - " << merged_track.End().Y()
543  //<< ", Z: " << merged_track.Start().Z() << " - " << merged_track.End().Z()
544  //<< std::endl;
545  merged_tracks->push_back(merged_track);
546 
547  // save merged hit list
548  novaddt::HitList merged_hit_list;
549  // get the hit list for each track
550  novaddt::HitList x_hit_list = *fohl.at(min_x_index);
551  for (auto const& hit : x_hit_list ) {
552  merged_hit_list.emplace_back(hit);
553  }
554  novaddt::HitList y_hit_list = *fohl.at(min_y_index);
555  for (auto const& hit : y_hit_list ) {
556  merged_hit_list.emplace_back(hit);
557  }
558  merged_hit_lists->push_back(merged_hit_list);
559 
561  } else {
562  //LOG_DEBUG("Merge2DTracks") << "Minimum: " << min_metric
563  //<< " > threshold: " << _metricThreshold
564  //<< std::endl;
565  break;
566  }
567  } // end of while
568  // loop over tracks and store the remaining ones as 2-D tracks
569  for(unsigned int i_track = 0; i_track < tracks.size(); ++i_track){
570  if(std::find(used_tracks.begin(), used_tracks.end(), i_track)!=used_tracks.end()) continue;
571  art::Ptr<novaddt::Track> track = tracks[i_track];
572  int view = 0;
573  TVector3 start;
574  TVector3 end;
575 
576  if( track->View() == daqchannelmap::X_VIEW ){
577  view = 1;
578  start.SetXYZ(track->StartV(), 0., track->StartZ());
579  end.SetXYZ( track->EndV(), 0., track->EndZ());
580  } else {
581  start.SetXYZ(0., track->StartV(), track->StartZ());
582  end.SetXYZ( 0., track->EndV(), track->EndZ());
583  }
584 
585  novaddt::Track3D merged_track = Track3D(false, view, start, end);
586 
587  //LOG_DEBUG("Merge2DTracks") << "\t\tmerged (2D) track[" << i_track
588  //<< "]\t Is3D: " << merged_track.Is3D()
589  //<< ", view: " << merged_track.View()
590  //<< ", X: " << merged_track.Start().X() << " - " << merged_track.End().X()
591  //<< ", Y: " << merged_track.Start().Y() << " - " << merged_track.End().Y()
592  //<< ", Z: " << merged_track.Start().Z() << " - " << merged_track.End().Z()
593  //<< std::endl;
594 
595  merged_tracks->push_back(merged_track);
596 
597  // save merged hit list
598  novaddt::HitList merged_hit_list;
599  // get the hit list for the track
600  novaddt::HitList hit_list = *fohl.at(fohl_index[i_track]);
601  for (auto const& hit : hit_list ) {
602  merged_hit_list.emplace_back(hit);
603  }
604  merged_hit_lists->push_back(merged_hit_list);
605 
607  } // end of loop on tracks
608  //assert(false);
609 }
610 //------------------------------------------------------------------------
612  art::Handle< std::vector<novaddt::Track> > all_tracks,
613  std::unique_ptr< std::vector<novaddt::Track3D> > & merged_tracks,
614  std::unique_ptr< std::vector<novaddt::HitList> > & merged_hit_lists,
616 {
617  //LOG_DEBUG("Merge2DTracks") << "\tMerging tracks" << std::endl;
618  // this method assumes that the Tracking module or an equivalent has made these tracks
619  // such tracks always come in pairs and the x track is always the first
620  // the fact that there is two has been asserted previously.
621  assert(tracks[0]->View() == daqchannelmap::X_VIEW);
622  assert(tracks[1]->View() != daqchannelmap::X_VIEW);
623 
624  // save merged track
625  TVector3 start;
626  TVector3 end;
627  // assume north going, that is the min Z is the start point
628  // find the minimum x-z
629  // assume that track is already north going
630  double start_x = tracks[0]->StartV();
631  double end_x = tracks[0]->EndV();
632  double min_xz_z = tracks[0]->StartZ();
633  double max_xz_z = tracks[0]->EndZ();
634  // but check for south going
635  if (tracks[0]->StartZ() > tracks[0]->EndZ() ){
636  start_x = tracks[0]->EndV();
637  end_x = tracks[0]->StartV();
638  min_xz_z = tracks[0]->EndZ();
639  max_xz_z = tracks[0]->StartZ();
640  }
641  // repeat for y-x
642  double start_y = tracks[1]->StartV();
643  double end_y = tracks[1]->EndV();
644  double min_yz_z = tracks[1]->StartZ();
645  double max_yz_z = tracks[1]->EndZ();
646  // but check for south going
647  if (tracks[1]->StartZ() > tracks[1]->EndZ() ){
648  start_y = tracks[1]->EndV();
649  end_y = tracks[1]->StartV();
650  min_yz_z = tracks[1]->EndZ();
651  max_yz_z = tracks[1]->StartZ();
652  }
653 
654  size_t zero_index, one_index;
655  bool zero_found = false, one_found = false;
656  for (size_t j=0; j<all_tracks->size(); ++j)
657  {
658  if (!zero_found && same_track(&(*tracks[0]),&all_tracks->at(j))) {
659  zero_index = j;
660  zero_found = true;
661  }
662  if (!one_found && same_track(&(*tracks[1]),&all_tracks->at(j))) {
663  one_index = j;
664  one_found = true;
665  }
666  if (zero_found && one_found) break;
667  }
668 
669  novaddt::HitList x_hit_list = *fohl.at(zero_index);
670  novaddt::HitList y_hit_list = *fohl.at(one_index);
671  double start_z,end_z;
672  // If the x hit list is empty set z from to y
673  if (!x_hit_list.size()){
674  start_z = min_yz_z;
675  end_z = max_yz_z;
676  // else if the y hit list is empty set z from to x
677  } else if (!y_hit_list.size()){
678  start_z = min_xz_z;
679  end_z = max_xz_z;
680  // else use the minimum and maximum of both
681  } else {
682  start_z = std::min(min_xz_z, min_yz_z);
683  end_z = std::max(max_xz_z, max_yz_z);
684  }
685 
686  start.SetXYZ(start_x,
687  start_y,
688  start_z);
689  end.SetXYZ( end_x,
690  end_y,
691  end_z);
692 
693  novaddt::Track3D merged_track = Track3D(true, 2, start, end);
694 
695  //LOG_DEBUG("Merge2DTracks") << "\t\tmerged track:"
696  //<< ", Is3D: " << merged_track.Is3D()
697  //<< ", view: " << merged_track.View()
698  //<< ", X: " << merged_track.Start().X() << " - " << merged_track.End().X()
699  //<< ", Y: " << merged_track.Start().Y() << " - " << merged_track.End().Y()
700  //<< ", Z: " << merged_track.Start().Z() << " - " << merged_track.End().Z()
701  //<< std::endl;
702  merged_tracks->push_back(merged_track);
703 
704  // save merged hit list
705  novaddt::HitList merged_hit_list;
706  // get the hit list for each track
707 
708  for (auto const& hit : x_hit_list ) {
709  merged_hit_list.emplace_back(hit);
710  }
711  for (auto const& hit : y_hit_list ) {
712  merged_hit_list.emplace_back(hit);
713  }
714  merged_hit_lists->push_back(merged_hit_list);
715 
717 }
718 //------------------------------------------------------------------------
720 {
721  double this_dX = std::fabs(x_track->EndV() - x_track->StartV());
722  double this_dXZ = std::fabs(x_track->EndZ() - x_track->StartZ());
723  double this_dY = std::fabs(y_track->EndV() - y_track->StartV());
724  double this_dYZ = std::fabs(y_track->EndZ() - y_track->StartZ());
725  double length = std::sqrt( (this_dX*this_dX) +
726  (this_dXZ*this_dXZ) +
727  (this_dY*this_dY) +
728  (this_dYZ*this_dYZ) );
729  return length;
730 }
731 //------------------------------------------------------------------------
733 {
734  double n_hits = 0;
735 
736  // get the hit list for each track
737  art::Ptr<novaddt::HitList> x_hit_list = fohl.at(i);
738  art::Ptr<novaddt::HitList> y_hit_list = fohl.at(j);
739 
740  n_hits+=double(x_hit_list->size());
741  n_hits+=double(y_hit_list->size());
742 
743  return n_hits;
744 }
745 //------------------------------------------------------------------------
747 {
748  std::cout << "--- novaddt::Merge2DTracks endJob " << std::endl;
749  std::cout << "\tNumber of slices seen: " << _nSlices << std::endl;
750  std::cout << "\tNumber of tracks seen: " << _nTracks << std::endl;
751  std::cout << "\tNumber of merged tracks made: " << _nMergedTracks << std::endl;
752  std::cout << "\tNumber of merged 3D tracks: " << _nMergedTracks_3D << std::endl;
753  std::cout << "\tNumber of merged 2D tracks: " << _nMergedTracks_2D << std::endl;
754 }
755 //------------------------------------------------------------------------
const float DELTA
T max(const caf::Proxy< T > &a, T b)
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double _metricThreshold
the maximum value of metric which can be considered matched
float const & StartZ() const
Definition: Track.h:68
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Definition: event.h:19
std::string _sliceModuleLabel
label of module making the HitList
uint8_t const & View() const
Definition: Track.h:65
int _nMergedTracks
number of merged tracks produced
DEFINE_ART_MODULE(TestTMapFile)
bool _alwaysMergeMin
A boolean switch to force the merge of the track pair with the minimum metric.
std::string _trackToSliceInstanceLabel
instance of track to slice association produced by track module
int _nMergedTracks_3D
number of merged tracks produced that are 3D
bool same_track(const Track *lhs, const Track *rhs)
double nHitsOfTrackPair(art::FindOneP< novaddt::HitList > fohl, unsigned int i, unsigned int j)
float const & EndV() const
Definition: Track.h:70
length
Definition: demo0.py:21
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
int _nTracks
number of tracks seen
bool filter(art::Event &e) override
void twoViewMerge(std::vector< art::Ptr< novaddt::Track > > tracks, art::Handle< std::vector< novaddt::Track > > all_tracks, std::unique_ptr< std::vector< novaddt::Track3D > > &merged_tracks, std::unique_ptr< std::vector< novaddt::HitList > > &hitListcol, art::FindOneP< novaddt::HitList > fohl)
Identifier for the X measuring view of the detector (top)
const double j
Definition: BetheBloch.cxx:29
int _nSlices
number of slices seen
Definition: View.py:1
const float DELTA_TDC
OStream cout
Definition: OStream.cxx:6
std::string _trackInstanceLabel
instance label making the Tracks
std::string _trackModuleLabel
label of module making the Tracks
std::unordered_set< size_t > unseenTracks
bool _simpleTwoViewMode
A boolean switch to a mode where we expect exactly two tracks and we combine these without checking m...
Definition: structs.h:12
double lengthOfTrackPair(art::Ptr< novaddt::Track > x_track, art::Ptr< novaddt::Track > y_track)
std::string _sliceInstanceLabel
instance label making the HitList
Merge2DTracks(fhicl::ParameterSet const &p)
float const & StartV() const
Definition: Track.h:67
unsigned int _minhits
the maximum value of metric which can be considered matched
assert(nhit_max >=nhit_nbins)
void simpleMerge(std::vector< art::Ptr< novaddt::Track > > tracks, art::Handle< std::vector< novaddt::Track > > all_tracks, std::unique_ptr< std::vector< novaddt::Track3D > > &merged_tracks, std::unique_ptr< std::vector< novaddt::HitList > > &hitListcol, art::FindOneP< novaddt::HitList > fohl)
int _nMergedTracks_2D
number of merged tracks produced that are 2D
float const & EndZ() const
Definition: Track.h:71
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
Definition: fwd.h:28
void fillDummyTracks(std::vector< art::Ptr< novaddt::Track > > tracks, std::unique_ptr< std::vector< novaddt::Track3D > > &merged_tracks)
enum BeamMode string