LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
reco_tool::CandHitMorphological Class Reference
Inheritance diagram for reco_tool::CandHitMorphological:
reco_tool::ICandidateHitFinder

Public Member Functions

 CandHitMorphological (const fhicl::ParameterSet &pset)
 
void findHitCandidates (const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
 
void MergeHitCandidates (const recob::Wire::RegionsOfInterest_t::datarange_t &, const HitCandidateVec &, MergeHitCandidateVec &) const override
 

Private Types

using MaxMinPair = std::pair< Waveform::const_iterator, Waveform::const_iterator >
 
using CandHitParams = std::tuple< Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator >
 
using CandHitParamsVec = std::vector< CandHitParams >
 
using HitCandidateVec = std::vector< HitCandidate >
 
using MergeHitCandidateVec = std::vector< HitCandidateVec >
 
using Waveform = std::vector< float >
 

Private Member Functions

void findHitCandidates (Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator, const size_t, float, HitCandidateVec &) const
 
void findHitCandidates (Waveform::const_iterator, Waveform::const_iterator, const size_t, int, float, HitCandidateVec &) const
 
bool getListOfHitCandidates (Waveform::const_iterator, Waveform::const_iterator, int, float, CandHitParamsVec &) const
 
Waveform::const_iterator findNearestMax (Waveform::const_iterator, Waveform::const_iterator) const
 
Waveform::const_iterator findNearestMin (Waveform::const_iterator, Waveform::const_iterator) const
 
Waveform::const_iterator findStartTick (Waveform::const_iterator, Waveform::const_iterator) const
 
Waveform::const_iterator findStopTick (Waveform::const_iterator, Waveform::const_iterator) const
 

Private Attributes

const size_t fPlane
 
const float fDilationThreshold
 
const float fDilationFraction
 
const float fErosionFraction
 
const int fMinDeltaTicks
 
const float fMinDeltaPeaks
 
const float fMinHitHeight
 
const size_t fNumInterveningTicks
 
const int fStructuringElement
 
const bool fOutputHistograms
 
const bool fOutputWaveforms
 
const float fFitNSigmaFromCenter
 
art::TFileDirectory * fHistDirectory
 
TH1F * fDStopStartHist
 
TH1F * fDMaxTickMinTickHist
 
TH1F * fDMaxDerivMinDerivHist
 
TH1F * fMaxErosionHist
 
TH1F * fMaxDilationHist
 
TH1F * fMaxDilEroRatHist
 
size_t fLastChannel
 
size_t fChannelCnt
 
std::unique_ptr< reco_tool::IWaveformToolfWaveformTool
 
const geo::GeometryCorefGeometry = lar::providerFrom<geo::Geometry>()
 

Detailed Description

Definition at line 25 of file CandHitMorphological_tool.cc.

Member Typedef Documentation

Definition at line 66 of file CandHitMorphological_tool.cc.

Definition at line 37 of file ICandidateHitFinder.h.

Definition at line 38 of file ICandidateHitFinder.h.

using reco_tool::ICandidateHitFinder::Waveform = std::vector<float>
inherited

Definition at line 40 of file ICandidateHitFinder.h.

Constructor & Destructor Documentation

reco_tool::CandHitMorphological::CandHitMorphological ( const fhicl::ParameterSet pset)
explicit

Definition at line 126 of file CandHitMorphological_tool.cc.

References art::errors::Configuration, dir, fChannelCnt, fDMaxDerivMinDerivHist, fDMaxTickMinTickHist, fDStopStartHist, fHistDirectory, fLastChannel, fMaxDilationHist, fMaxDilEroRatHist, fMaxErosionHist, fOutputHistograms, fOutputWaveforms, fPlane, fWaveformTool, art::ServiceHandle< T, SCOPE >::get(), fhicl::ParameterSet::get(), and art::Globals::instance().

127  : fPlane(pset.get<size_t>("Plane", 0))
128  , fDilationThreshold(pset.get<float>("DilationThreshold", 4.))
129  , fDilationFraction(pset.get<float>("DilationFraction", 0.75))
130  , fErosionFraction(pset.get<float>("ErosionFraction", 0.2))
131  , fMinDeltaTicks(pset.get<int>("MinDeltaTicks", 0))
132  , fMinDeltaPeaks(pset.get<float>("MinDeltaPeaks", 0.025))
133  , fMinHitHeight(pset.get<float>("MinHitHeight", 1.0))
134  , fNumInterveningTicks(pset.get<size_t>("NumInterveningTicks", 6))
135  , fStructuringElement(pset.get<int>("StructuringElement", 20))
136  , fOutputHistograms(pset.get<bool>("OutputHistograms", false))
137  , fOutputWaveforms(pset.get<bool>("OutputWaveforms", false))
138  , fFitNSigmaFromCenter(pset.get<float>("FitNSigmaFromCenter", 5.))
139  {
140 
141  if (art::Globals::instance()->nthreads() > 1u) {
142  if (fOutputHistograms) {
144  << "Cannot fill histograms when multiple threads configured, please set "
145  "fOutputHistograms to false or change number of threads to 1\n";
146  }
147 
148  if (fOutputWaveforms) {
150  << "Cannot write output waveforms when multiple threads configured, please set "
151  "fOutputHistograms to false or change number of threads to 1\n";
152  }
153  }
154  // Recover the baseline tool
155  fWaveformTool =
156  art::make_tool<reco_tool::IWaveformTool>(pset.get<fhicl::ParameterSet>("WaveformAlgs"));
157 
158  // Set the last channel to some nonsensical value
159  fLastChannel = std::numeric_limits<size_t>::max();
160  fChannelCnt = 0;
161 
162  // If asked, define the global histograms
163  if (fOutputHistograms) {
164  // Access ART's TFileService, which will handle creating and writing
165  // histograms and n-tuples for us.
167 
168  fHistDirectory = tfs.get();
169 
170  // Make a directory for these histograms
171  art::TFileDirectory dir = fHistDirectory->mkdir(Form("HitPlane_%1zu", fPlane));
172 
174  dir.make<TH1F>(Form("DStopStart_%1zu", fPlane), ";Delta Stop/Start;", 100, 0., 100.);
176  dir.make<TH1F>(Form("DMaxTMinT_%1zu", fPlane), ";Delta Max/Min Tick;", 100, 0., 100.);
178  dir.make<TH1F>(Form("DMaxDMinD_%1zu", fPlane), ";Delta Max/Min Deriv;", 200, 0., 100.);
180  dir.make<TH1F>(Form("MaxErosion_%1zu", fPlane), ";Max Erosion;", 200, -50., 150.);
182  dir.make<TH1F>(Form("MaxDilation_%1zu", fPlane), ";Max Dilation;", 200, -50., 150.);
184  dir.make<TH1F>(Form("MaxDilEroRat_%1zu", fPlane), ";Max Dil/Ero;", 200, -1., 1.);
185  }
186 
187  return;
188  }
T * get() const
Definition: ServiceHandle.h:69
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
T get(std::string const &key) const
Definition: ParameterSet.h:314
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
TDirectory * dir
Definition: macro.C:5
static Globals * instance()
Definition: Globals.cc:17

Member Function Documentation

void reco_tool::CandHitMorphological::findHitCandidates ( const recob::Wire::RegionsOfInterest_t::datarange_t &  dataRange,
const size_t  roiStartTick,
const size_t  channel,
const size_t  eventCount,
HitCandidateVec hitCandidateVec 
) const
overridevirtual

Implements reco_tool::ICandidateHitFinder.

Definition at line 190 of file CandHitMorphological_tool.cc.

References util::abs(), geo::GeometryCore::ChannelToWire(), dir, fChannelCnt, fDilationThreshold, fDMaxDerivMinDerivHist, fDMaxTickMinTickHist, fDStopStartHist, fFitNSigmaFromCenter, fGeometry, fHistDirectory, fLastChannel, fMaxDilationHist, fMaxDilEroRatHist, fMaxErosionHist, fOutputHistograms, fOutputWaveforms, fStructuringElement, and fWaveformTool.

Referenced by findHitCandidates().

196  {
197  // In this case we want to find hit candidates based on the derivative of of the input waveform
198  // We get this from our waveform algs too...
199  Waveform rawDerivativeVec;
200  Waveform derivativeVec;
201 
202  // Recover the actual waveform
203  const Waveform& waveform = dataRange.data();
204 
205  fWaveformTool->firstDerivative(waveform, rawDerivativeVec);
206  fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
207 
208  // Now we get the erosion/dilation vectors
209  Waveform erosionVec;
210  Waveform dilationVec;
211  Waveform averageVec;
212  Waveform differenceVec;
213 
214  reco_tool::HistogramMap histogramMap;
215 
216  // Compute the morphological filter vectors
217  fWaveformTool->getErosionDilationAverageDifference(waveform,
219  histogramMap,
220  erosionVec,
221  dilationVec,
222  averageVec,
223  differenceVec);
224 
225  // Now find the hits
226  findHitCandidates(derivativeVec.begin(),
227  derivativeVec.end(),
228  erosionVec.begin(),
229  erosionVec.end(),
230  dilationVec.begin(),
231  dilationVec.end(),
232  roiStartTick,
234  hitCandidateVec);
235 
236  // Limit start and stop tick to the neighborhood of the peak
237  if (fFitNSigmaFromCenter > 0) {
238  for (auto& hc : hitCandidateVec) {
239  auto startCand = hc.hitCenter - fFitNSigmaFromCenter * hc.hitSigma;
240  if (startCand >= 0) hc.startTick = std::max(hc.startTick, size_t(startCand));
241  hc.stopTick =
242  std::min(hc.stopTick, size_t(hc.hitCenter + fFitNSigmaFromCenter * hc.hitSigma));
243  }
244  }
245 
246  // Reset the hit height from the input waveform
247  for (auto& hitCandidate : hitCandidateVec) {
248  size_t centerIdx = hitCandidate.hitCenter;
249 
250  hitCandidate.hitHeight = waveform.at(centerIdx);
251  }
252 
253  // Keep track of histograms if requested
254  if (fOutputWaveforms) {
255  // Recover the details...
256  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
257  size_t plane = wids[0].Plane;
258  size_t cryo = wids[0].Cryostat;
259  size_t tpc = wids[0].TPC;
260  size_t wire = wids[0].Wire;
261 
262  if (channel != fLastChannel) fChannelCnt = 0;
263 
264  // Make a directory for these histograms
265  art::TFileDirectory dir = fHistDirectory->mkdir(
266  Form("Event%04zu/c%1zuT%1zuP%1zu/Wire_%05zu", eventCount, cryo, tpc, plane, wire));
267 
268  size_t waveformSize = waveform.size();
269  size_t waveStart = dataRange.begin_index();
270 
271  TProfile* waveHist =
272  dir.make<TProfile>(Form("HWfm_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
273  "Waveform",
274  waveformSize,
275  0,
276  waveformSize,
277  -500.,
278  500.);
279  TProfile* derivHist =
280  dir.make<TProfile>(Form("HDer_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
281  "Derivative",
282  waveformSize,
283  0,
284  waveformSize,
285  -500.,
286  500.);
287  TProfile* erosionHist =
288  dir.make<TProfile>(Form("HEro_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
289  "Erosion",
290  waveformSize,
291  0,
292  waveformSize,
293  -500.,
294  500.);
295  TProfile* dilationHist =
296  dir.make<TProfile>(Form("HDil_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
297  "Dilation",
298  waveformSize,
299  0,
300  waveformSize,
301  -500.,
302  500.);
303  TProfile* candHitHist =
304  dir.make<TProfile>(Form("HCan_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
305  "Cand Hits",
306  waveformSize,
307  0,
308  waveformSize,
309  -500.,
310  500.);
311  TProfile* maxDerivHist =
312  dir.make<TProfile>(Form("HMax_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
313  "Maxima",
314  waveformSize,
315  0,
316  waveformSize,
317  -500.,
318  500.);
319  TProfile* strtStopHist =
320  dir.make<TProfile>(Form("HSSS_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
321  "Start/Stop",
322  waveformSize,
323  0,
324  waveformSize,
325  -500.,
326  500.);
327 
328  // Fill wave/derivative
329  for (size_t idx = 0; idx < waveform.size(); idx++) {
330  waveHist->Fill(roiStartTick + idx, waveform.at(idx));
331  derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
332  erosionHist->Fill(roiStartTick + idx, erosionVec.at(idx));
333  dilationHist->Fill(roiStartTick + idx, dilationVec.at(idx));
334  }
335 
336  // Fill hits
337  for (const auto& hitCandidate : hitCandidateVec) {
338  candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
339  maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
340  maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
341  strtStopHist->Fill(hitCandidate.startTick, waveform.at(hitCandidate.startTick));
342  strtStopHist->Fill(hitCandidate.stopTick, waveform.at(hitCandidate.stopTick));
343  }
344 
345  fLastChannel = channel;
346  fChannelCnt++;
347  }
348 
349  if (fOutputHistograms) {
350  // Fill hits
351  for (const auto& hitCandidate : hitCandidateVec) {
352  fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
353  fDMaxTickMinTickHist->Fill(hitCandidate.minTick - hitCandidate.maxTick, 1.);
354  fDMaxDerivMinDerivHist->Fill(hitCandidate.maxDerivative - hitCandidate.minDerivative, 1.);
355  }
356 
357  // Get the max dilation/erosion
358  Waveform::const_iterator maxDilationItr =
359  std::max_element(dilationVec.begin(), dilationVec.end());
360  Waveform::const_iterator maxErosionItr =
361  std::max_element(erosionVec.begin(), erosionVec.end());
362 
363  float dilEroRat(1.);
364 
365  if (std::abs(*maxDilationItr) > 0.) dilEroRat = *maxErosionItr / *maxDilationItr;
366 
367  fMaxErosionHist->Fill(*maxErosionItr, 1.);
368  fMaxDilationHist->Fill(*maxDilationItr, 1.);
369  fMaxDilEroRatHist->Fill(dilEroRat, 1.);
370  }
371 
372  return;
373  }
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
constexpr auto abs(T v)
Returns the absolute value of the argument.
intermediate_table::const_iterator const_iterator
void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
std::map< int, TProfile * > HistogramMap
Definition: IWaveformTool.h:42
TDirectory * dir
Definition: macro.C:5
void reco_tool::CandHitMorphological::findHitCandidates ( Waveform::const_iterator  derivStartItr,
Waveform::const_iterator  derivStopItr,
Waveform::const_iterator  erosionStartItr,
Waveform::const_iterator  erosionStopItr,
Waveform::const_iterator  dilationStartItr,
Waveform::const_iterator  dilationStopItr,
const size_t  roiStartTick,
float  dilationThreshold,
HitCandidateVec hitCandidateVec 
) const
private

Definition at line 375 of file CandHitMorphological_tool.cc.

References fDilationThreshold, fErosionFraction, findHitCandidates(), fMinDeltaPeaks, and fMinDeltaTicks.

384  {
385  // This function aims to use the erosion/dilation vectors to find candidate hit regions
386  // Once armed with a region then the "standard" differential approach is used to return the candidate peaks
387 
388  // Don't do anything if not enough ticks
389  int ticksInInputWaveform = std::distance(derivStartItr, derivStopItr);
390 
391  if (ticksInInputWaveform < fMinDeltaTicks) return;
392 
393  // This function is recursive, we start by finding the largest element of the dilation vector
394  Waveform::const_iterator maxItr = std::max_element(dilationStartItr, dilationStopItr);
395  float maxVal = *maxItr;
396 
397  // Check that the peak is of reasonable height...
398  if (maxVal < dilationThreshold) return;
399 
400  int maxBin = std::distance(dilationStartItr, maxItr);
401 
402  // The candidate hit region we want lies between the two nearest minima to the maximum we just found
403  // subject to the condition that the erosion vector has return to less than zero
404  Waveform::const_iterator firstDerItr = derivStartItr + maxBin;
405  Waveform::const_iterator erosionItr = erosionStartItr + maxBin;
406 
407  float firstDerivValue = -1.;
408  float erosionCutValue = fErosionFraction * maxVal;
409 
410  // Search for starting point
411  while (firstDerItr != derivStartItr) {
412  // Look for the turnover point
413  if (*erosionItr-- < erosionCutValue) {
414  // We are looking for the zero crossing signifying a minimum value in the waveform
415  // (so the previous derivative < 0 while current is > 0)
416  // We are moving "backwards" so the current value <= 0, the previous value > 0
417  if (*firstDerItr * firstDerivValue <= 0. && firstDerivValue > 0.) break;
418  }
419 
420  firstDerivValue = *firstDerItr--;
421  }
422 
423  // Set the start bin
424  int hitRegionStart = std::distance(derivStartItr, firstDerItr);
425 
426  // Now go the other way
427  Waveform::const_iterator lastDerItr = derivStartItr + maxBin;
428 
429  // Reset the local variables
430  float lastDerivValue = 1.;
431 
432  erosionItr = erosionStartItr + maxBin;
433 
434  // Search for starting point
435  while (lastDerItr != derivStopItr) {
436  if (*erosionItr++ <= erosionCutValue) {
437  // We are again looking for the zero crossing signifying a minimum value in the waveform
438  // This time we are moving forward, so test is that previous value < 0, new value >= 0
439  if (*lastDerItr * lastDerivValue <= 0. && lastDerivValue < 0.) break;
440  }
441 
442  lastDerivValue = *lastDerItr++;
443  }
444 
445  // Set the stop bin
446  int hitRegionStop = std::distance(derivStartItr, lastDerItr);
447 
448  // Recursive call to find any hits in front of where we are now
449  if (hitRegionStart > fMinDeltaTicks)
450  findHitCandidates(derivStartItr,
451  derivStartItr + hitRegionStart,
452  erosionStartItr,
453  erosionStartItr + hitRegionStart,
454  dilationStartItr,
455  dilationStartItr + hitRegionStart,
456  roiStartTick,
458  hitCandidateVec);
459 
460  // Call the differential hit finding to get the actual hits within the region
461  findHitCandidates(derivStartItr + hitRegionStart,
462  derivStartItr + hitRegionStop,
463  roiStartTick + hitRegionStart,
466  hitCandidateVec);
467 
468  // Now call ourselves again to find any hits trailing the region we just identified
469  if (std::distance(lastDerItr, derivStopItr) > fMinDeltaTicks)
470  findHitCandidates(derivStartItr + hitRegionStop,
471  derivStopItr,
472  erosionStartItr + hitRegionStop,
473  erosionStopItr,
474  dilationStartItr + hitRegionStop,
475  dilationStopItr,
476  roiStartTick + hitRegionStop,
478  hitCandidateVec);
479 
480  return;
481  }
intermediate_table::const_iterator const_iterator
void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
void reco_tool::CandHitMorphological::findHitCandidates ( Waveform::const_iterator  startItr,
Waveform::const_iterator  stopItr,
const size_t  roiStartTick,
int  dTicksThreshold,
float  dPeakThreshold,
HitCandidateVec hitCandidateVec 
) const
private

Definition at line 483 of file CandHitMorphological_tool.cc.

References getListOfHitCandidates(), reco_tool::ICandidateHitFinder::HitCandidate::hitCenter, reco_tool::ICandidateHitFinder::HitCandidate::hitHeight, reco_tool::ICandidateHitFinder::HitCandidate::hitSigma, art::left(), reco_tool::ICandidateHitFinder::HitCandidate::maxDerivative, reco_tool::ICandidateHitFinder::HitCandidate::maxTick, reco_tool::ICandidateHitFinder::HitCandidate::minDerivative, reco_tool::ICandidateHitFinder::HitCandidate::minTick, art::right(), reco_tool::ICandidateHitFinder::HitCandidate::startTick, and reco_tool::ICandidateHitFinder::HitCandidate::stopTick.

489  {
490  // Search for candidate hits...
491  // Strategy is to get the list of all possible max/min pairs of the input derivative vector and then
492  // look for candidate hits in that list
493  CandHitParamsVec candHitParamsVec;
494 
496  startItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec)) {
497  // We've been given a list of candidate hits so now convert to hits
498  // Version one... simply convert all the candidates
499  for (const auto& tuple : candHitParamsVec) {
500  // Create a new hit candidate and store away
501  HitCandidate hitCandidate;
502 
503  Waveform::const_iterator candStartItr = std::get<0>(tuple);
504  Waveform::const_iterator maxItr = std::get<1>(tuple);
505  Waveform::const_iterator minItr = std::get<2>(tuple);
506  Waveform::const_iterator candStopItr = std::get<3>(tuple);
507 
508  Waveform::const_iterator peakItr =
509  std::min_element(maxItr, minItr, [](const auto& left, const auto& right) {
510  return std::fabs(left) < std::fabs(right);
511  });
512 
513  // Check balance
514  if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
515  peakItr--;
516  else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
517  peakItr++;
518 
519  // Special handling of the start tick for multiple hits
520  size_t hitCandidateStartTick = roiStartTick + std::distance(startItr, candStartItr);
521 
522  if (!hitCandidateVec.empty()) {
523  int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
524 
525  if (deltaTicks > 0) {
526  hitCandidateStartTick -= deltaTicks / 2;
527  hitCandidateVec.back().stopTick += deltaTicks / 2;
528  }
529  }
530 
531  hitCandidate.startTick = hitCandidateStartTick;
532  hitCandidate.stopTick = roiStartTick + std::distance(startItr, candStopItr);
533  hitCandidate.maxTick = roiStartTick + std::distance(startItr, maxItr);
534  hitCandidate.minTick = roiStartTick + std::distance(startItr, minItr);
535  hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
536  hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
537  hitCandidate.hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
538  hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
539  hitCandidate.hitHeight = hitCandidate.hitSigma *
540  (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
541 
542  hitCandidateVec.push_back(hitCandidate);
543  }
544  }
545 
546  // // The idea will be to find the largest deviation in the input derivative waveform as the starting point. Depending
547  // // on if a maximum or minimum, we search forward or backward to find the minimum or maximum that our extremum
548  // // corresponds to.
549  // std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair = std::minmax_element(startItr, stopItr);
550  //
551  // Waveform::const_iterator maxItr = minMaxPair.second;
552  // Waveform::const_iterator minItr = minMaxPair.first;
553  //
554  // // Use the larger of the two as the starting point and recover the nearest max or min
555  // if (std::fabs(*maxItr) > std::fabs(*minItr)) minItr = findNearestMin(maxItr, stopItr);
556  // else maxItr = findNearestMax(minItr,startItr);
557  //
558  // int deltaTicks = std::distance(maxItr,minItr);
559  // float range = *maxItr - *minItr;
560  //
561  // // At some point small rolling oscillations on the waveform need to be ignored...
562  // if (deltaTicks >= dTicksThreshold && range > dPeakThreshold)
563  // {
564  // // Need to back up to find zero crossing, this will be the starting point of our
565  // // candidate hit but also the endpoint of the pre sub-waveform we'll search next
566  // Waveform::const_iterator newEndItr = findStartTick(maxItr, startItr);
567  //
568  // int startTick = std::distance(startItr,newEndItr);
569  //
570  // // Now need to go forward to again get close to zero, this will then be the end point
571  // // of our candidate hit and the starting point for the post sub-waveform to search
572  // Waveform::const_iterator newStartItr = findStopTick(minItr, stopItr);
573  //
574  // int stopTick = std::distance(startItr,newStartItr);
575  //
576  // // Find hits in the section of the waveform leading up to this candidate hit
577  // if (startTick > 2)
578  // {
579  // // Special handling for merged hits
580  // if (*(newEndItr-1) > 0.) {dTicksThreshold = 2; dPeakThreshold = 0.; }
581  // else {dTicksThreshold = fMinDeltaTicks; dPeakThreshold = fMinDeltaPeaks;}
582  //
583  // findHitCandidates(startItr,newEndItr+1,roiStartTick,dTicksThreshold,dPeakThreshold,hitCandidateVec);
584  // }
585  //
586  // // Create a new hit candidate and store away
587  // HitCandidate_t hitCandidate;
588  //
589  // Waveform::const_iterator peakItr = std::min_element(maxItr,minItr,[](const auto& left, const auto& right){return std::fabs(left) < std::fabs(right);});
590  //
591  // // Check balance
592  // if (2 * std::distance(peakItr,minItr) < std::distance(maxItr,peakItr)) peakItr--;
593  // else if (2 * std::distance(maxItr,peakItr) < std::distance(peakItr,minItr)) peakItr++;
594  //
595  // // Special handling of the start tick for multiple hits
596  // size_t hitCandidateStartTick = roiStartTick + startTick;
597  //
598  // if (!hitCandidateVec.empty())
599  // {
600  // int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
601  //
602  // if (deltaTicks > 0)
603  // {
604  // hitCandidateStartTick -= deltaTicks / 2;
605  // hitCandidateVec.back().stopTick += deltaTicks / 2;
606  // }
607  // }
608  //
609  // hitCandidate.startTick = hitCandidateStartTick;
610  // hitCandidate.stopTick = roiStartTick + stopTick;
611  // hitCandidate.maxTick = roiStartTick + std::distance(startItr,maxItr);
612  // hitCandidate.minTick = roiStartTick + std::distance(startItr,minItr);
613  // hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
614  // hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
615  // hitCandidate.hitCenter = roiStartTick + std::distance(startItr,peakItr) + 0.5;
616  // hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
617  // hitCandidate.hitHeight = hitCandidate.hitSigma * (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
618  //
619  // hitCandidateVec.push_back(hitCandidate);
620  //
621  // // Finally, search the section of the waveform following this candidate for more hits
622  // if (std::distance(newStartItr,stopItr) > 2)
623  // {
624  // // Special handling for merged hits
625  // if (*(newStartItr+1) < 0.) {dTicksThreshold = 2; dPeakThreshold = 0.; }
626  // else {dTicksThreshold = fMinDeltaTicks; dPeakThreshold = fMinDeltaPeaks;}
627  //
628  // findHitCandidates(newStartItr,stopItr,roiStartTick + stopTick,dTicksThreshold,dPeakThreshold,hitCandidateVec);
629  // }
630  // }
631 
632  return;
633  }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
std::vector< CandHitParams > CandHitParamsVec
intermediate_table::const_iterator const_iterator
bool getListOfHitCandidates(Waveform::const_iterator, Waveform::const_iterator, int, float, CandHitParamsVec &) const
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitMorphological::findNearestMax ( Waveform::const_iterator  minItr,
Waveform::const_iterator  startItr 
) const
private

Definition at line 754 of file CandHitMorphological_tool.cc.

Referenced by getListOfHitCandidates().

757  {
758  // Set the internal loop variable...
759  Waveform::const_iterator lastItr = minItr;
760 
761  // One extra condition to watch for here, make sure we can actually back up!
762  if (std::distance(startItr, minItr) > 0) {
763  // Similar to searching for a maximum, we loop backward over ticks looking for the waveform to start decreasing
764  while ((lastItr - 1) != startItr) {
765  if (*(lastItr - 1) < *lastItr) break;
766 
767  lastItr--;
768  }
769  }
770 
771  return lastItr;
772  }
intermediate_table::const_iterator const_iterator
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitMorphological::findNearestMin ( Waveform::const_iterator  maxItr,
Waveform::const_iterator  stopItr 
) const
private

Definition at line 735 of file CandHitMorphological_tool.cc.

Referenced by getListOfHitCandidates().

738  {
739  // reset the min iterator and search forward to find the nearest minimum
740  Waveform::const_iterator lastItr = maxItr;
741 
742  // The strategy is simple...
743  // We are at a maximum so we search forward until we find the lowest negative point
744  while ((lastItr + 1) != stopItr) {
745  if (*(lastItr + 1) > *lastItr) break;
746 
747  lastItr++;
748  }
749 
750  // The minimum will be the last iterator value...
751  return lastItr;
752  }
intermediate_table::const_iterator const_iterator
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitMorphological::findStartTick ( Waveform::const_iterator  maxItr,
Waveform::const_iterator  startItr 
) const
private

Definition at line 774 of file CandHitMorphological_tool.cc.

Referenced by getListOfHitCandidates().

777  {
778  Waveform::const_iterator lastItr = maxItr;
779 
780  // If we can't back up then there is nothing to do
781  if (std::distance(startItr, lastItr) > 0) {
782  // In theory, we are starting at a maximum and want to find the "start" of the candidate peak
783  // Ideally we would look to search backward to the point where the (derivative) waveform crosses zero again.
784  // However, the complication is that we need to watch for the case where two peaks are merged together and
785  // we might run through another peak before crossing zero...
786  // So... loop until we hit the startItr...
787  Waveform::const_iterator loopItr = lastItr - 1;
788 
789  while (loopItr != startItr) {
790  // Ideal world case, we cross zero... but we might encounter a minimum... or an inflection point
791  if (*loopItr < 0. || !(*loopItr < *lastItr)) break;
792 
793  lastItr = loopItr--;
794  }
795  }
796  else
797  lastItr = startItr;
798 
799  return lastItr;
800  }
intermediate_table::const_iterator const_iterator
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitMorphological::findStopTick ( Waveform::const_iterator  minItr,
Waveform::const_iterator  stopItr 
) const
private

Definition at line 802 of file CandHitMorphological_tool.cc.

References DEFINE_ART_CLASS_TOOL.

Referenced by getListOfHitCandidates().

805  {
806  Waveform::const_iterator lastItr = minItr;
807 
808  // If we can't go forward then there is really nothing to do
809  if (std::distance(minItr, stopItr) > 1) {
810  // Pretty much the same strategy as for finding the start tick...
811  Waveform::const_iterator loopItr = lastItr + 1;
812 
813  while (loopItr != stopItr) {
814  // Ideal case that we have crossed zero coming from a minimum... but watch for a maximum as well
815  if (*loopItr > 0. || !(*loopItr > *lastItr)) break;
816 
817  lastItr = loopItr++;
818  }
819  }
820 
821  return lastItr;
822  }
intermediate_table::const_iterator const_iterator
bool reco_tool::CandHitMorphological::getListOfHitCandidates ( Waveform::const_iterator  startItr,
Waveform::const_iterator  stopItr,
int  dTicksThreshold,
float  dPeakThreshold,
CandHitParamsVec candHitParamsVec 
) const
private

Definition at line 635 of file CandHitMorphological_tool.cc.

References findNearestMax(), findNearestMin(), findStartTick(), findStopTick(), and fMinDeltaTicks.

Referenced by findHitCandidates().

640  {
641  // We'll check if any of our candidates meet the requirements so declare the result here
642  bool foundCandidate(false);
643 
644  int dTicks = std::distance(startItr, stopItr);
645 
646  // Search for candidate hits...
647  // But only if enough ticks
648  if (dTicks < fMinDeltaTicks) return foundCandidate;
649 
650  // Generally, the mission is simple... the goal is to find all possible combinations of maximum/minimum pairs in
651  // the input (presumed) derivative waveform. We can do this with a divice and conquer approach where we start by
652  // finding the largerst max or min and start from there
653  MaxMinPair minMaxPair = std::minmax_element(startItr, stopItr);
654 
655  Waveform::const_iterator maxItr = minMaxPair.second;
656  Waveform::const_iterator minItr = minMaxPair.first;
657 
658  // Use the larger of the two as the starting point and recover the nearest max or min
659  if (std::fabs(*maxItr) > std::fabs(*minItr))
660  minItr = findNearestMin(maxItr, stopItr);
661  else
662  maxItr = findNearestMax(minItr, startItr);
663 
664  int deltaTicks = std::distance(maxItr, minItr);
665  float range = *maxItr - *minItr;
666 
667  if (deltaTicks < 2) return foundCandidate;
668 
669  // Check if this particular max/min pair would meet the requirements...
670  if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) foundCandidate = true;
671 
672  // Need to back up to find zero crossing, this will be the starting point of our
673  // candidate hit but also the endpoint of the pre sub-waveform we'll search next
674  Waveform::const_iterator candStartItr = findStartTick(maxItr, startItr);
675 
676  // Now need to go forward to again get close to zero, this will then be the end point
677  // of our candidate hit and the starting point for the post sub-waveform to search
678  Waveform::const_iterator candStopItr = findStopTick(minItr, stopItr);
679 
680  // Call ourself to find hit candidates preceding this one
681  bool prevTicks = getListOfHitCandidates(
682  startItr, candStartItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
683 
684  // The above call will have populated the list of candidate max/min pairs preceding this one, so now add our contribution
685  candHitParamsVec.emplace_back(candStartItr, maxItr, minItr, candStopItr);
686 
687  // Now catch any that might follow this one
688  bool postTicks = getListOfHitCandidates(
689  candStopItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
690 
691  return foundCandidate || prevTicks || postTicks;
692  }
Waveform::const_iterator findNearestMin(Waveform::const_iterator, Waveform::const_iterator) const
Waveform::const_iterator findStopTick(Waveform::const_iterator, Waveform::const_iterator) const
intermediate_table::const_iterator const_iterator
Waveform::const_iterator findStartTick(Waveform::const_iterator, Waveform::const_iterator) const
bool getListOfHitCandidates(Waveform::const_iterator, Waveform::const_iterator, int, float, CandHitParamsVec &) const
std::pair< Waveform::const_iterator, Waveform::const_iterator > MaxMinPair
Waveform::const_iterator findNearestMax(Waveform::const_iterator, Waveform::const_iterator) const
void reco_tool::CandHitMorphological::MergeHitCandidates ( const recob::Wire::RegionsOfInterest_t::datarange_t &  ,
const HitCandidateVec hitCandidateVec,
MergeHitCandidateVec mergedHitsVec 
) const
overridevirtual

Implements reco_tool::ICandidateHitFinder.

Definition at line 694 of file CandHitMorphological_tool.cc.

References fMinHitHeight, and fNumInterveningTicks.

698  {
699  // If nothing on the input end then nothing to do
700  if (hitCandidateVec.empty()) return;
701 
702  // The idea is to group hits that "touch" so they can be part of common fit, those that
703  // don't "touch" are fit independently. So here we build the output vector to achieve that
704  // Get a container for the hits...
705  HitCandidateVec groupedHitVec;
706 
707  // Initialize the end of the last hit which we'll set to the first input hit's stop
708  size_t lastStopTick = hitCandidateVec.front().stopTick;
709 
710  // Step through the input hit candidates and group them by proximity
711  for (const auto& hitCandidate : hitCandidateVec) {
712  // Small pulse height hits should not be considered?
713  if (hitCandidate.hitHeight > fMinHitHeight) {
714  // Check condition that we have a new grouping
715  if (hitCandidate.startTick > lastStopTick + fNumInterveningTicks &&
716  !groupedHitVec.empty()) {
717  mergedHitsVec.emplace_back(groupedHitVec);
718 
719  groupedHitVec.clear();
720  }
721 
722  // Add the current hit to the current group
723  groupedHitVec.emplace_back(hitCandidate);
724 
725  lastStopTick = hitCandidate.stopTick;
726  }
727  }
728 
729  // Check end condition
730  if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
731 
732  return;
733  }
std::vector< HitCandidate > HitCandidateVec

Member Data Documentation

size_t reco_tool::CandHitMorphological::fChannelCnt
mutableprivate

Definition at line 116 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

const float reco_tool::CandHitMorphological::fDilationFraction
private

Definition at line 88 of file CandHitMorphological_tool.cc.

const float reco_tool::CandHitMorphological::fDilationThreshold
private

Definition at line 87 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates().

TH1F* reco_tool::CandHitMorphological::fDMaxDerivMinDerivHist
private

Definition at line 106 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

TH1F* reco_tool::CandHitMorphological::fDMaxTickMinTickHist
private

Definition at line 105 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

TH1F* reco_tool::CandHitMorphological::fDStopStartHist
private

Definition at line 104 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

const float reco_tool::CandHitMorphological::fErosionFraction
private

Definition at line 89 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates().

const float reco_tool::CandHitMorphological::fFitNSigmaFromCenter
private

Definition at line 99 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates().

const geo::GeometryCore* reco_tool::CandHitMorphological::fGeometry = lar::providerFrom<geo::Geometry>()
private

Definition at line 121 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates().

art::TFileDirectory* reco_tool::CandHitMorphological::fHistDirectory
private

Definition at line 101 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

size_t reco_tool::CandHitMorphological::fLastChannel
mutableprivate

Definition at line 115 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

TH1F* reco_tool::CandHitMorphological::fMaxDilationHist
private

Definition at line 108 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

TH1F* reco_tool::CandHitMorphological::fMaxDilEroRatHist
private

Definition at line 109 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

TH1F* reco_tool::CandHitMorphological::fMaxErosionHist
private

Definition at line 107 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

const float reco_tool::CandHitMorphological::fMinDeltaPeaks
private

Definition at line 91 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates().

const int reco_tool::CandHitMorphological::fMinDeltaTicks
private

Definition at line 90 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates(), and getListOfHitCandidates().

const float reco_tool::CandHitMorphological::fMinHitHeight
private

Definition at line 92 of file CandHitMorphological_tool.cc.

Referenced by MergeHitCandidates().

const size_t reco_tool::CandHitMorphological::fNumInterveningTicks
private

Definition at line 93 of file CandHitMorphological_tool.cc.

Referenced by MergeHitCandidates().

const bool reco_tool::CandHitMorphological::fOutputHistograms
private

Definition at line 95 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

const bool reco_tool::CandHitMorphological::fOutputWaveforms
private

Definition at line 97 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().

const size_t reco_tool::CandHitMorphological::fPlane
private

Definition at line 86 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological().

const int reco_tool::CandHitMorphological::fStructuringElement
private

Definition at line 94 of file CandHitMorphological_tool.cc.

Referenced by findHitCandidates().

std::unique_ptr<reco_tool::IWaveformTool> reco_tool::CandHitMorphological::fWaveformTool
private

Definition at line 119 of file CandHitMorphological_tool.cc.

Referenced by CandHitMorphological(), and findHitCandidates().


The documentation for this class was generated from the following file: