LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PointIdAlg.cxx
Go to the documentation of this file.
1 // Class: PointIdAlg
3 // Authors: D.Stefan (Dorota.Stefan@ncbj.gov.pl), from DUNE, CERN/NCBJ, since May 2016
4 // R.Sulej (Robert.Sulej@cern.ch), from DUNE, FNAL/NCBJ, since May 2016
5 // P.Plonski, from DUNE, WUT, since May 2016
6 // D.Smith, from LArIAT, BU, 2017: real data dump
8 
10 
11 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom<>()
12 #include "larcorealg/Geometry/Exceptions.h" // geo::InvalidWireIDError
13 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::InvalidChannelID
21 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
22 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
24 
31 #include "cetlib/search_path.h"
32 #include "cetlib_except/exception.h"
34 
35 #include "tensorflow/core/public/session.h"
36 
37 #include "TMath.h"
38 
39 #include <algorithm>
40 #include <cmath>
41 #include <map>
42 #include <string>
43 #include <sys/stat.h>
44 #include <unordered_map>
45 #include <vector>
46 
47 // ------------------------------------------------------
48 // -------------------ModelInterface---------------------
49 // ------------------------------------------------------
50 
51 std::vector<std::vector<float>> nnet::ModelInterface::Run(
52  std::vector<std::vector<std::vector<float>>> const& inps,
53  int samples)
54 {
55  if ((samples == 0) || inps.empty() || inps.front().empty() || inps.front().front().empty())
56  return std::vector<std::vector<float>>();
57 
58  if ((samples == -1) || (samples > (int)inps.size())) { samples = inps.size(); }
59 
60  std::vector<std::vector<float>> results;
61  for (int i = 0; i < samples; ++i) {
62  results.push_back(Run(inps[i]));
63  }
64  return results;
65 }
66 
67 std::string nnet::ModelInterface::findFile(const char* fileName) const
68 {
69  std::string fname_out;
70  cet::search_path sp("FW_SEARCH_PATH");
71  if (!sp.find_file(fileName, fname_out)) {
72  struct stat buffer;
73  if (stat(fileName, &buffer) == 0) { fname_out = fileName; }
74  else {
75  throw art::Exception(art::errors::NotFound) << "Could not find the model file " << fileName;
76  }
77  }
78  return fname_out;
79 }
80 
81 // ------------------------------------------------------
82 // ----------------KerasModelInterface-------------------
83 // ------------------------------------------------------
84 
86  : m(nnet::ModelInterface::findFile(modelFileName).c_str())
87 {
88  mf::LogInfo("KerasModelInterface") << "Keras model loaded.";
89 }
90 // ------------------------------------------------------
91 
92 std::vector<float> nnet::KerasModelInterface::Run(std::vector<std::vector<float>> const& inp2d)
93 {
94  std::vector<std::vector<std::vector<float>>> inp3d;
95  inp3d.push_back(inp2d); // lots of copy, should add 2D to keras...
96 
97  keras::DataChunk* sample = new keras::DataChunk2D();
98  sample->set_data(inp3d); // and more copy...
99  std::vector<float> out = m.compute_output(sample);
100  delete sample;
101  return out;
102 }
103 
104 // ------------------------------------------------------
105 // -----------------TfModelInterface---------------------
106 // ------------------------------------------------------
107 
109 {
110  g = tf::Graph::create(nnet::ModelInterface::findFile(modelFileName).c_str(),
111  {"cnn_output", "_netout"});
112  if (!g) { throw art::Exception(art::errors::Unknown) << "TF model failed."; }
113 
114  mf::LogInfo("TfModelInterface") << "TF model loaded.";
115 }
116 // ------------------------------------------------------
117 
118 std::vector<std::vector<float>> nnet::TfModelInterface::Run(
119  std::vector<std::vector<std::vector<float>>> const& inps,
120  int samples)
121 {
122  if ((samples == 0) || inps.empty() || inps.front().empty() || inps.front().front().empty())
123  return std::vector<std::vector<float>>();
124 
125  if ((samples == -1) || (samples > (long long int)inps.size())) { samples = inps.size(); }
126 
127  long long int rows = inps.front().size(), cols = inps.front().front().size();
128 
129  std::vector<tensorflow::Tensor> _x;
130  _x.push_back(
131  tensorflow::Tensor(tensorflow::DT_FLOAT, tensorflow::TensorShape({samples, rows, cols, 1})));
132  auto input_map = _x[0].tensor<float, 4>();
133  for (long long int s = 0; s < samples; ++s) {
134  const auto& sample = inps[s];
135  for (long long int r = 0; r < rows; ++r) {
136  const auto& row = sample[r];
137  for (long long int c = 0; c < cols; ++c) {
138  input_map(s, r, c, 0) = row[c];
139  }
140  }
141  }
142 
143  return g->runx(_x);
144 }
145 // ------------------------------------------------------
146 
147 std::vector<float> nnet::TfModelInterface::Run(std::vector<std::vector<float>> const& inp2d)
148 {
149  long long int rows = inp2d.size(), cols = inp2d.front().size();
150 
151  std::vector<tensorflow::Tensor> _x;
152  _x.push_back(
153  tensorflow::Tensor(tensorflow::DT_FLOAT, tensorflow::TensorShape({1, rows, cols, 1})));
154  auto input_map = _x[0].tensor<float, 4>();
155  for (long long int r = 0; r < rows; ++r) {
156  const auto& row = inp2d[r];
157  for (long long int c = 0; c < cols; ++c) {
158  input_map(0, r, c, 0) = row[c];
159  }
160  }
161 
162  auto out = g->runx(_x);
163  if (!out.empty())
164  return out.front();
165  else
166  return std::vector<float>();
167 }
168 
169 // ------------------------------------------------------
170 // --------------------PointIdAlg------------------------
171 // ------------------------------------------------------
172 
174  : img::DataProviderAlg(config)
175  , fNNet(0)
176  , fPatchSizeW(config.PatchSizeW())
177  , fPatchSizeD(config.PatchSizeD())
178  , fCurrentWireIdx(99999)
179  , fCurrentScaledDrift(99999)
180 {
182  fNNetOutputs = config.NNetOutputs();
183 
184  deleteNNet();
185 
186  if ((fNNetModelFilePath.length() > 5) &&
187  (fNNetModelFilePath.compare(fNNetModelFilePath.length() - 5, 5, ".nnet") == 0)) {
189  }
190  else if ((fNNetModelFilePath.length() > 3) &&
191  (fNNetModelFilePath.compare(fNNetModelFilePath.length() - 3, 3, ".pb") == 0)) {
193  }
194  else {
195  mf::LogError("PointIdAlg") << "File name extension not supported.";
196  }
197 
198  if (!fNNet) { throw cet::exception("nnet::PointIdAlg") << "Loading model from file failed."; }
199 
200  resizePatch();
201 }
202 // ------------------------------------------------------
203 
205 {
206  deleteNNet();
207 }
208 // ------------------------------------------------------
209 
211 {
213  for (auto& r : fWireDriftPatch)
214  r.resize(fPatchSizeD);
215 }
216 // ------------------------------------------------------
217 
218 float nnet::PointIdAlg::predictIdValue(unsigned int wire, float drift, size_t outIdx) const
219 {
220  float result = 0.;
221 
222  if (!bufferPatch(wire, drift)) {
223  mf::LogError("PointIdAlg") << "Patch buffering failed.";
224  return result;
225  }
226 
227  if (fNNet) {
228  auto out = fNNet->Run(fWireDriftPatch);
229  if (!out.empty()) { result = out[outIdx]; }
230  else {
231  mf::LogError("PointIdAlg") << "Problem with applying model to input.";
232  }
233  }
234 
235  return result;
236 }
237 // ------------------------------------------------------
238 
239 std::vector<float> nnet::PointIdAlg::predictIdVector(unsigned int wire, float drift) const
240 {
241  std::vector<float> result;
242 
243  if (!bufferPatch(wire, drift)) {
244  mf::LogError("PointIdAlg") << "Patch buffering failed.";
245  return result;
246  }
247 
248  if (fNNet) {
249  result = fNNet->Run(fWireDriftPatch);
250  if (result.empty()) { mf::LogError("PointIdAlg") << "Problem with applying model to input."; }
251  }
252 
253  return result;
254 }
255 // ------------------------------------------------------
256 
257 std::vector<std::vector<float>> nnet::PointIdAlg::predictIdVectors(
258  std::vector<std::pair<unsigned int, float>> points) const
259 {
260  if (points.empty() || !fNNet) { return std::vector<std::vector<float>>(); }
261 
262  std::vector<std::vector<std::vector<float>>> inps(
263  points.size(), std::vector<std::vector<float>>(fPatchSizeW, std::vector<float>(fPatchSizeD)));
264  for (size_t i = 0; i < points.size(); ++i) {
265  unsigned int wire = points[i].first;
266  float drift = points[i].second;
267  if (!bufferPatch(wire, drift, inps[i])) {
268  throw cet::exception("PointIdAlg") << "Patch buffering failed" << std::endl;
269  }
270  }
271 
272  return fNNet->Run(inps);
273 }
274 // ------------------------------------------------------
275 
276 bool nnet::PointIdAlg::isSamePatch(unsigned int wire1,
277  float drift1,
278  unsigned int wire2,
279  float drift2) const
280 {
281  if (fDownscaleFullView) {
282  size_t sd1 = (size_t)(drift1 / fDriftWindow);
283  size_t sd2 = (size_t)(drift2 / fDriftWindow);
284  if ((wire1 == wire2) && (sd1 == sd2)) return true; // the same position
285  }
286  else {
287  if ((wire1 == wire2) && ((size_t)drift1 == (size_t)drift2)) return true; // the same position
288  }
289 
290  return false; // not the same position
291 }
292 
293 bool nnet::PointIdAlg::isCurrentPatch(unsigned int wire, float drift) const
294 {
295  if (fDownscaleFullView) {
296  size_t sd = (size_t)(drift / fDriftWindow);
297  if ((fCurrentWireIdx == wire) && (fCurrentScaledDrift == sd))
298  return true; // still within the current position
299  }
300  else {
301  if ((fCurrentWireIdx == wire) && (fCurrentScaledDrift == drift))
302  return true; // still within the current position
303  }
304 
305  return false; // not a current position
306 }
307 // ------------------------------------------------------
308 
309 std::vector<float> nnet::PointIdAlg::flattenData2D(std::vector<std::vector<float>> const& patch)
310 {
311  std::vector<float> flat;
312  if (patch.empty() || patch.front().empty()) {
313  mf::LogError("DataProviderAlg") << "Patch is empty.";
314  return flat;
315  }
316 
317  flat.resize(patch.size() * patch.front().size());
318 
319  for (size_t w = 0, i = 0; w < patch.size(); ++w) {
320  auto const& wire = patch[w];
321  for (size_t d = 0; d < wire.size(); ++d, ++i) {
322  flat[i] = wire[d];
323  }
324  }
325 
326  return flat;
327 }
328 // ------------------------------------------------------
329 
330 bool nnet::PointIdAlg::isInsideFiducialRegion(unsigned int wire, float drift) const
331 {
332  size_t marginW = fPatchSizeW / 8; // fPatchSizeX/2 will make patch always completely filled
333  size_t marginD = fPatchSizeD / 8;
334 
335  size_t scaledDrift = (size_t)(drift / fDriftWindow);
336  if ((wire >= marginW) && (wire < fAlgView.fNWires - marginW) && (scaledDrift >= marginD) &&
337  (scaledDrift < fAlgView.fNScaledDrifts - marginD))
338  return true;
339  else
340  return false;
341 }
342 // ------------------------------------------------------
343 
344 // ------------------------------------------------------
345 // ------------------TrainingDataAlg---------------------
346 // ------------------------------------------------------
347 
349  : img::DataProviderAlg(config)
350  , fEdepTot(0)
351  , fWireProducerLabel(config.WireLabel())
352  , fHitProducerLabel(config.HitLabel())
353  , fTrackModuleLabel(config.TrackLabel())
354  , fSimulationProducerLabel(config.SimulationLabel())
355  , fSimChannelProducerLabel(config.SimChannelLabel())
356  , fSaveVtxFlags(config.SaveVtxFlags())
357  , fAdcDelay(config.AdcDelayTicks())
358  , fEventsPerBin(100, 0)
359 {
360  // If no sim channel producer is set then make it the same as the simulation label
362 
364 }
365 // ------------------------------------------------------
366 
368 // ------------------------------------------------------
369 
371  detinfo::DetectorClocksData const& clock_data,
372  detinfo::DetectorPropertiesData const& det_prop,
373  size_t wires,
374  size_t drifts)
375 {
376  auto view = img::DataProviderAlg::resizeView(clock_data, det_prop, wires, drifts);
377 
378  fWireDriftEdep.resize(wires);
379  for (auto& w : fWireDriftEdep) {
380  w.resize(view.fNCachedDrifts);
381  std::fill(w.begin(), w.end(), 0.0F);
382  }
383 
384  fWireDriftPdg.resize(wires);
385  for (auto& w : fWireDriftPdg) {
386  w.resize(view.fNCachedDrifts);
387  std::fill(w.begin(), w.end(), 0);
388  }
389  return view;
390 }
391 // ------------------------------------------------------
392 
393 bool nnet::TrainingDataAlg::setWireEdepsAndLabels(std::vector<float> const& edeps,
394  std::vector<int> const& pdgs,
395  size_t wireIdx)
396 {
397  if ((wireIdx >= fWireDriftEdep.size()) || (edeps.size() != pdgs.size())) { return false; }
398 
399  size_t dstep = 1;
400  if (fDownscaleFullView) { dstep = fDriftWindow; }
401 
402  if (edeps.size() / dstep > fAlgView.fNCachedDrifts) { return false; }
403 
404  auto& wEdep = fWireDriftEdep[wireIdx];
405  auto& wPdg = fWireDriftPdg[wireIdx];
406 
407  for (size_t i = 0; i < fAlgView.fNCachedDrifts; ++i) {
408  size_t i0 = i * dstep;
409  size_t i1 = (i + 1) * dstep;
410 
411  int best_pdg = pdgs[i0] & nnet::TrainingDataAlg::kPdgMask;
412  int vtx_flags = pdgs[i0] & nnet::TrainingDataAlg::kVtxMask;
413  int type_flags = pdgs[i0] & nnet::TrainingDataAlg::kTypeMask;
414  float max_edep = edeps[i0];
415  for (size_t k = i0 + 1; k < i1; ++k) {
416  float ek = edeps[k];
417  if (ek > max_edep) {
418  max_edep = ek;
419  best_pdg = pdgs[k] & nnet::TrainingDataAlg::kPdgMask; // remember best matching pdg
420  }
421  type_flags |= pdgs[k] & nnet::TrainingDataAlg::kTypeMask; // accumulate track type flags
422  vtx_flags |= pdgs[k] & nnet::TrainingDataAlg::kVtxMask; // accumulate all vtx flags
423  }
424 
425  wEdep[i] = max_edep;
426 
427  best_pdg |= type_flags;
428  if (fSaveVtxFlags) best_pdg |= vtx_flags;
429  wPdg[i] = best_pdg;
430  }
431 
432  return true;
433 }
434 // ------------------------------------------------------
435 
437  detinfo::DetectorClocksData const& clockData,
438  detinfo::DetectorPropertiesData const& detProp,
439  const TLorentzVector& tvec,
440  unsigned int plane) const
441 {
443  wd.Wire = 0;
444  wd.Drift = 0;
445  wd.TPC = -1;
446  wd.Cryo = -1;
447 
448  try {
449  geo::Point_t vtx{tvec.X(), tvec.Y(), tvec.Z()};
450  if (fGeometry->FindTPCAtPosition(vtx).isValid) {
451  geo::TPCID tpcid = fGeometry->FindTPCAtPosition(vtx);
452  unsigned int tpc = tpcid.TPC, cryo = tpcid.Cryostat;
453 
454  // correct for the time offset
455  float dx = tvec.T() * 1.e-3 * detProp.DriftVelocity();
456  int driftDir = fGeometry->TPC(tpcid).DetectDriftDirection();
457  if (driftDir == 1) { dx *= -1; }
458  else if (driftDir != -1) {
459  throw cet::exception("nnet::TrainingDataAlg") << "drift direction is not X." << std::endl;
460  }
461  vtx.SetX(tvec.X() + dx);
462 
463  wd.Wire = fGeometry->NearestWireID(vtx, geo::PlaneID{tpcid, plane}).Wire;
464  wd.Drift = detProp.ConvertXToTicks(vtx.X(), plane, tpc, cryo);
465  wd.TPC = tpc;
466  wd.Cryo = cryo;
467  }
468  }
469  catch (const geo::InvalidWireError& e) {
470  mf::LogWarning("TrainingDataAlg")
471  << "Vertex projection out of wire planes, just skipping this vertex.";
472  }
473  catch (...) {
474  mf::LogWarning("TrainingDataAlg") << "Vertex projection out of wire planes, skip MC vertex.";
475  }
476  return wd;
477 }
478 // ------------------------------------------------------
479 
481  const simb::MCParticle& particle,
482  const std::unordered_map<int, const simb::MCParticle*>& particleMap) const
483 {
484  const float minElectronLength2 = 2.5 * 2.5;
485  const float maxDeltaLength2 = 0.15 * 0.15;
486 
487  int pdg = abs(particle.PdgCode());
488  if (pdg != 11) return false; // should be applied only to electrons
489 
490  size_t nSec = particle.NumberDaughters();
491  for (size_t d = 0; d < nSec; ++d) {
492  auto d_search = particleMap.find(particle.Daughter(d));
493  if (d_search != particleMap.end()) {
494  auto const& daughter = *((*d_search).second);
495  int d_pdg = abs(daughter.PdgCode());
496  if (d_pdg != 22) { return false; } // not the end of the shower
497  }
498  }
499 
500  float trkLength2 = 0;
501  auto const* p = &particle;
502  bool branching = false;
503  while (!branching) {
504  trkLength2 += particleRange2(*p);
505  auto m_search = particleMap.find(p->Mother());
506  if (m_search != particleMap.end()) {
507  p = (*m_search).second;
508  int m_pdg = abs(p->PdgCode());
509  if (m_pdg == 11) {
510  nSec = p->NumberDaughters();
511  size_t ne = 0;
512  for (size_t d = 0; d < nSec; ++d) {
513  auto d_search = particleMap.find(p->Daughter(d));
514  if (d_search != particleMap.end()) {
515  auto const& daughter = *((*d_search).second);
516  int d_pdg = abs(daughter.PdgCode());
517  if (d_pdg == 11) {
518  if (particleRange2(daughter) > maxDeltaLength2) { ne++; }
519  }
520  }
521  }
522  if (ne > 1) { branching = true; }
523  }
524  else
525  break;
526  }
527  else
528  break;
529  }
530 
531  return (trkLength2 > minElectronLength2);
532 }
533 
535  const simb::MCParticle& particle,
536  const std::unordered_map<int, const simb::MCParticle*>& particleMap) const
537 {
538  bool hasElectron = false, hasNuMu = false, hasNuE = false;
539 
540  int pdg = abs(particle.PdgCode());
541  //if ((pdg == 13) && (particle.EndProcess() == "FastScintillation" || particle.EndProcess() == "Decay" || particle.EndProcess() == "muMinusCaptureAtRest")) // potential muon decay at rest
542  if ((pdg == 13) && (particle.EndProcess() == "FastScintillation" ||
543  particle.EndProcess() == "Decay")) // potential muon decay at rest
544  {
545  unsigned int nSec = particle.NumberDaughters();
546  for (size_t d = 0; d < nSec; ++d) {
547  auto d_search = particleMap.find(particle.Daughter(d));
548  if (d_search != particleMap.end()) {
549  auto const& daughter = *((*d_search).second);
550  int d_pdg = abs(daughter.PdgCode());
551  if (d_pdg == 11)
552  hasElectron = true;
553  else if (d_pdg == 14)
554  hasNuMu = true;
555  else if (d_pdg == 12)
556  hasNuE = true;
557  }
558  }
559  }
560 
561  return (hasElectron && hasNuMu && hasNuE);
562 }
563 
565  std::unordered_map<size_t, std::unordered_map<int, int>>& wireToDriftToVtxFlags,
566  detinfo::DetectorClocksData const& clockData,
567  detinfo::DetectorPropertiesData const& detProp,
568  const std::unordered_map<int, const simb::MCParticle*>& particleMap,
569  unsigned int plane) const
570 {
571  for (auto const& p : particleMap) {
572  auto const& particle = *p.second;
573 
574  double ekStart = 1000. * (particle.E() - particle.Mass());
575  double ekEnd = 1000. * (particle.EndE() - particle.Mass());
576  int pdg = abs(particle.PdgCode());
577  int flagsStart = nnet::TrainingDataAlg::kNone;
578  int flagsEnd = nnet::TrainingDataAlg::kNone;
579 
580  switch (pdg) {
581  case 22: // gamma
582  if ((particle.EndProcess() == "conv") && (ekStart > 40.0)) // conversion, gamma > 40MeV
583  {
584  flagsEnd = nnet::TrainingDataAlg::kConv;
585  }
586  break;
587 
588  case 11: // e+/-
589  if (isElectronEnd(particle, particleMap)) { flagsEnd = nnet::TrainingDataAlg::kElectronEnd; }
590  break;
591 
592  case 13: // mu+/-
593  if (isMuonDecaying(particle, particleMap)) { flagsEnd = nnet::TrainingDataAlg::kDecay; }
594  break;
595 
596  case 111: // pi0
597  flagsStart = nnet::TrainingDataAlg::kPi0;
598  break;
599 
600  case 321: // K+/-
601  case 211: // pi+/-
602  case 2212: // proton
603  if (ekStart > 50.0) {
604  if (particle.Mother() != 0) {
605  auto search = particleMap.find(particle.Mother());
606  if (search != particleMap.end()) {
607  auto const& mother = *((*search).second);
608  int m_pdg = abs(mother.PdgCode());
609  unsigned int nSec = mother.NumberDaughters();
610  unsigned int nVisible = 0;
611  if (nSec > 1) {
612  for (size_t d = 0; d < nSec; ++d) {
613  auto d_search = particleMap.find(mother.Daughter(d));
614  if (d_search != particleMap.end()) {
615  auto const& daughter = *((*d_search).second);
616  int d_pdg = abs(daughter.PdgCode());
617  if (((d_pdg == 2212) || (d_pdg == 211) || (d_pdg == 321)) &&
618  (1000. * (daughter.E() - daughter.Mass()) > 50.0)) {
619  ++nVisible;
620  }
621  }
622  }
623  }
624  // hadron with Ek > 50MeV (so well visible) and
625  // produced by another hadron (but not neutron, so not single track from nothing) or
626  // at least secondary hadrons with Ek > 50MeV (so this is a good kink or V-like)
627  if (((m_pdg != pdg) && (m_pdg != 2112)) || ((m_pdg != 2112) && (nVisible > 0)) ||
628  ((m_pdg == 2112) && (nVisible > 1))) {
629  flagsStart = nnet::TrainingDataAlg::kHadr;
630  }
631  }
632  }
633 
634  if (particle.EndProcess() == "FastScintillation") // potential decay at rest
635  {
636  unsigned int nSec = particle.NumberDaughters();
637  for (size_t d = 0; d < nSec; ++d) {
638  auto d_search = particleMap.find(particle.Daughter(d));
639  if (d_search != particleMap.end()) {
640  auto const& daughter = *((*d_search).second);
641  int d_pdg = abs(daughter.PdgCode());
642  if ((pdg == 321) && (d_pdg == 13)) {
644  break;
645  }
646  if ((pdg == 211) && (d_pdg == 13)) {
648  break;
649  }
650  }
651  }
652  }
653 
654  if ((particle.EndProcess() == "Decay") && (ekEnd > 200.0)) // decay in flight
655  {
656  unsigned int nSec = particle.NumberDaughters();
657  for (size_t d = 0; d < nSec; ++d) {
658  auto d_search = particleMap.find(particle.Daughter(d));
659  if (d_search != particleMap.end()) {
660  auto const& daughter = *((*d_search).second);
661  int d_pdg = abs(daughter.PdgCode());
662  if ((pdg == 321) && (d_pdg == 13)) {
663  flagsEnd = nnet::TrainingDataAlg::kHadr;
664  break;
665  }
666  if ((pdg == 211) && (d_pdg == 13)) {
667  flagsEnd = nnet::TrainingDataAlg::kHadr;
668  break;
669  }
670  }
671  }
672  }
673  }
674  break;
675 
676  default: continue;
677  }
678 
679  if (particle.Process() == "primary") { flagsStart |= nnet::TrainingDataAlg::kNuPri; }
680 
681  if (flagsStart != nnet::TrainingDataAlg::kNone) {
682  auto wd = getProjection(clockData, detProp, particle.Position(), plane);
683 
684  if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
685  wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= flagsStart;
686  }
687  }
688  if (flagsEnd != nnet::TrainingDataAlg::kNone) {
689  auto wd = getProjection(clockData, detProp, particle.EndPosition(), plane);
690  if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
691  wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= flagsEnd;
692  }
693  }
694 
695  //TY: check elastic/inelastic scattering
696  if (pdg == 321 || pdg == 211 || pdg == 2212) {
697  simb::MCTrajectory truetraj = particle.Trajectory();
698  auto thisTrajectoryProcessMap1 = truetraj.TrajectoryProcesses();
699  if (thisTrajectoryProcessMap1.size()) {
700  for (auto const& couple1 : thisTrajectoryProcessMap1) {
701  if ((truetraj.KeyToProcess(couple1.second)).find("Elastic") != std::string::npos) {
702  auto wd = getProjection(clockData, detProp, truetraj.at(couple1.first).first, plane);
703  if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
704  wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= nnet::TrainingDataAlg::kElastic;
705  }
706  }
707  if ((truetraj.KeyToProcess(couple1.second)).find("Inelastic") != std::string::npos) {
708  auto wd = getProjection(clockData, detProp, truetraj.at(couple1.first).first, plane);
709  if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
710  wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= nnet::TrainingDataAlg::kInelastic;
711  }
712  }
713  }
714  }
715  }
716  }
717 }
718 // ------------------------------------------------------
719 
721  detinfo::DetectorClocksData const& clockData,
722  detinfo::DetectorPropertiesData const& detProp,
723  unsigned int plane,
724  unsigned int tpc,
725  unsigned int cryo)
726 {
727 
729  std::vector<art::Ptr<recob::Wire>> Wirelist;
730 
731  if (event.getByLabel(fWireProducerLabel, wireHandle)) art::fill_ptr_vector(Wirelist, wireHandle);
732 
733  if (!setWireDriftData(clockData, detProp, *wireHandle, plane, tpc, cryo)) {
734  mf::LogError("TrainingDataAlg") << "Wire data not set.";
735  return false;
736  }
737 
738  // Hit info
740  std::vector<art::Ptr<recob::Hit>> Hitlist;
741 
742  if (event.getByLabel(fHitProducerLabel, HitHandle)) art::fill_ptr_vector(Hitlist, HitHandle);
743 
744  // Track info
746  std::vector<art::Ptr<recob::Track>> Tracklist;
747 
748  if (event.getByLabel(fTrackModuleLabel, TrackHandle))
749  art::fill_ptr_vector(Tracklist, TrackHandle);
750 
751  art::FindManyP<recob::Track> ass_trk_hits(HitHandle, event, fTrackModuleLabel);
752 
753  // Loop over wires (sorry about hard coded value) to fill in 1) pdg and 2) charge depo
754  for (size_t widx = 0; widx < 240; ++widx) {
755 
756  std::vector<float> labels_deposit(fAlgView.fNDrifts, 0); // full-drift-length buffers
757  std::vector<int> labels_pdg(fAlgView.fNDrifts, 0);
758 
759  // First, the charge depo
760  for (size_t subwidx = 0; subwidx < Wirelist.size(); ++subwidx) {
761  if (widx + 240 == Wirelist[subwidx]->Channel()) {
762  labels_deposit = Wirelist[subwidx]->Signal();
763  break;
764  }
765  }
766 
767  // Second, the pdg code
768  // This code finds the angle of the track and records
769  // events based on its angle to try to get an isometric sample
770  // instead of just a bunch of straight tracks
771 
772  // Meta code:
773  // For each hit:
774  // find farthest hit from point
775  // then find farthest hit from THAT one
776  // should be start and end of track, then just use trig
777 
778  for (size_t iHit = 0; iHit < Hitlist.size(); ++iHit) {
779 
780  if (Hitlist[iHit]->Channel() != widx + 240) { continue; }
781  if (Hitlist[iHit]->View() != 1) { continue; }
782 
783  // Make sure there is a track association
784  if (ass_trk_hits.at(iHit).size() == 0) { continue; }
785 
786  // Not sure about this
787  // Cutting on length to not just get a bunch of shower stubs
788  // Might add a lot of bias though
789  if (ass_trk_hits.at(iHit)[0]->Length() < 5) { continue; }
790 
791  // Search for farest hit from this one
792  int far_index = 0;
793  double far_dist = 0;
794 
795  for (size_t jHit = 0; jHit < Hitlist.size(); ++jHit) {
796  if (jHit == iHit) { continue; }
797  if (Hitlist[jHit]->View() != 1) { continue; }
798 
799  if (ass_trk_hits.at(jHit).size() == 0) { continue; }
800  if (ass_trk_hits.at(jHit)[0]->ID() != ass_trk_hits.at(iHit)[0]->ID()) { continue; }
801 
802  double dist = sqrt((Hitlist[iHit]->Channel() - Hitlist[jHit]->Channel()) *
803  (Hitlist[iHit]->Channel() - Hitlist[jHit]->Channel()) +
804  (Hitlist[iHit]->PeakTime() - Hitlist[jHit]->PeakTime()) *
805  (Hitlist[iHit]->PeakTime() - Hitlist[jHit]->PeakTime()));
806 
807  if (far_dist < dist) {
808  far_dist = dist;
809  far_index = jHit;
810  }
811  }
812 
813  // Search for the other end of the track
814  int other_end = 0;
815  int other_dist = 0;
816 
817  for (size_t jHit = 0; jHit < Hitlist.size(); ++jHit) {
818  if (jHit == iHit or int(jHit) == far_index) { continue; }
819  if (Hitlist[jHit]->View() != 1) { continue; }
820 
821  if (ass_trk_hits.at(jHit).size() == 0) { continue; }
822  if (ass_trk_hits.at(jHit)[0]->ID() != ass_trk_hits.at(iHit)[0]->ID()) { continue; }
823 
824  double dist = sqrt((Hitlist[far_index]->Channel() - Hitlist[jHit]->Channel()) *
825  (Hitlist[far_index]->Channel() - Hitlist[jHit]->Channel()) +
826  (Hitlist[far_index]->PeakTime() - Hitlist[jHit]->PeakTime()) *
827  (Hitlist[far_index]->PeakTime() - Hitlist[jHit]->PeakTime()));
828 
829  if (other_dist < dist) {
830  other_dist = dist;
831  other_end = jHit;
832  }
833  }
834 
835  // We have the end points now
836  double del_wire = double(Hitlist[other_end]->Channel() - Hitlist[far_index]->Channel());
837  double del_time = double(Hitlist[other_end]->PeakTime() - Hitlist[far_index]->PeakTime());
838  double hypo = sqrt(del_wire * del_wire + del_time * del_time);
839 
840  if (hypo == 0) { continue; } // Should never happen, but doing it anyway
841 
842  double cosser = TMath::Abs(del_wire / hypo);
843  double norm_ang = TMath::ACos(cosser) * 2 / TMath::Pi();
844 
845  // Using fEventsPerBin to keep track of number of hits per angle (normalized to 0 to 1)
846 
847  int binner = int(norm_ang * fEventsPerBin.size());
848  if (binner >= (int)fEventsPerBin.size()) {
849  binner = fEventsPerBin.size() - 1;
850  } // Dealing with rounding errors
851 
852  // So we should get a total of 5000 * 100 = 50,000 if we use the whole set
853  if (fEventsPerBin[binner] > 5000) { continue; }
854  fEventsPerBin[binner]++;
855 
856  // If survives everything, saves the pdg
857  labels_pdg[Hitlist[iHit]->PeakTime()] = 211; // Same as pion for now
858  }
859 
860  setWireEdepsAndLabels(labels_deposit, labels_pdg, widx);
861 
862  } // for each Wire
863 
864  return true;
865 }
866 
868  detinfo::DetectorClocksData const& clockData,
869  detinfo::DetectorPropertiesData const& detProp,
870  unsigned int plane,
871  unsigned int tpc,
872  unsigned int cryo)
873 {
875  event.getValidHandle<std::vector<recob::Wire>>(fWireProducerLabel);
876 
877  if (!setWireDriftData(clockData, detProp, *wireHandle, plane, tpc, cryo)) {
878  mf::LogError("TrainingDataAlg") << "Wire data not set.";
879  return false;
880  }
881 
882  if (!fSaveSimInfo || event.isRealData()) {
883  mf::LogInfo("TrainingDataAlg") << "Skip MC simulation info.";
884  return true;
885  }
886 
888  double electronsToGeV = 1. / larParameters->GeVToElectrons();
889 
890  auto particleHandle =
891  event.getValidHandle<std::vector<simb::MCParticle>>(fSimulationProducerLabel);
892 
893  auto simChannelHandle =
894  event.getValidHandle<std::vector<sim::SimChannel>>(fSimChannelProducerLabel);
895 
896  std::unordered_map<int, const simb::MCParticle*> particleMap;
897  for (auto const& particle : *particleHandle) {
898  particleMap[particle.TrackId()] = &particle;
899  }
900 
901  std::unordered_map<size_t, std::unordered_map<int, int>> wireToDriftToVtxFlags;
902  if (fSaveVtxFlags) collectVtxFlags(wireToDriftToVtxFlags, clockData, detProp, particleMap, plane);
903 
904  fEdepTot = 0;
905 
906  std::map<int, int> trackToPDG;
907  for (size_t widx = 0; widx < fAlgView.fNWires; ++widx) {
908  auto wireChannelNumber = fAlgView.fWireChannels[widx];
909  if (wireChannelNumber == raw::InvalidChannelID) continue;
910 
911  std::vector<float> labels_deposit(fAlgView.fNDrifts, 0); // full-drift-length buffers,
912  std::vector<int> labels_pdg(labels_deposit.size(), 0); // both of the same size,
913  int labels_size = labels_deposit.size(); // cached as int for comparisons below
914 
915  std::map<int, std::map<int, double>> timeToTrackToCharge;
916  for (auto const& channel : *simChannelHandle) {
917  if (channel.Channel() != wireChannelNumber) continue;
918 
919  auto const& timeSlices = channel.TDCIDEMap();
920  for (auto const& timeSlice : timeSlices) {
921  int time = timeSlice.first;
922 
923  auto const& energyDeposits = timeSlice.second;
924  for (auto const& energyDeposit : energyDeposits) {
925  int pdg = 0;
926  int tid = energyDeposit.trackID;
927  if (tid < 0) // negative tid means it is EM activity, and -tid is the mother
928  {
929  pdg = 11;
930  tid = -tid;
931 
932  auto search = particleMap.find(tid);
933  if (search == particleMap.end()) {
934  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
935  continue;
936  }
937  auto const& mother = *((*search).second); // mother particle of this EM
938  int mPdg = abs(mother.PdgCode());
939  if ((mPdg == 13) || (mPdg == 211) || (mPdg == 2212)) {
940  if (energyDeposit.numElectrons > 10)
941  pdg |= nnet::TrainingDataAlg::kDelta; // tag delta ray
942  }
943  }
944  else {
945  auto search = particleMap.find(tid);
946  if (search == particleMap.end()) {
947  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
948  continue;
949  }
950  auto const& particle = *((*search).second);
951  pdg = abs(particle.PdgCode());
952 
953  if (particle.Process() == "primary") {
954  if (pdg == 11) {
955  pdg |= nnet::TrainingDataAlg::kPriEl; // tag primary
956  }
957  else if (pdg == 13) {
958  pdg |= nnet::TrainingDataAlg::kPriMu; // tag primary
959  }
960  }
961 
962  /*
963  auto msearch = particleMap.find(particle.Mother());
964  if (msearch != particleMap.end()) {
965  auto const& mother = *((*msearch).second);
966  if (pdg == 11) // electron, check if it is Michel
967  {
968  if (nnet::TrainingDataAlg::isMuonDecaying(mother, particleMap)) {
969  std::cout<<particle.Process()<<std::endl;
970  pdg |= nnet::TrainingDataAlg::kMichel; // tag Michel
971  }
972  }
973  }
974  */
975  if (pdg == 11) { // electron, check if it is Michel or delta ray
976  if (particle.Process() == "Decay") {
977  pdg |= nnet::TrainingDataAlg::kMichel; // tag Michel
978  }
979  else if (particle.Process() == "muIoni") {
980  pdg |= nnet::TrainingDataAlg::kDelta; // tag delta ray
981  }
982  }
983  }
984 
985  trackToPDG[energyDeposit.trackID] = pdg;
986 
987  double energy = energyDeposit.numElectrons * electronsToGeV;
988  timeToTrackToCharge[time][energyDeposit.trackID] += energy;
989  fEdepTot += energy;
990 
991  } // loop over energy deposits
992  } // loop over time slices
993  } // for each SimChannel
994 
996  for (auto const& ttc : timeToTrackToCharge) {
997  float max_deposit = 0.0;
998  int max_pdg = 0;
999  for (auto const& tc : ttc.second) {
1000 
1001  if (tc.second > max_deposit) {
1002  max_deposit = tc.second;
1003  max_pdg = trackToPDG[tc.first];
1004  }
1005  }
1006 
1007  int tick_idx = clockData.TPCTDC2Tick(ttc.first) + fAdcDelay;
1008 
1009  if (tick_idx < labels_size && tick_idx >= 0) {
1010  labels_deposit[tick_idx] = max_deposit;
1011  labels_pdg[tick_idx] = max_pdg & type_pdg_mask;
1012  }
1013  }
1014 
1015  for (auto const& drift_flags : wireToDriftToVtxFlags[widx]) {
1016  int drift = drift_flags.first, flags = drift_flags.second;
1017  if ((drift >= 0) && (drift < labels_size)) { labels_pdg[drift] |= flags; }
1018  }
1019  setWireEdepsAndLabels(labels_deposit, labels_pdg, widx);
1020  } // for each Wire
1021 
1022  return true;
1023 }
1024 // ------------------------------------------------------
1025 
1026 bool nnet::TrainingDataAlg::findCrop(float max_e_cut,
1027  unsigned int& w0,
1028  unsigned int& w1,
1029  unsigned int& d0,
1030  unsigned int& d1) const
1031 {
1032  if (fWireDriftEdep.empty() || fWireDriftEdep.front().empty()) return false;
1033 
1034  float max_cut = 0.25 * max_e_cut;
1035 
1036  w0 = 0;
1037  float cut = 0;
1038  while (w0 < fWireDriftEdep.size()) {
1039  for (auto const d : fWireDriftEdep[w0])
1040  cut += d;
1041  if (cut < max_cut)
1042  w0++;
1043  else
1044  break;
1045  }
1046  w1 = fWireDriftEdep.size() - 1;
1047  cut = 0;
1048  while (w1 > w0) {
1049  for (auto const d : fWireDriftEdep[w1])
1050  cut += d;
1051  if (cut < max_cut)
1052  w1--;
1053  else
1054  break;
1055  }
1056  w1++;
1057 
1058  d0 = 0;
1059  cut = 0;
1060  while (d0 < fWireDriftEdep.front().size()) {
1061  for (size_t i = w0; i < w1; ++i)
1062  cut += fWireDriftEdep[i][d0];
1063  if (cut < max_cut)
1064  d0++;
1065  else
1066  break;
1067  }
1068  d1 = fWireDriftEdep.front().size() - 1;
1069  cut = 0;
1070  while (d1 > d0) {
1071  for (size_t i = w0; i < w1; ++i)
1072  cut += fWireDriftEdep[i][d1];
1073  if (cut < max_cut)
1074  d1--;
1075  else
1076  break;
1077  }
1078  d1++;
1079 
1080  unsigned int margin = 20;
1081  if ((w1 - w0 > 8) && (d1 - d0 > 8)) {
1082  if (w0 < margin)
1083  w0 = 0;
1084  else
1085  w0 -= margin;
1086 
1087  if (w1 > fWireDriftEdep.size() - margin)
1088  w1 = fWireDriftEdep.size();
1089  else
1090  w1 += margin;
1091 
1092  if (d0 < margin)
1093  d0 = 0;
1094  else
1095  d0 -= margin;
1096 
1097  if (d1 > fWireDriftEdep.front().size() - margin)
1098  d1 = fWireDriftEdep.front().size();
1099  else
1100  d1 += margin;
1101 
1102  return true;
1103  }
1104  else
1105  return false;
1106 }
1107 // ------------------------------------------------------
bool isInsideFiducialRegion(unsigned int wire, float drift) const
Definition: PointIdAlg.cxx:330
void collectVtxFlags(std::unordered_map< size_t, std::unordered_map< int, int >> &wireToDriftToVtxFlags, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::unordered_map< int, const simb::MCParticle * > &particleMap, unsigned int plane) const
Definition: PointIdAlg.cxx:564
TRandom r
Definition: spectrum.C:23
Store parameters for running LArG4.
geo::GeometryCore const * fGeometry
static std::unique_ptr< Graph > create(const char *graph_file_name, const std::vector< std::string > &outputs={}, bool use_bundle=false, int ninputs=1, int noutputs=1)
Definition: tf_graph.h:29
int PdgCode() const
Definition: MCParticle.h:213
bool setWireEdepsAndLabels(std::vector< float > const &edeps, std::vector< int > const &pdgs, size_t wireIdx)
Definition: PointIdAlg.cxx:393
Utilities related to art service access.
bool isMuonDecaying(const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const
Definition: PointIdAlg.cxx:534
std::vector< std::vector< float > > Run(std::vector< std::vector< std::vector< float >>> const &inps, int samples=-1) override
Definition: PointIdAlg.cxx:118
std::vector< std::vector< float > > fWireDriftPatch
Definition: PointIdAlg.h:154
bool setDataEventData(const art::Event &event, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: PointIdAlg.cxx:720
unsigned int fAdcDelay
Definition: PointIdAlg.h:346
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
art::InputTag fTrackModuleLabel
Definition: PointIdAlg.h:340
Declaration of signal hit object.
virtual DataProviderAlgView resizeView(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, size_t wires, size_t drifts)
const value_type & at(const size_type i) const
Definition: MCTrajectory.h:175
bool setEventData(const art::Event &event, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: PointIdAlg.cxx:867
virtual void set_data(std::vector< std::vector< std::vector< float >>> const &)
Definition: keras_model.h:53
bool bufferPatch(size_t wire, float drift, std::vector< std::vector< float >> &patch) const
Definition: PointIdAlg.h:158
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:210
std::string KeyToProcess(unsigned char const &key) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
DataProviderAlg(const fhicl::ParameterSet &pset)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.cc:53
WireDrift getProjection(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TLorentzVector &tvec, unsigned int plane) const
Definition: PointIdAlg.cxx:436
double TPCTDC2Tick(double const tdc) const
Given electronics clock count [tdc] returns TPC time-tick.
std::vector< float > compute_output(keras::DataChunk *dc)
Definition: keras_model.cc:421
std::vector< std::string > fNNetOutputs
Definition: PointIdAlg.h:151
~PointIdAlg() override
Definition: PointIdAlg.cxx:204
bool isCurrentPatch(unsigned int wire, float drift) const
Definition: PointIdAlg.cxx:293
int NumberDaughters() const
Definition: MCParticle.h:218
int Daughter(const int i) const
Definition: MCParticle.cxx:118
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
art::InputTag fWireProducerLabel
Definition: PointIdAlg.h:338
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
bool isElectronEnd(const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const
Definition: PointIdAlg.cxx:480
std::string const & label() const noexcept
Definition: InputTag.cc:79
std::vector< float > predictIdVector(unsigned int wire, float drift) const
calculate multi-class probabilities for [wire, drift] point
Definition: PointIdAlg.cxx:239
Collection of exceptions for Geometry system.
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:31
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
std::string EndProcess() const
Definition: MCParticle.h:217
fhicl::Sequence< std::string > NNetOutputs
Definition: PointIdAlg.h:106
std::vector< raw::ChannelID_t > fWireChannels
Definition: EmTrack.h:40
double energy
Definition: plottest35.C:25
virtual std::vector< float > Run(std::vector< std::vector< float >> const &inp2d)=0
std::vector< size_t > fEventsPerBin
Definition: PointIdAlg.h:348
Float_t d
Definition: plot.C:235
double ConvertXToTicks(double X, int p, int t, int c) const
KerasModelInterface(const char *modelFileName)
Definition: PointIdAlg.cxx:85
Provides recob::Track data product.
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
art::InputTag fSimChannelProducerLabel
Definition: PointIdAlg.h:342
bool findCrop(float max_e_cut, unsigned int &w0, unsigned int &w1, unsigned int &d0, unsigned int &d1) const
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
size_t fPatchSizeD
Definition: PointIdAlg.h:155
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
std::vector< float > Run(std::vector< std::vector< float >> const &inp2d) override
Definition: PointIdAlg.cxx:92
ProcessMap const & TrajectoryProcesses() const
Definition: MCTrajectory.h:188
Definition of data types for geometry description.
keras::KerasModel m
Definition: PointIdAlg.h:81
DataProviderAlgView fAlgView
std::vector< std::vector< int > > fWireDriftPdg
Definition: PointIdAlg.h:336
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
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)
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
fhicl::Atom< std::string > NNetModelFile
Definition: PointIdAlg.h:104
~TrainingDataAlg() override
TfModelInterface(const char *modelFileName)
Definition: PointIdAlg.cxx:108
std::vector< std::vector< float > > fWireDriftEdep
Definition: PointIdAlg.h:335
Contains all timing reference information for the detector.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
std::string findFile(const char *fileName) const
Definition: PointIdAlg.cxx:67
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
Definition: TPCGeo.cxx:149
unsigned int Cryo() const
Pool sum of pixels in a patch around the wire/drift pixel.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
static float particleRange2(const simb::MCParticle &particle)
Definition: PointIdAlg.h:321
img::DataProviderAlgView resizeView(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, size_t wires, size_t drifts) override
Definition: PointIdAlg.cxx:370
object containing MC truth information necessary for making RawDigits and doing back tracking ...
size_t fCurrentScaledDrift
Definition: PointIdAlg.h:157
Exception thrown on invalid wire number.
Definition: Exceptions.h:39
Declaration of basic channel signal object.
nnet::ModelInterface * fNNet
Definition: PointIdAlg.h:152
size_t fPatchSizeW
Definition: PointIdAlg.h:155
art::InputTag fSimulationProducerLabel
Definition: PointIdAlg.h:341
std::string fNNetModelFilePath
Definition: PointIdAlg.h:150
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
bool isSamePatch(unsigned int wire1, float drift1, unsigned int wire2, float drift2) const
test if two wire/drift coordinates point to the same patch
Definition: PointIdAlg.cxx:276
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
std::vector< std::vector< float > > predictIdVectors(std::vector< std::pair< unsigned int, float >> points) const
Definition: PointIdAlg.cxx:257
art::InputTag fHitProducerLabel
Definition: PointIdAlg.h:339
double GeVToElectrons() const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Event finding and building.
TrainingDataAlg(const fhicl::ParameterSet &pset)
Definition: PointIdAlg.h:252
size_t fCurrentWireIdx
Definition: PointIdAlg.h:157
unsigned int TPC() const
float predictIdValue(unsigned int wire, float drift, size_t outIdx=0) const
calculate single-value prediction (2-class probability) for [wire, drift] point
Definition: PointIdAlg.cxx:218
PointIdAlg(const fhicl::ParameterSet &pset)
Definition: PointIdAlg.h:114
static std::vector< float > flattenData2D(std::vector< std::vector< float >> const &patch)
Definition: PointIdAlg.cxx:309