LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MCMatchAlg.cxx
Go to the documentation of this file.
1 #include "MCMatchAlg.h"
2 
3 #include "TString.h"
11 
12 #include <string>
13 
14 namespace btutil {
15 
17  {
18  _view_to_plane.clear();
19  }
20 
22  const std::vector<unsigned int>& g4_trackid_v,
23  const std::vector<sim::SimChannel>& simch_v,
24  const std::vector<std::vector<art::Ptr<recob::Hit>>>& cluster_v)
25  {
26  fBTAlgo.Reset(g4_trackid_v, simch_v);
27 
28  return BuildMap(clockData, cluster_v);
29  }
30 
32  const std::vector<std::vector<unsigned int>>& g4_trackid_v,
33  const std::vector<sim::SimChannel>& simch_v,
34  const std::vector<std::vector<art::Ptr<recob::Hit>>>& cluster_v)
35  {
36 
37  fBTAlgo.Reset(g4_trackid_v, simch_v);
38 
39  return BuildMap(clockData, cluster_v);
40  }
41 
43  const std::vector<std::vector<art::Ptr<recob::Hit>>>& cluster_v)
44  {
45  size_t num_mcobj = fBTAlgo.NumParts();
46  size_t num_cluster = cluster_v.size();
47  //auto geo = ::larutil::Geometry::GetME();
49 
50  //
51  // Perform back-tracking
52  //
53  // (1) Load cluster/hit data product
54  // (2) Loop over all clusters and find charge fraction to each MC object
55  // (3) For each MC object, find a cluster per plane with the highest charge
56 
57  // Loop over clusters & get charge fraction per MC object
58  _summed_mcq.clear();
59  _cluster_mcq_v.clear();
60  _cluster_plane_id.clear();
61 
62  _summed_mcq.resize(num_mcobj + 1, std::vector<double>(geo->Nplanes(), 0));
63  _cluster_mcq_v.reserve(num_cluster);
64  _cluster_plane_id.reserve(num_cluster);
65 
66  for (auto const& hit_v : cluster_v) {
67 
68  size_t plane = geo->Nplanes();
69 
70  // Create hit list
71  std::vector<WireRange_t> wr_v;
72  wr_v.reserve(hit_v.size());
73 
74  for (auto const& h : hit_v) {
75 
76  WireRange_t wr;
77  wr.ch = h->Channel();
78  wr.start = h->StartTick();
79  wr.end = h->EndTick();
80  wr_v.push_back(wr);
81  if (plane == geo->Nplanes()) plane = h->WireID().Plane;
82  //if(plane==geo->Nplanes()) plane = h->View();
83  }
84 
85  _cluster_plane_id.push_back(plane);
86 
87  auto mcq_v = fBTAlgo.MCQ(clockData, wr_v);
88 
89  for (size_t i = 0; i < mcq_v.size(); ++i)
90 
91  _summed_mcq[i][plane] += mcq_v[i];
92 
93  _cluster_mcq_v.push_back(mcq_v);
94  }
95 
96  //
97  // Find the best matched pair (and its charge) per MC object
98  //
99  _bmatch_id.clear();
100  _bmatch_id.resize(num_mcobj, std::vector<int>(geo->Nplanes(), -1));
101 
102  std::vector<std::vector<double>> bmatch_mcq(num_mcobj, std::vector<double>(geo->Nplanes(), 0));
103 
104  for (size_t mc_index = 0; mc_index < num_mcobj; ++mc_index) {
105 
106  for (size_t cluster_index = 0; cluster_index < num_cluster; ++cluster_index) {
107 
108  //if((_cluster_mcq_v[cluster_index].back()) < 0) continue;
109 
110  auto plane = _cluster_plane_id.at(cluster_index);
111 
112  double q = _cluster_mcq_v[cluster_index][mc_index];
113 
114  if (bmatch_mcq[mc_index][plane] < q) {
115 
116  bmatch_mcq[mc_index][plane] = q;
117  _bmatch_id[mc_index][plane] = cluster_index;
118  }
119  }
120  }
121 
122  return true;
123  }
124 
125  double MCMatchAlg::ClusterCorrectness(const size_t cluster_index, const size_t mc_index) const
126  {
127 
128  if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");
129 
130  if (cluster_index >= _cluster_mcq_v.size())
131  throw MCBTException(Form(
132  "Input cluster index (%zu) out of range (%zu)!", cluster_index, _cluster_mcq_v.size()));
133 
134  if (mc_index >= _bmatch_id.size())
135  throw MCBTException(
136  Form("Input MC object index (%zu) out of range (%zu)!", mc_index, _bmatch_id.size()));
137 
138  auto plane = _cluster_plane_id.at(cluster_index);
139 
140  auto best_cluster_index = _bmatch_id.at(mc_index).at(plane);
141 
142  if (best_cluster_index < 0) return -1;
143 
144  return _cluster_mcq_v.at(cluster_index).at(mc_index) /
145  _cluster_mcq_v.at(best_cluster_index).at(mc_index);
146  }
147 
148  std::pair<size_t, double> MCMatchAlg::ShowerCorrectness(
149  const std::vector<unsigned int> cluster_indices) const
150  {
151 
152  if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");
153 
154  if (cluster_indices.size() > _cluster_mcq_v.size())
155  throw MCBTException(
156  Form("Input cluster indices length (%zu) > # of registered clusters (%zu)!",
157  cluster_indices.size(),
158  _cluster_mcq_v.size()));
159 
160  if (!cluster_indices.size()) throw MCBTException("Input cluster indices empty!");
161 
162  // Compute efficiency per MC
163  std::vector<double> match_eff(_bmatch_id.size(), 1);
164 
165  for (auto const& cluster_index : cluster_indices) {
166 
167  for (size_t mc_index = 0; mc_index < _bmatch_id.size(); ++mc_index) {
168 
169  double correctness = ClusterCorrectness(cluster_index, mc_index);
170 
171  if (correctness >= 0) match_eff.at(mc_index) *= correctness;
172  }
173  }
174 
175  std::pair<size_t, double> result(0, -1);
176 
177  // Find the best qratio
178  for (size_t mc_index = 0; mc_index < match_eff.size(); ++mc_index) {
179 
180  if (match_eff.at(mc_index) < result.second) continue;
181 
182  result.second = match_eff.at(mc_index);
183 
184  result.first = mc_index;
185  }
186  return result;
187  }
188 
189  std::pair<double, double> MCMatchAlg::ClusterEP(const size_t cluster_index,
190  const size_t mc_index) const
191  {
192  if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");
193 
194  if (cluster_index >= _cluster_mcq_v.size())
195  throw MCBTException(Form(
196  "Input cluster index (%zu) out of range (%zu)!", cluster_index, _cluster_mcq_v.size()));
197 
198  if (mc_index >= _bmatch_id.size())
199  throw MCBTException(
200  Form("Input MC object index (%zu) out of range (%zu)!", mc_index, _bmatch_id.size()));
201 
202  // efficiency = this cluster's mcq for specified mc obj / total mcq by this mcshower in all clusters
203  // purity = this cluster's mcq for this mcshower / total mcq from all mc obj in this cluster
204 
205  std::pair<double, double> result;
206 
207  auto plane = _cluster_plane_id[cluster_index];
208 
209  result.first = _cluster_mcq_v[cluster_index][mc_index] / _summed_mcq[mc_index][plane];
210 
211  double cluster_mcq_total = 0;
212  for (auto const& q : _cluster_mcq_v[cluster_index])
213  cluster_mcq_total += q;
214 
215  result.second = _cluster_mcq_v[cluster_index][mc_index] / cluster_mcq_total;
216 
217  return result;
218  }
219 
220  const std::vector<int>& MCMatchAlg::BestClusters(const size_t mc_index) const
221  {
222  if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");
223 
224  if (mc_index >= _bmatch_id.size())
225  throw MCBTException(
226  Form("Input MC object index (%zu) out of range (%zu)!", mc_index, _bmatch_id.size()));
227 
228  return _bmatch_id[mc_index];
229  }
230 
231  std::pair<double, double> MCMatchAlg::BestClusterEP(const size_t mc_index,
232  const size_t plane_id) const
233  {
234 
235  auto c_index_v = BestClusters(mc_index);
236 
237  if (c_index_v.size() <= plane_id)
238  throw MCBTException(Form(
239  "Plane ID %zu exceeds # of planes recorded in data (%zu)...", plane_id, c_index_v.size()));
240 
241  std::pair<double, double> result(0, 0);
242 
243  if (c_index_v[plane_id] < 0) return result;
244 
245  return ClusterEP(c_index_v[plane_id], mc_index);
246  }
247 
248 }
unsigned int ch
Definition: MCBTAlg.h:34
Class def header for a class MCMatchAlg.
std::pair< double, double > BestClusterEP(const size_t mcshower_index, const size_t plane_id) const
Definition: MCMatchAlg.cxx:231
void Reset(const std::vector< unsigned int > &g4_trackid_v, const std::vector< sim::SimChannel > &simch_v)
Definition: MCBTAlg.cxx:23
Declaration of signal hit object.
Class def header for a class MCBTAlg.
bool BuildMap(detinfo::DetectorClocksData const &clockData, const std::vector< unsigned int > &g4_trackid_v, const std::vector< sim::SimChannel > &simch_v, const std::vector< std::vector< art::Ptr< recob::Hit >>> &cluster_v)
Constructs needed information for Reco=>MC matching.
Definition: MCMatchAlg.cxx:21
std::vector< std::vector< double > > _summed_mcq
Definition: MCMatchAlg.h:103
std::vector< unsigned char > _cluster_plane_id
Definition: MCMatchAlg.h:106
std::vector< std::vector< double > > _cluster_mcq_v
Definition: MCMatchAlg.h:104
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
double ClusterCorrectness(const size_t cluster_index, const size_t mcshower_index) const
Definition: MCMatchAlg.cxx:125
MCBTAlg fBTAlgo
MCBTAlg instance.
Definition: MCMatchAlg.h:99
std::pair< double, double > ClusterEP(const size_t cluster_index, const size_t mcshower_index) const
For a specified cluster, compute cluster efficiency and purity in terms of specified MC object...
Definition: MCMatchAlg.cxx:189
std::vector< double > MCQ(detinfo::DetectorClocksData const &clockData, const WireRange_t &hit) const
Definition: MCBTAlg.cxx:101
size_t NumParts() const
Definition: MCBTAlg.h:115
Definition of data types for geometry description.
const std::vector< int > & BestClusters(const size_t mcshower_index) const
Definition: MCMatchAlg.cxx:220
std::vector< size_t > _view_to_plane
Definition: MCMatchAlg.h:101
Contains all timing reference information for the detector.
MCMatchAlg()
Default constructor.
Definition: MCMatchAlg.cxx:16
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
std::pair< size_t, double > ShowerCorrectness(const std::vector< unsigned int > cluster_indices) const
Definition: MCMatchAlg.cxx:148
std::vector< std::vector< int > > _bmatch_id
Definition: MCMatchAlg.h:107
Namespace collecting geometry-related classes utilities.
art framework interface to geometry description
Class def header for exception classes in MCComp package.