LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
EmTrackClusterId2out_module.cc
Go to the documentation of this file.
1 // Class: EmTrackClusterId2out
3 // Module Type: producer
4 // File: EmTrackClusterId2out_module.cc
5 // Authors: dorota.stefan@cern.ch pplonski86@gmail.com robert.sulej@cern.ch
6 //
7 // Module applies CNN to 2D image made of deconvoluted wire waveforms in order
8 // to distinguish EM-like activity from track-like objects. New clusters of
9 // hits are produced to include also unclustered hits and tag everything in
10 // a common way.
11 // NOTE: This module uses 2-output CNN models, see EmTrackClusterId and
12 // EmTrackMichelClusterId for usage of 3 and 4-output models.
13 //
15 
26 
39 #include "cetlib_except/exception.h"
40 #include "fhiclcpp/ParameterSet.h"
41 #include "fhiclcpp/types/Atom.h"
42 #include "fhiclcpp/types/Comment.h"
43 #include "fhiclcpp/types/Name.h"
45 #include "fhiclcpp/types/Table.h"
47 
48 #include <iostream>
49 #include <memory>
50 #include <unordered_map>
51 #include <utility>
52 #include <vector>
53 
54 namespace nnet {
55 
57  public:
58  // these types to be replaced with use of feature proposed in redmine #12602
59  typedef std::unordered_map<unsigned int, std::vector<size_t>> view_keymap;
60  typedef std::unordered_map<unsigned int, view_keymap> tpc_view_keymap;
61  typedef std::unordered_map<unsigned int, tpc_view_keymap> cryo_tpc_view_keymap;
62 
63  struct Config {
64  using Name = fhicl::Name;
66 
69  Comment("number of samples processed in one batch")};
70 
72  Name("WireLabel"),
73  Comment("tag of deconvoluted ADC on wires (recob::Wire)")};
74 
76  Comment("tag of hits to be EM/track tagged")};
77 
79  Name("ClusterModuleLabel"),
80  Comment("tag of clusters to be used as a source of EM/track tagged new clusters (incl. "
81  "single-hit clusters ) using accumulated results from hits")};
82 
84  Name("TrackModuleLabel"),
85  Comment("tag of 3D tracks to be EM/track tagged using accumulated results from hits in the "
86  "best 2D projection")};
87 
89  Name("Views"),
90  Comment("tag clusters in selected views only, or in all views if empty list")};
91  };
93  explicit EmTrackClusterId2out(Parameters const& p);
94 
99 
100  private:
101  void produce(art::Event& e) override;
102 
103  bool isViewSelected(int view) const;
104 
105  size_t fBatchSize;
107  anab::MVAWriter<2> fMVAWriter; // <-------------- using 2-output CNN model
108 
114 
115  std::vector<int> fViews;
116 
117  art::InputTag fNewClustersTag; // input tag for the clusters produced by this module
118  };
119  // ------------------------------------------------------
120 
122  : EDProducer{config}
123  , fBatchSize(config().BatchSize())
124  , fPointIdAlg(config().PointIdAlg())
125  , fMVAWriter(producesCollector(), "emtrack")
126  , fWireProducerLabel(config().WireLabel())
127  , fHitModuleLabel(config().HitModuleLabel())
128  , fClusterModuleLabel(config().ClusterModuleLabel())
129  , fTrackModuleLabel(config().TrackModuleLabel())
130  , fViews(config().Views())
131  ,
132 
133  fNewClustersTag(config.get_PSet().get<std::string>("module_label"),
134  "",
136  {
138 
139  if (!fClusterModuleLabel.label().empty()) {
140  produces<std::vector<recob::Cluster>>();
141  produces<art::Assns<recob::Cluster, recob::Hit>>();
142 
144  fDoClusters = true;
145  }
146  else {
147  fDoClusters = false;
148  }
149 
150  if (!fTrackModuleLabel.label().empty()) {
152  fDoTracks = true;
153  }
154  else {
155  fDoTracks = false;
156  }
157  }
158  // ------------------------------------------------------
159 
161  {
162  mf::LogVerbatim("EmTrackClusterId2out")
163  << "next event: " << evt.run() << " / " << evt.id().event();
164 
165  auto wireHandle = evt.getValidHandle<std::vector<recob::Wire>>(fWireProducerLabel);
166 
167  unsigned int cryo, tpc, view;
168 
169  // ******************* get and sort hits ********************
170  auto hitListHandle = evt.getValidHandle<std::vector<recob::Hit>>(fHitModuleLabel);
171  std::vector<art::Ptr<recob::Hit>> hitPtrList;
172  art::fill_ptr_vector(hitPtrList, hitListHandle);
173 
175  for (auto const& h : hitPtrList) {
176  view = h->WireID().Plane;
177  if (!isViewSelected(view)) continue;
178 
179  cryo = h->WireID().Cryostat;
180  tpc = h->WireID().TPC;
181 
182  hitMap[cryo][tpc][view].push_back(h.key());
183  }
184 
185  // ********************* classify hits **********************
186  auto hitID = fMVAWriter.initOutputs<recob::Hit>(
187  fHitModuleLabel, hitPtrList.size(), fPointIdAlg.outputLabels());
188 
189  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
190  auto const detProp =
192 
193  std::vector<char> hitInFA(
194  hitPtrList.size(),
195  0); // tag hits in fid. area as 1, use 0 for hits close to the projectrion edges
196  for (auto const& pcryo : hitMap) {
197  cryo = pcryo.first;
198  for (auto const& ptpc : pcryo.second) {
199  tpc = ptpc.first;
200  for (auto const& pview : ptpc.second) {
201  view = pview.first;
202  if (!isViewSelected(view)) continue; // should not happen, hits were selected
203 
204  fPointIdAlg.setWireDriftData(clockData, detProp, *wireHandle, view, tpc, cryo);
205 
206  // (1) do all hits in this plane ------------------------------------------------
207  for (size_t idx = 0; idx < pview.second.size(); idx += fBatchSize) {
208  std::vector<std::pair<unsigned int, float>> points;
209  std::vector<size_t> keys;
210  for (size_t k = 0; k < fBatchSize; ++k) {
211  if (idx + k >= pview.second.size()) { break; } // careful about the tail
212 
213  size_t h = pview.second[idx + k]; // h is the Ptr< recob::Hit >::key()
214  const recob::Hit& hit = *(hitPtrList[h]);
215  points.emplace_back(hit.WireID().Wire, hit.PeakTime());
216  keys.push_back(h);
217  }
218 
219  auto batch_out = fPointIdAlg.predictIdVectors(points);
220  if (points.size() != batch_out.size()) {
221  throw cet::exception("EmTrackClusterId") << "hits processing failed" << std::endl;
222  }
223 
224  for (size_t k = 0; k < points.size(); ++k) {
225  size_t h = keys[k];
226  fMVAWriter.setOutput(hitID, h, batch_out[k]);
227  if (fPointIdAlg.isInsideFiducialRegion(points[k].first, points[k].second)) {
228  hitInFA[h] = 1;
229  }
230  }
231  } // hits done ------------------------------------------------------------------
232  }
233  }
234  }
235 
236  // (2) do clusters when hits are ready in all planes ----------------------------------------
237  if (fDoClusters) {
238  // **************** prepare for new clusters ****************
239  auto clusters = std::make_unique<std::vector<recob::Cluster>>();
240  auto clu2hit = std::make_unique<art::Assns<recob::Cluster, recob::Hit>>();
241 
242  // ************** get and sort input clusters ***************
243  auto cluListHandle = evt.getValidHandle<std::vector<recob::Cluster>>(fClusterModuleLabel);
244  std::vector<art::Ptr<recob::Cluster>> cluPtrList;
245  art::fill_ptr_vector(cluPtrList, cluListHandle);
246 
248  for (auto const& c : cluPtrList) {
249  view = c->Plane().Plane;
250  if (!isViewSelected(view)) continue;
251 
252  cryo = c->Plane().Cryostat;
253  tpc = c->Plane().TPC;
254 
255  cluMap[cryo][tpc][view].push_back(c.key());
256  }
257 
258  auto cluID =
260 
261  unsigned int cidx = 0; // new clusters index
262  art::FindManyP<recob::Hit> hitsFromClusters(cluListHandle, evt, fClusterModuleLabel);
263  std::vector<bool> hitUsed(hitPtrList.size(), false); // tag hits used in clusters
264  for (auto const& pcryo : cluMap) {
265  cryo = pcryo.first;
266  for (auto const& ptpc : pcryo.second) {
267  tpc = ptpc.first;
268  for (auto const& pview : ptpc.second) {
269  view = pview.first;
270  if (!isViewSelected(view)) continue; // should not happen, clusters were pre-selected
271 
272  for (size_t c : pview.second) // c is the Ptr< recob::Cluster >::key()
273  {
274  auto v = hitsFromClusters.at(c);
275  if (v.empty()) continue;
276 
277  for (auto const& hit : v) {
278  if (hitUsed[hit.key()]) {
279  mf::LogWarning("EmTrackClusterId2out") << "hit already used in another cluster";
280  }
281  hitUsed[hit.key()] = true;
282  }
283 
284  auto vout = fMVAWriter.getOutput<recob::Hit>(
285  v, [&](art::Ptr<recob::Hit> const& ptr) { return (float)hitInFA[ptr.key()]; });
286 
287  mf::LogVerbatim("EmTrackClusterId2out")
288  << "cluster in tpc:" << tpc << " view:" << view << " size:" << v.size()
289  << " p:" << vout[0];
290 
291  clusters->emplace_back(recob::Cluster(0.0F,
292  0.0F,
293  0.0F,
294  0.0F,
295  0.0F,
296  0.0F,
297  0.0F,
298  0.0F,
299  0.0F,
300  0.0F,
301  0.0F,
302  0.0F,
303  0.0F,
304  0.0F,
305  0.0F,
306  0.0F,
307  0.0F,
308  0.0F,
309  v.size(),
310  0.0F,
311  0.0F,
312  cidx,
313  (geo::View_t)view,
314  v.front()->WireID().planeID()));
315  util::CreateAssn(evt, *clusters, v, *clu2hit);
316  cidx++;
317 
318  fMVAWriter.addOutput(cluID, vout); // add copy of the input cluster
319  }
320 
321  // (2b) make single-hit clusters --------------------------------------------
322  for (size_t h : hitMap[cryo][tpc][view]) // h is the Ptr< recob::Hit >::key()
323  {
324  if (hitUsed[h]) continue;
325 
326  auto vout = fMVAWriter.getOutput<recob::Hit>(h);
327 
328  mf::LogVerbatim("EmTrackClusterId2out")
329  << "single hit in tpc:" << tpc << " view:" << view
330  << " wire:" << hitPtrList[h]->WireID().Wire
331  << " drift:" << hitPtrList[h]->PeakTime() << " p:" << vout[0];
332 
333  art::PtrVector<recob::Hit> cluster_hits;
334  cluster_hits.push_back(hitPtrList[h]);
335  clusters->emplace_back(recob::Cluster(0.0F,
336  0.0F,
337  0.0F,
338  0.0F,
339  0.0F,
340  0.0F,
341  0.0F,
342  0.0F,
343  0.0F,
344  0.0F,
345  0.0F,
346  0.0F,
347  0.0F,
348  0.0F,
349  0.0F,
350  0.0F,
351  0.0F,
352  0.0F,
353  1,
354  0.0F,
355  0.0F,
356  cidx,
357  (geo::View_t)view,
358  hitPtrList[h]->WireID().planeID()));
359  util::CreateAssn(evt, *clusters, cluster_hits, *clu2hit);
360  cidx++;
361 
362  fMVAWriter.addOutput(cluID, vout); // add single-hit cluster tagging unclutered hit
363  }
364  mf::LogVerbatim("EmTrackClusterId2out")
365  << "...produced " << cidx - pview.second.size() << " single-hit clusters.";
366  }
367  }
368  }
369 
370  evt.put(std::move(clusters));
371  evt.put(std::move(clu2hit));
372  } // all clusters done ----------------------------------------------------------------------
373 
374  // (3) do tracks when all hits in all cryo/tpc/plane are done -------------------------------
375  if (fDoTracks) {
376  auto trkListHandle = evt.getValidHandle<std::vector<recob::Track>>(fTrackModuleLabel);
377  art::FindManyP<recob::Hit> hitsFromTracks(trkListHandle, evt, fTrackModuleLabel);
378  std::vector<std::vector<art::Ptr<recob::Hit>>> trkHitPtrList(trkListHandle->size());
379  for (size_t t = 0; t < trkListHandle->size(); ++t) {
380  auto v = hitsFromTracks.at(t);
381  size_t nh[3] = {0, 0, 0};
382  for (auto const& hptr : v) {
383  ++nh[hptr->View()];
384  }
385  size_t best_view = 2; // collection
386  if ((nh[0] >= nh[1]) && (nh[0] > 2 * nh[2])) best_view = 0; // ind1
387  if ((nh[1] >= nh[0]) && (nh[1] > 2 * nh[2])) best_view = 1; // ind2
388 
389  size_t k = 0;
390  while (!isViewSelected(best_view)) {
391  best_view = (best_view + 1) % 3;
392  if (++k > 3) {
393  throw cet::exception("EmTrackClusterId2out")
394  << "No views selected at all?" << std::endl;
395  }
396  }
397 
398  for (auto const& hptr : v) {
399  if (hptr->View() == best_view) trkHitPtrList[t].emplace_back(hptr);
400  }
401  }
402 
403  auto trkID = fMVAWriter.initOutputs<recob::Track>(
404  fTrackModuleLabel, trkHitPtrList.size(), fPointIdAlg.outputLabels());
405  for (size_t t = 0; t < trkHitPtrList.size(); ++t) // t is the Ptr< recob::Track >::key()
406  {
407  auto vout =
408  fMVAWriter.getOutput<recob::Hit>(trkHitPtrList[t], [&](art::Ptr<recob::Hit> const& ptr) {
409  return (float)hitInFA[ptr.key()];
410  });
411  fMVAWriter.setOutput(trkID, t, vout);
412  }
413  }
414  // tracks done ------------------------------------------------------------------------------
415 
416  fMVAWriter.saveOutputs(evt);
417  }
418  // ------------------------------------------------------
419 
421  {
422  if (fViews.empty())
423  return true;
424  else {
425  for (auto k : fViews)
426  if (k == view) { return true; }
427  return false;
428  }
429  }
430  // ------------------------------------------------------
431 
433 
434 }
bool isInsideFiducialRegion(unsigned int wire, float drift) const
Definition: PointIdAlg.cxx:330
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< std::string > const & outputLabels() const
network output labels
Definition: PointIdAlg.h:121
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
fhicl::Atom< art::InputTag > ClusterModuleLabel
void setOutput(FVector_ID id, size_t key, std::array< float, N > const &values)
Definition: MVAWriter.h:214
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Set of hits with a 2D structure.
Definition: Cluster.h:69
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
void produce(art::Event &e) override
fhicl::Atom< art::InputTag > TrackModuleLabel
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
std::string const & label() const noexcept
Definition: InputTag.cc:79
FVector_ID initOutputs(std::string const &dataTag, size_t dataSize, std::vector< std::string > const &names=std::vector< std::string >(N,""))
EmTrackClusterId2out & operator=(EmTrackClusterId2out const &)=delete
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
std::array< float, N > getOutput(std::vector< art::Ptr< T >> const &items) const
Definition: MVAWriter.h:251
Definition: EmTrack.h:40
Provides recob::Track data product.
Declaration of cluster object.
Definition of data types for geometry description.
void saveOutputs(art::Event &evt)
Check consistency and save all the results in the event.
Definition: MVAWriter.h:404
Detector simulation of raw signals on wires.
ProducesCollector & producesCollector() noexcept
void addOutput(FVector_ID id, std::array< float, N > const &values)
Definition: MVAWriter.h:231
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
std::unordered_map< unsigned int, tpc_view_keymap > cryo_tpc_view_keymap
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
bool setWireDriftData(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const std::vector< recob::Wire > &wires, unsigned int plane, unsigned int tpc, unsigned int cryo)
Utility object to perform functions of association.
std::unordered_map< unsigned int, std::vector< size_t > > view_keymap
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
EventNumber_t event() const
Definition: EventID.h:116
Declaration of basic channel signal object.
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
constexpr PlaneID const & planeID() const
Definition: geo_types.h:620
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.cc:29
std::unordered_map< unsigned int, view_keymap > tpc_view_keymap
std::vector< std::vector< float > > predictIdVectors(std::vector< std::pair< unsigned int, float >> points) const
Definition: PointIdAlg.cxx:257
EventID id() const
Definition: Event.cc:23
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33