LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ClusterMatchAlg.cxx
Go to the documentation of this file.
1 //
3 // ClusterMatchAlg source
4 //
6 
7 #include "ClusterMatchAlg.h"
8 
9 #include "TString.h"
10 #include "TTree.h"
11 
14 #include "art_root_io/TFileService.h"
17 #include "fhiclcpp/ParameterSet.h"
19 
28 
29 namespace cluster {
30 
31  ClusterMatchAlg::ClusterMatchAlg(fhicl::ParameterSet const& pset) : _ModName_MCTruth("")
32  {
33  _debug_mode = pset.get<bool>("DebugMode");
34  _store_sps = pset.get<bool>("StoreSpacePoint");
35  _num_sps_cut = pset.get<size_t>("CutParam_NumSpacePoint");
36  _overlay_tratio_cut = pset.get<double>("CutParam_OverlayTimeFraction");
37  _qratio_cut = pset.get<double>("CutParam_SumChargeRatio");
38  std::vector<size_t> algo_list = pset.get<std::vector<size_t>>("MatchAlgoList");
39 
40  _sps_algo = new trkf::SpacePointAlg(pset.get<fhicl::ParameterSet>("SpacePointAlg"));
41  _match_tree = 0;
42  _cluster_tree = 0;
43  _save_cluster_info = true;
44  _det_params_prepared = false;
45 
46  for (size_t i = 0; i < (size_t)(kMATCH_METHOD_MAX); ++i)
47 
48  _match_methods[i] = false;
49 
50  for (auto const v : algo_list) {
51 
52  if (v >= (size_t)(kMATCH_METHOD_MAX))
53 
54  mf::LogError("ClusterMatchAlg") << Form("Invalid algorithm enum: %zu", v);
55 
56  else
57  _match_methods[v] = true;
58  }
59 
60  ReportConfig();
61 
63  }
64 
66  {
67  std::ostringstream msg;
68  msg << std::endl
69  << " ClusterMatchAlg Configuration: " << std::endl
70  << "---------------------------------------------" << std::endl;
71  msg << " Debug Mode ... " << (_debug_mode ? "enabled!" : "disabled!") << std::endl;
72  msg << " RoughZ ....... " << (_match_methods[kRoughZ] ? "enabled!" : "disabled!") << std::endl;
73  msg << " RoughT ....... " << (_match_methods[kRoughT] ? "enabled!" : "disabled!") << std::endl;
74  msg << " SpacePoint ... " << (_match_methods[kSpacePoint] ? "enabled!" : "disabled!")
75  << std::endl;
76  msg << " SumCharge .... " << (_match_methods[kSumCharge] ? "enabled!" : "disabled!")
77  << std::endl;
78  msg << std::endl;
79  msg << " Overlay-Time Fraction Cut : " << _overlay_tratio_cut << std::endl
80  << " Charge-Ratio Diff. Cut : " << _qratio_cut << std::endl
81  << " Minimum # of SpacePoint : " << _num_sps_cut << std::endl
82  << std::endl;
83  msg << "---------------------------------------------" << std::endl;
84 
85  mf::LogWarning("ClusterMatchAlg") << msg.str();
86  }
87 
89  {
90  _ucluster_v.clear();
91  _vcluster_v.clear();
92  _wcluster_v.clear();
93 
94  _uhits_v.clear();
95  _vhits_v.clear();
96  _whits_v.clear();
97  }
98 
100  {
101  _matched_uclusters_v.clear();
102  _matched_vclusters_v.clear();
103  _matched_wclusters_v.clear();
104  _matched_sps_v.clear();
105  }
106 
108  {
109 
110  _mc_E = 0;
111  _mc_Px = 0;
112  _mc_Py = 0;
113  _mc_Pz = 0;
114  _mc_Vx = 0;
115  _mc_Vy = 0;
116  _mc_Vz = 0;
117  _pdgid = 0;
118  _tot_u = 0;
119  _tot_v = 0;
120  _tot_w = 0;
121  _tot_pass_z = 0;
122  _tot_pass_t = 0;
123  _tot_pass_sps = 0;
124  _tot_pass_qsum = 0;
125  _qratio_v.clear();
126  _uv_tratio_v.clear();
127  _vw_tratio_v.clear();
128  _wu_tratio_v.clear();
129  _u_nhits_v.clear();
130  _v_nhits_v.clear();
131  _w_nhits_v.clear();
132  _nsps.clear();
133 
134  _view_v.clear();
135  _charge_v.clear();
136  _nhits_v.clear();
137  _tstart_min_v.clear();
138  _tstart_max_v.clear();
139  _tpeak_min_v.clear();
140  _tpeak_max_v.clear();
141  _tend_min_v.clear();
142  _tend_max_v.clear();
143  }
144 
146  {
149  ClearTTreeInfo();
150  }
151 
153  {
154  if (!_match_tree) {
156  _match_tree = fileService->make<TTree>("match_tree", "");
157  _match_tree->Branch("mc_E", &_mc_E, "mc_E/D");
158  _match_tree->Branch("mc_Px", &_mc_Px, "mc_Px/D");
159  _match_tree->Branch("mc_Py", &_mc_Py, "mc_Py/D");
160  _match_tree->Branch("mc_Pz", &_mc_Pz, "mc_Pz/D");
161  _match_tree->Branch("mc_Vx", &_mc_Vx, "mc_Vx/D");
162  _match_tree->Branch("mc_Vy", &_mc_Vy, "mc_Vy/D");
163  _match_tree->Branch("mc_Vz", &_mc_Vz, "mc_Vz/D");
164 
165  _match_tree->Branch("pdgid", &_pdgid, "pdgid/I");
166  _match_tree->Branch("tot_u", &_tot_u, "tot_u/s");
167  _match_tree->Branch("tot_v", &_tot_v, "tot_v/s");
168  _match_tree->Branch("tot_w", &_tot_w, "tot_w/s");
169  _match_tree->Branch("tot_pass_t", &_tot_pass_t, "tot_pass_t/s");
170  _match_tree->Branch("tot_pass_z", &_tot_pass_z, "tot_pass_z/s");
171  _match_tree->Branch("tot_pass_sps", &_tot_pass_sps, "tot_pass_sps/s");
172  _match_tree->Branch("tot_pass_qsum", &_tot_pass_qsum, "tot_pass_qsum/s");
173 
174  _match_tree->Branch("uv_tratio_v", "std::vector<double>", &_uv_tratio_v);
175  _match_tree->Branch("vw_tratio_v", "std::vector<double>", &_vw_tratio_v);
176  _match_tree->Branch("wu_tratio_v", "std::vector<double>", &_wu_tratio_v);
177 
178  _match_tree->Branch("qratio_v", "std::vector<double>", &_qratio_v);
179  _match_tree->Branch("u_nhits_v", "std::vector<UShort_t>", &_u_nhits_v);
180  _match_tree->Branch("v_nhits_v", "std::vector<UShort_t>", &_v_nhits_v);
181  _match_tree->Branch("w_nhits_v", "std::vector<UShort_t>", &_w_nhits_v);
182  _match_tree->Branch("nsps", "std::vector<UShort_t>", &_nsps);
183  }
186  _cluster_tree = fileService->make<TTree>("cluster_tree", "");
187  _cluster_tree->Branch("view_v", "std::vector<uint16_t>", &_view_v);
188  _cluster_tree->Branch("charge_v", "std::vector<double>", &_charge_v);
189  _cluster_tree->Branch("nhits_v", "std::vector<uint16_t>", &_nhits_v);
190  _cluster_tree->Branch("tstart_min_v", "std::vector<double>", &_tstart_min_v);
191  _cluster_tree->Branch("tstart_max_v", "std::vector<double>", &_tstart_max_v);
192  _cluster_tree->Branch("tpeak_min_v", "std::vector<double>", &_tpeak_min_v);
193  _cluster_tree->Branch("tpeak_max_v", "std::vector<double>", &_tpeak_max_v);
194  _cluster_tree->Branch("tend_min_v", "std::vector<double>", &_tend_min_v);
195  _cluster_tree->Branch("tend_max_v", "std::vector<double>", &_tend_max_v);
196  }
197  }
198 
200  {
201  if (!_ModName_MCTruth.size()) return;
202 
203  std::vector<const simb::MCTruth*> mciArray;
204 
205  try {
206 
207  evt.getView(_ModName_MCTruth, mciArray);
208  }
209  catch (art::Exception const& e) {
210 
211  if (e.categoryCode() != art::errors::ProductNotFound) throw;
212  }
213 
214  for (size_t i = 0; i < mciArray.size(); ++i) {
215 
216  if (i == 1) {
217  mf::LogWarning("ClusterMatchAlg") << " Ignoring > 2nd MCTruth in MC generator...";
218  break;
219  }
220  const simb::MCTruth* mci_ptr(mciArray.at(i));
221 
222  for (size_t j = 0; j < (size_t)(mci_ptr->NParticles()); ++j) {
223 
224  if (j == 1) {
225  mf::LogWarning("ClusterMatchAlg") << " Ignoring > 2nd MCParticle in MC generator...";
226  break;
227  }
228 
229  const simb::MCParticle part(mci_ptr->GetParticle(j));
230 
231  _pdgid = part.PdgCode();
232  _mc_E = part.E();
233  _mc_Px = part.Px();
234  _mc_Py = part.Py();
235  _mc_Pz = part.Pz();
236  _mc_Vx = part.Vx();
237  _mc_Vy = part.Vy();
238  _mc_Vz = part.Vz();
239  }
240  }
241  }
242 
244  {
245  if (!_det_params_prepared) {
246  // Total number of planes
248  _tot_planes = geo->Nplanes();
249 
250  // Ask DetectorPrperties about time-offset among different wire planes ... used to correct timing
251  // difference among different wire planes in the following loop.
252  _time_offset_uplane = det_prop.GetXTicksOffset(geo::kU, 0, 0);
253  _time_offset_vplane = det_prop.GetXTicksOffset(geo::kV, 0, 0);
255  if (_tot_planes > 2) _time_offset_wplane = det_prop.GetXTicksOffset(geo::kW, 0, 0);
256  _det_params_prepared = true;
257  }
258  }
259 
261  const recob::Cluster& in_cluster,
262  const std::vector<art::Ptr<recob::Hit>>& in_hit_v)
263  {
264  PrepareDetParams(det_prop);
265  cluster_match_info ci(in_cluster.ID());
266  ci.view = in_cluster.View();
267 
269  FillHitInfo(ci, hit_ptrv, in_hit_v);
270 
271  // Save created art::PtrVector & cluster_match_info struct object
272  switch (ci.view) {
273  case geo::kU:
274  _uhits_v.push_back(hit_ptrv);
275  _ucluster_v.push_back(ci);
277  break;
278  case geo::kV:
279  _vhits_v.push_back(hit_ptrv);
280  _vcluster_v.push_back(ci);
282  break;
283  case geo::kW:
284  _whits_v.push_back(hit_ptrv);
285  _wcluster_v.push_back(ci);
287  break;
288  default:
289  mf::LogError("ClusterMatchAlg") << Form("Found an invalid plane ID: %d", in_cluster.View());
290  }
291  }
292 
294  art::PtrVector<recob::Hit>& out_hit_v,
295  const std::vector<art::Ptr<recob::Hit>>& in_hit_v)
296  {
297 
298  out_hit_v.reserve(in_hit_v.size());
299 
300  double time_offset = 0;
301  if (ci.view == geo::kU)
302  time_offset = _time_offset_uplane;
303  else if (ci.view == geo::kV)
304  time_offset = _time_offset_vplane;
305  else if (ci.view == geo::kW)
306  time_offset = _time_offset_wplane;
307 
308  // Loop over hits in this cluster
309  for (auto const hit : in_hit_v) {
310 
311  unsigned int wire = hit->WireID().Wire;
312  double tstart = hit->PeakTimePlusRMS(-1.) - time_offset;
313  double tpeak = hit->PeakTime() - time_offset;
314  double tend = hit->PeakTimePlusRMS(+1.) - time_offset;
315 
316  ci.sum_charge += hit->Integral();
317 
318  ci.wire_max = (ci.wire_max < wire) ? wire : ci.wire_max;
319  ci.wire_min = (ci.wire_min > wire) ? wire : ci.wire_min;
320 
321  ci.start_time_max = (ci.start_time_max < tstart) ? tstart : ci.start_time_max;
322  ci.peak_time_max = (ci.peak_time_max < tpeak) ? tpeak : ci.peak_time_max;
323  ci.end_time_max = (ci.end_time_max < tend) ? tend : ci.end_time_max;
324 
325  ci.start_time_min = (ci.start_time_min > tstart) ? tstart : ci.start_time_min;
326  ci.peak_time_min = (ci.peak_time_min > tpeak) ? tpeak : ci.peak_time_min;
327  ci.end_time_min = (ci.end_time_min > tend) ? tend : ci.end_time_min;
328 
329  out_hit_v.push_back(hit);
330  }
331 
332  ci.nhits = in_hit_v.size();
333  }
334 
336  {
337  if (_cluster_tree) {
338  _view_v.push_back(ci.view);
339  _charge_v.push_back(ci.sum_charge);
340  _nhits_v.push_back(ci.nhits);
341  _tstart_min_v.push_back(ci.start_time_min);
342  _tstart_max_v.push_back(ci.start_time_max);
343  _tpeak_min_v.push_back(ci.peak_time_min);
344  _tpeak_max_v.push_back(ci.peak_time_max);
345  _tend_min_v.push_back(ci.end_time_min);
346  _tend_max_v.push_back(ci.end_time_max);
347  }
348  }
349 
350  //########################################################################################
352  const cluster_match_info& ci2,
353  const geo::View_t v1,
354  const geo::View_t v2) const
355  {
357  double y, z_min, z_max;
358  y = z_min = z_max = -1;
359  constexpr geo::TPCID tpcid{0, 0};
360  geo::PlaneID const plane_1{tpcid, v1};
361  geo::PlaneID const plane_2{tpcid, v2};
362  geo_h->IntersectionPoint(
363  geo::WireID{plane_1, ci1.wire_min}, geo::WireID{plane_2, ci2.wire_min}, y, z_min);
364  geo_h->IntersectionPoint(
365  geo::WireID{plane_1, ci1.wire_max}, geo::WireID{plane_2, ci2.wire_max}, y, z_max);
366  return z_max > z_min;
367  }
368 
370  const cluster_match_info& ci2)
371  {
372  double time_overlay = std::min(ci1.end_time_max, ci2.end_time_max) -
373  std::max(ci1.start_time_min, ci2.start_time_min);
374 
375  double overlay_tratio =
376  time_overlay /
377  (ci1.end_time_max - ci1.start_time_min + ci2.end_time_max - ci2.start_time_min) * 2.;
378 
379  if ((ci1.view == geo::kU && ci2.view == geo::kV) ||
380  (ci1.view == geo::kV && ci2.view == geo::kU))
381  _uv_tratio_v.push_back(overlay_tratio);
382  else if ((ci1.view == geo::kV && ci2.view == geo::kW) ||
383  (ci1.view == geo::kW && ci2.view == geo::kV))
384  _vw_tratio_v.push_back(overlay_tratio);
385  else if ((ci1.view == geo::kW && ci2.view == geo::kU) ||
386  (ci1.view == geo::kU && ci2.view == geo::kW))
387  _wu_tratio_v.push_back(overlay_tratio);
388 
389  return (overlay_tratio > _overlay_tratio_cut);
390  }
391 
393  {
394  double qratio = (uc.sum_charge) / (vc.sum_charge);
395 
396  // For quality check log
397  _qratio_v.push_back(qratio);
398 
399  return ((1 - _qratio_cut) < qratio && (qratio) < (1 + _qratio_cut));
400  }
401 
403  detinfo::DetectorPropertiesData const& det_prop,
404  const size_t uindex,
405  const size_t vindex,
406  const size_t windex,
407  std::vector<recob::SpacePoint>& sps_v)
408  {
409  bool use_wplane = _tot_planes > 2;
410 
411  if (uindex >= _ucluster_v.size() || vindex >= _vcluster_v.size() ||
412  (use_wplane && (windex >= _wcluster_v.size()))) {
413 
414  mf::LogError("ClusterMatchAlg")
415  << std::endl
416  << Form(
417  "Requested to cluster-index (U,V,W) = (%zu,%zu,%zu) where max-length is (%zu,%zu,%zu)",
418  uindex,
419  vindex,
420  windex,
421  _ucluster_v.size(),
422  _vcluster_v.size(),
423  _wcluster_v.size())
424  << std::endl;
425  return false;
426  }
427 
428  // Define a time range in which hits are used for spacepoint finding ... here "peak time" is the relevant one
429  double trange_min =
430  std::min(_ucluster_v.at(uindex).peak_time_min, _vcluster_v.at(vindex).peak_time_min);
431  if (use_wplane) trange_min = std::min(trange_min, _wcluster_v.at(windex).peak_time_min);
432 
433  double trange_max =
434  std::max(_ucluster_v.at(uindex).peak_time_max, _vcluster_v.at(vindex).peak_time_max);
435  if (use_wplane) trange_max = std::max(trange_max, _wcluster_v.at(windex).peak_time_max);
436 
437  // Space-point algorithm applies additional dT
438  trange_min -= _sps_algo->maxDT();
439  trange_max += _sps_algo->maxDT();
440 
441  // Make PtrVector<recob::Hit> for relevant Hits
442  art::PtrVector<recob::Hit> hit_group;
443  size_t max_size = _uhits_v.at(uindex).size() + _vhits_v.at(vindex).size();
444  if (use_wplane) max_size += _whits_v.at(windex).size();
445  hit_group.reserve(max_size);
446  // Loop over hits in U-plane
447  for (auto const hit : _uhits_v.at(uindex)) {
448  if (hit->PeakTime() < trange_min) continue;
449  if (hit->PeakTime() > trange_max) continue;
450  hit_group.push_back(hit);
451  }
452  // Check if any hit found in this plane
453  size_t u_nhits = hit_group.size();
454  if (!u_nhits && !_debug_mode) return false;
455  // Loop over hits in V-plane
456  for (auto const hit : _vhits_v.at(vindex)) {
457  if (hit->PeakTime() < trange_min) continue;
458  if (hit->PeakTime() > trange_max) continue;
459  hit_group.push_back(hit);
460  }
461  // Check if any hit found in this plane
462  size_t v_nhits = hit_group.size() - u_nhits;
463  if (!(v_nhits) && !_debug_mode) return false;
464 
465  // Loop over hits in W-plane
466  if (use_wplane) {
467  for (auto const hit : _whits_v.at(windex)) {
468  if (hit->PeakTime() < trange_min) continue;
469  if (hit->PeakTime() > trange_max) continue;
470  hit_group.push_back(hit);
471  }
472  }
473  // Check if any hit found in this plane
474  size_t w_nhits = hit_group.size() - u_nhits - v_nhits;
475  if (!(w_nhits) && use_wplane && !_debug_mode) return false;
476 
477  // Run SpacePoint finder algo
478  if (u_nhits && v_nhits && (!use_wplane || (w_nhits && use_wplane))) {
480  _sps_algo->makeSpacePoints(clock_data, det_prop, hit_group, sps_v);
481  }
482 
483  size_t nsps = sps_v.size();
484  _u_nhits_v.push_back(u_nhits);
485  _v_nhits_v.push_back(v_nhits);
486  if (use_wplane) _w_nhits_v.push_back(w_nhits);
487  _nsps.push_back(nsps);
488 
489  if (nsps < _num_sps_cut) return false;
490  return true;
491  }
492 
493  std::vector<std::vector<unsigned int>> ClusterMatchAlg::GetMatchedClusters() const
494  {
495  std::vector<std::vector<unsigned int>> result;
496  result.push_back(_matched_uclusters_v);
497  result.push_back(_matched_vclusters_v);
498  result.push_back(_matched_wclusters_v);
499  return result;
500  }
501 
503  detinfo::DetectorPropertiesData const& detProp)
504  {
505  std::ostringstream msg;
506  msg << Form("Received (U,V,W) = (%zu,%zu,%zu) clusters...",
507  _uhits_v.size(),
508  _vhits_v.size(),
509  _whits_v.size())
510  << std::endl;
511  _tot_u = _ucluster_v.size();
512  _tot_v = _vcluster_v.size();
513  _tot_w = _wcluster_v.size();
514 
515  if (!(_tot_u + _tot_v + _tot_w)) {
516 
517  mf::LogError(__PRETTY_FUNCTION__)
518  << "No input cluster info found! Aborting the function call...";
519 
520  return;
521  }
522 
523  // Initialization
524  PrepareTTree();
526 
527  bool overlay_2d = false;
528  for (size_t uci_index = 0; uci_index < _ucluster_v.size(); ++uci_index) {
529 
530  for (size_t vci_index = 0; vci_index < _vcluster_v.size(); ++vci_index) {
531 
532  overlay_2d = true;
533 
534  // Apply cuts
535  // Rough z-position overlay cut
536  if (_match_methods[kRoughZ]) {
537 
538  if (Match_RoughZ(_ucluster_v.at(uci_index), _vcluster_v.at(vci_index), geo::kU, geo::kV))
539  _tot_pass_z++;
540  else if (!_debug_mode)
541  continue;
542  else
543  overlay_2d = false;
544  }
545 
546  // Sum charge cut
547  if (_match_methods[kSumCharge]) {
548 
549  if (Match_SumCharge(_ucluster_v.at(uci_index), _vcluster_v.at(vci_index)))
550  _tot_pass_qsum++;
551  else if (!_debug_mode)
552  continue;
553  else
554  overlay_2d = false;
555  }
556 
557  // Rough time overlap cut
558  if (_match_methods[kRoughT]) {
559 
560  if (Match_RoughTime(_ucluster_v.at(uci_index), _vcluster_v.at(vci_index)))
561  _tot_pass_t++;
562  else if (!_debug_mode)
563  continue;
564  else
565  overlay_2d = false;
566  }
567 
568  // SpacePoint cut
569  std::vector<recob::SpacePoint> sps_v;
571 
572  if (Match_SpacePoint(clockData, detProp, uci_index, vci_index, 0, sps_v))
573  _tot_pass_sps++;
574  else if (!_debug_mode)
575  continue;
576  else
577  overlay_2d = false;
578  }
579 
580  if (overlay_2d) {
581  _matched_uclusters_v.push_back((unsigned int)(_ucluster_v[uci_index].cluster_index));
582  _matched_vclusters_v.push_back((unsigned int)(_vcluster_v[vci_index].cluster_index));
583  if (_store_sps) _matched_sps_v.push_back(sps_v);
584  }
585  } // end of ... _vcluster_v loop
586  } // end of ... _ucluster_v loop
587 
588  // Report
589  msg << std::endl
590  << Form("Found %zu matched cluster pairs...", _matched_uclusters_v.size()) << std::endl;
591  for (size_t i = 0; i < _matched_uclusters_v.size(); ++i) {
592 
593  if (i == 0) msg << "Listing matched clusters (U,V)..." << std::endl;
594 
595  msg << Form("Pair %-2zu: (%-3d, %-3d)", i, _matched_uclusters_v[i], _matched_vclusters_v[i])
596  << std::endl;
597  }
598  msg << std::endl;
599  mf::LogWarning("ClusterMatchAlg") << msg.str();
600 
601  if (_match_tree) _match_tree->Fill();
602  if (_cluster_tree) _cluster_tree->Fill();
603 
604  // Clear input event data holders
607  ClearTTreeInfo();
608  }
609 
611  detinfo::DetectorPropertiesData const& det_prop)
612  {
613  std::ostringstream msg;
614  msg << Form("Received (U,V,W) = (%zu,%zu,%zu) clusters...",
615  _uhits_v.size(),
616  _vhits_v.size(),
617  _whits_v.size())
618  << std::endl;
619  _tot_u = _ucluster_v.size();
620  _tot_v = _vcluster_v.size();
621  _tot_w = _wcluster_v.size();
622 
623  if (!(_tot_u + _tot_v + _tot_w)) {
624 
625  mf::LogError(__PRETTY_FUNCTION__)
626  << "No input cluster info found! Aborting the function call...";
627 
628  return;
629  }
630 
631  // Clear match information
632  PrepareTTree();
634 
635  bool overlay_2d = true;
636  bool overlay_3d = true;
637  // Loop over all possible u-v-w cluster combination
638  for (size_t uci_index = 0; uci_index < _ucluster_v.size(); ++uci_index) {
639 
640  for (size_t vci_index = 0; vci_index < _vcluster_v.size(); ++vci_index) {
641 
642  // Apply cuts that can be done with U&V planes here
643  overlay_2d = true;
644 
645  // Rough z-position overlay cut
646  if (_match_methods[kRoughZ]) {
647 
648  if (Match_RoughZ(_ucluster_v.at(uci_index), _vcluster_v.at(vci_index), geo::kU, geo::kV))
649  _tot_pass_z++;
650  else if (!_debug_mode)
651  continue;
652  else
653  overlay_2d = false;
654  }
655 
656  // Sum charge cut
657  if (_match_methods[kSumCharge]) {
658 
659  if (Match_SumCharge(_ucluster_v.at(uci_index), _vcluster_v.at(vci_index)))
660  _tot_pass_qsum++;
661  else if (!_debug_mode)
662  continue;
663  else
664  overlay_2d = false;
665  }
666 
667  for (size_t wci_index = 0; wci_index < _wcluster_v.size(); ++wci_index) {
668 
669  overlay_3d = overlay_2d;
670  // Apply cuts that requires 3 planes here
671 
672  // Rough time overlap cut
673  if (_match_methods[kRoughT]) {
674 
675  bool rough_time_match =
676  Match_RoughTime(_ucluster_v.at(uci_index), _vcluster_v.at(vci_index));
677  if (!_debug_mode && !rough_time_match) continue;
678 
679  rough_time_match =
680  (Match_RoughTime(_vcluster_v.at(vci_index), _wcluster_v.at(wci_index)) &&
681  rough_time_match);
682  if (!_debug_mode && !rough_time_match) continue;
683 
684  rough_time_match =
685  (Match_RoughTime(_wcluster_v.at(wci_index), _ucluster_v.at(uci_index)) &&
686  rough_time_match);
687 
688  overlay_3d = overlay_3d && rough_time_match;
689  if (rough_time_match)
690  _tot_pass_t++;
691  else if (!_debug_mode)
692  continue;
693  }
694 
695  // SpacePoint cut
696  std::vector<recob::SpacePoint> sps_v;
698 
699  if (Match_SpacePoint(clock_data, det_prop, uci_index, vci_index, wci_index, sps_v))
700  _tot_pass_sps++;
701  else if (!_debug_mode)
702  continue;
703  else
704  overlay_3d = false;
705  }
706 
707  if (overlay_3d) {
708  _matched_uclusters_v.push_back((unsigned int)(_ucluster_v[uci_index].cluster_index));
709  _matched_vclusters_v.push_back((unsigned int)(_vcluster_v[vci_index].cluster_index));
710  _matched_wclusters_v.push_back((unsigned int)(_wcluster_v[wci_index].cluster_index));
711  if (_store_sps) _matched_sps_v.push_back(sps_v);
712  }
713  } // end of ... _wcluster_v loop
714  } // end of ... _vcluster_v loop
715  } // end of ... _ucluster_v loop
716 
717  // Report
718  msg << std::endl
719  << Form("Found %zu matched cluster pairs...", _matched_uclusters_v.size()) << std::endl;
720  for (size_t i = 0; i < _matched_uclusters_v.size(); ++i) {
721 
722  if (i == 0) msg << "Listing matched clusters (U,V,W)..." << std::endl;
723 
724  msg << Form("Pair %-2zu: (%-3d, %-3d, %-3d)",
725  i,
729  << std::endl;
730  }
731  msg << std::endl;
732  mf::LogWarning("ClusterMatchAlg") << msg.str();
733 
734  if (_match_tree) _match_tree->Fill();
735  if (_cluster_tree) _cluster_tree->Fill();
736 
738  ClearTTreeInfo();
739  }
740 
741 } // namespace match
std::vector< uint16_t > _w_nhits_v
Use summed charge comparison ... see Match_SumCharge() description.
void ClearMatchInputInfo()
Method to clear input cluster information.
void reserve(size_type n)
Definition: PtrVector.h:337
std::vector< std::vector< unsigned int > > GetMatchedClusters() const
Method to retrieve matched cluster combinations. The format is [wire_plane][cluster_index].
trkf::SpacePointAlg * _sps_algo
SpacePointFinder algorithm pointer.
double peak_time_max
Maximum "peak time" among all hits in this cluster.
std::vector< art::PtrVector< recob::Hit > > _vhits_v
Local Hit pointer vector container ... V-plane.
void FillMCInfo(const art::Event &evt)
Internal method to fill MCTruth information when available.
Utilities related to art service access.
bool Match_RoughTime(const cluster_match_info &ci1, const cluster_match_info &ci2)
Checks min/max hit timing among two clusters and make sure there is an overlap.
double end_time_min
Minimum "end time" among all hits in this cluster.
std::vector< double > _vw_tratio_v
void AppendClusterInfo(detinfo::DetectorPropertiesData const &det_prop, const recob::Cluster &in_cluster, const std::vector< art::Ptr< recob::Hit >> &in_hit_v)
Method to fill cluster information to be used for matching.
unsigned short wire_min
Minimum wire number in this cluster.
std::vector< art::PtrVector< recob::Hit > > _whits_v
Local Hit pointer vector container ... W-plane.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Planes which measure V.
Definition: geo_types.h:136
double GetXTicksOffset(int p, int t, int c) const
double end_time_max
Maximum "end time" among all hits in this cluster.
bool Match_SpacePoint(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const size_t uindex, const size_t vindex, const size_t windex, std::vector< recob::SpacePoint > &sps_v)
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
std::vector< double > _wu_tratio_v
std::vector< double > _qratio_v
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
void ClearEventInfo()
Method to clear event-wise information.
void ClearMatchOutputInfo()
Method to clear output matched cluster information.
bool _store_sps
Boolean to enable storage of SpacePoint vector.
bool Match_RoughZ(const cluster_match_info &ci1, const cluster_match_info &ci2, const geo::View_t v1, const geo::View_t v2) const
Set of hits with a 2D structure.
Definition: Cluster.h:69
Use SpacePoint finder algorithm ... see Match_SpacePoint() description.
std::vector< cluster_match_info > _wcluster_v
Local cluster data container... W-plane.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
Particle class.
unsigned short _tot_pass_sps
std::vector< double > _tpeak_max_v
std::vector< double > _charge_v
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::vector< uint16_t > _nhits_v
Planes which measure U.
Definition: geo_types.h:135
TString part[npart]
Definition: Style.C:32
void PrepareTTree()
Internal method to create output TTree for quality checking of the algorithm.
std::vector< double > _tend_max_v
std::vector< double > _tend_min_v
std::vector< uint16_t > _nsps
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
double peak_time_min
Minimum "peak time" among all hits in this cluster.
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
T get(std::string const &key) const
Definition: ParameterSet.h:314
double sum_charge
Summed charge among all hits in this cluster.
Rough-Time comparison method ... see Match_RoughTime() description.
std::string _ModName_MCTruth
MCTruth producer&#39;s module name.
std::vector< art::PtrVector< recob::Hit > > _uhits_v
Local Hit pointer vector container ... U-plane.
double maxDT() const noexcept
Definition: SpacePointAlg.h:88
std::vector< uint16_t > _u_nhits_v
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
Declaration of cluster object.
void AppendClusterTreeVariables(const cluster_match_info &ci)
Internal method to fill cluster-info tree.
std::vector< double > _tstart_max_v
std::vector< cluster_match_info > _vcluster_v
Local cluster data container... V-plane.
std::vector< double > _tstart_min_v
std::vector< double > _uv_tratio_v
bool _debug_mode
Boolean to enable debug mode (call all enabled matching methods)
ClusterMatchAlg(fhicl::ParameterSet const &pset)
Default constructor with fhicl parameters.
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
void ClearTTreeInfo()
Method to clear TTree variables.
geo::View_t View() const
Returns the view for this cluster.
Definition: Cluster.h:714
void clearHitMap() const
bool Match_SumCharge(const cluster_match_info &uc, const cluster_match_info &vc)
std::vector< uint16_t > _view_v
unsigned short _tot_pass_qsum
Contains all timing reference information for the detector.
ID_t ID() const
Identifier of this cluster.
Definition: Cluster.h:711
unsigned short wire_max
Maximum wire number in this cluster.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< uint16_t > _v_nhits_v
std::vector< std::vector< recob::SpacePoint > > _matched_sps_v
Local SpacePoint vector container.
std::vector< unsigned int > _matched_uclusters_v
U plane matched clusters&#39; index.
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
double start_time_max
Maximum "start time" among all hits in this cluster.
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:137
bool _match_methods[kMATCH_METHOD_MAX]
Boolean list for enabled algorithms.
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
TCEvent evt
Definition: DataStructs.cxx:8
void MatchTwoPlanes(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop)
Two plane version of cluster matching method.
Event generator information.
Definition: MCTruth.h:32
void ReportConfig() const
Method to report the current configuration.
Float_t e
Definition: plot.C:35
Algorithm for generating space points from hits.
Namespace collecting geometry-related classes utilities.
void PrepareDetParams(detinfo::DetectorPropertiesData const &clockData)
Internal method, called only once, to fill detector-wise information.
std::vector< unsigned int > _matched_wclusters_v
W plane matched clusters&#39; index.
std::vector< unsigned int > _matched_vclusters_v
V plane matched clusters&#39; index.
void MatchThreePlanes(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop)
void FillHitInfo(cluster_match_info &ci, art::PtrVector< recob::Hit > &out_hit_v, const std::vector< art::Ptr< recob::Hit >> &in_hit_v)
Rough-Z comparison method ... see Match_RoughZ() description.
art framework interface to geometry description
std::vector< double > _tpeak_min_v
std::vector< cluster_match_info > _ucluster_v
Local cluster data container... U-plane.
double start_time_min
Minimum "start time" among all hits in this cluster.