LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
CandHitMorphological_tool.cc
Go to the documentation of this file.
1 // MT note: This implementation is not thread-safe.
6 
10 
12 #include "art/Utilities/Globals.h"
15 #include "art_root_io/TFileService.h"
16 #include "cetlib_except/exception.h"
18 
19 #include "TProfile.h"
20 
21 #include <cmath>
22 
23 namespace reco_tool {
24 
26  public:
27  explicit CandHitMorphological(const fhicl::ParameterSet& pset);
28 
29  void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&,
30  const size_t,
31  const size_t,
32  const size_t,
33  HitCandidateVec&) const override;
34 
35  void MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&,
36  const HitCandidateVec&,
37  MergeHitCandidateVec&) const override;
38 
39  private:
40  // Internal functions
41  //< Top level hit finding using erosion/dilation vectors
43  Waveform::const_iterator, //< derivative
45  Waveform::const_iterator, //< erosion
47  Waveform::const_iterator, //< dilation
48  const size_t,
49  float,
50  HitCandidateVec&) const;
51 
52  //< Fine grain hit finding within candidate peak regions using derivative method
55  const size_t,
56  int,
57  float,
58  HitCandidateVec&) const;
59 
60  //< For a given range, return the list of max/min pairs
61  using MaxMinPair = std::pair<Waveform::const_iterator, Waveform::const_iterator>;
62  using CandHitParams = std::tuple<Waveform::const_iterator,
65  Waveform::const_iterator>;
66  using CandHitParamsVec = std::vector<CandHitParams>;
67 
68  bool getListOfHitCandidates(Waveform::const_iterator,
69  Waveform::const_iterator,
70  int,
71  float,
72  CandHitParamsVec&) const;
73 
74  // Finding the nearest maximum/minimum from current point
75  Waveform::const_iterator findNearestMax(Waveform::const_iterator,
76  Waveform::const_iterator) const;
77  Waveform::const_iterator findNearestMin(Waveform::const_iterator,
78  Waveform::const_iterator) const;
79 
80  // handle finding the "start" and "stop" of a candidate hit
81  Waveform::const_iterator findStartTick(Waveform::const_iterator,
82  Waveform::const_iterator) const;
83  Waveform::const_iterator findStopTick(Waveform::const_iterator, Waveform::const_iterator) const;
84 
85  // some fhicl control variables
86  const size_t fPlane; //< Identifies the plane this tool is meant to operate on
87  const float fDilationThreshold; //< Dilation threshold
88  const float fDilationFraction; //< Fraction of max dilation to set for min dilation
89  const float fErosionFraction; //< Fraction of max dilation value to set min erosion
90  const int fMinDeltaTicks; //< minimum ticks from max to min to consider
91  const float fMinDeltaPeaks; //< minimum maximum to minimum peak distance
92  const float fMinHitHeight; //< Drop candidate hits with height less than this
93  const size_t fNumInterveningTicks; //< Number ticks between candidate hits to merge
94  const int fStructuringElement; //< Window size for morphologcial filter
95  const bool fOutputHistograms; //< If true will generate summary style histograms
96  const bool
97  fOutputWaveforms; //< If true will output waveform related info <<< very big output file!
98  const float
99  fFitNSigmaFromCenter; //< Limit ticks to fit to NSigma from hit center; not applied if zero or negative
100 
101  art::TFileDirectory* fHistDirectory;
102 
103  // Global histograms
104  TH1F* fDStopStartHist; //< Basically keeps track of the length of hit regions
105  TH1F* fDMaxTickMinTickHist; //< This will be a measure of the width of candidate hits
106  TH1F* fDMaxDerivMinDerivHist; //< This is the difference peak to peak of derivative for cand hit
107  TH1F* fMaxErosionHist; //< Keep track of the maximum erosion
108  TH1F* fMaxDilationHist; //< Keep track of the maximum dilation
109  TH1F* fMaxDilEroRatHist; //< Ratio of the maxima of the two
110 
111  //MT note: The mutable data members are only used in the histogram filling functions
112  //and histogram filling can only be done in single-threaded mode.
113  //Will need to consider design changes if this behavior changes.
114  mutable size_t
115  fLastChannel; //< Kludge to keep track of last channel when histogramming in effect
116  mutable size_t fChannelCnt; //< Counts the number of times a channel is used (assumed in order)
117 
118  //< All of the real work is done in the waveform tool
119  std::unique_ptr<reco_tool::IWaveformTool> fWaveformTool;
120 
123  };
124 
125  //----------------------------------------------------------------------
126  // Constructor.
128  : fPlane(pset.get<size_t>("Plane", 0))
129  , fDilationThreshold(pset.get<float>("DilationThreshold", 4.))
130  , fDilationFraction(pset.get<float>("DilationFraction", 0.75))
131  , fErosionFraction(pset.get<float>("ErosionFraction", 0.2))
132  , fMinDeltaTicks(pset.get<int>("MinDeltaTicks", 0))
133  , fMinDeltaPeaks(pset.get<float>("MinDeltaPeaks", 0.025))
134  , fMinHitHeight(pset.get<float>("MinHitHeight", 1.0))
135  , fNumInterveningTicks(pset.get<size_t>("NumInterveningTicks", 6))
136  , fStructuringElement(pset.get<int>("StructuringElement", 20))
137  , fOutputHistograms(pset.get<bool>("OutputHistograms", false))
138  , fOutputWaveforms(pset.get<bool>("OutputWaveforms", false))
139  , fFitNSigmaFromCenter(pset.get<float>("FitNSigmaFromCenter", 5.))
140  {
141 
142  if (art::Globals::instance()->nthreads() > 1u) {
143  if (fOutputHistograms) {
145  << "Cannot fill histograms when multiple threads configured, please set "
146  "fOutputHistograms to false or change number of threads to 1\n";
147  }
148 
149  if (fOutputWaveforms) {
151  << "Cannot write output waveforms when multiple threads configured, please set "
152  "fOutputHistograms to false or change number of threads to 1\n";
153  }
154  }
155  // Recover the baseline tool
156  fWaveformTool =
157  art::make_tool<reco_tool::IWaveformTool>(pset.get<fhicl::ParameterSet>("WaveformAlgs"));
158 
159  // Set the last channel to some nonsensical value
160  fLastChannel = std::numeric_limits<size_t>::max();
161  fChannelCnt = 0;
162 
163  // If asked, define the global histograms
164  if (fOutputHistograms) {
165  // Access ART's TFileService, which will handle creating and writing
166  // histograms and n-tuples for us.
168 
169  fHistDirectory = tfs.get();
170 
171  // Make a directory for these histograms
172  art::TFileDirectory dir = fHistDirectory->mkdir(Form("HitPlane_%1zu", fPlane));
173 
175  dir.make<TH1F>(Form("DStopStart_%1zu", fPlane), ";Delta Stop/Start;", 100, 0., 100.);
177  dir.make<TH1F>(Form("DMaxTMinT_%1zu", fPlane), ";Delta Max/Min Tick;", 100, 0., 100.);
179  dir.make<TH1F>(Form("DMaxDMinD_%1zu", fPlane), ";Delta Max/Min Deriv;", 200, 0., 100.);
181  dir.make<TH1F>(Form("MaxErosion_%1zu", fPlane), ";Max Erosion;", 200, -50., 150.);
183  dir.make<TH1F>(Form("MaxDilation_%1zu", fPlane), ";Max Dilation;", 200, -50., 150.);
185  dir.make<TH1F>(Form("MaxDilEroRat_%1zu", fPlane), ";Max Dil/Ero;", 200, -1., 1.);
186  }
187 
188  return;
189  }
190 
192  const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange,
193  const size_t roiStartTick,
194  const size_t channel,
195  const size_t eventCount,
196  HitCandidateVec& hitCandidateVec) const
197  {
198  // In this case we want to find hit candidates based on the derivative of of the input waveform
199  // We get this from our waveform algs too...
200  Waveform rawDerivativeVec;
201  Waveform derivativeVec;
202 
203  // Recover the actual waveform
204  const Waveform& waveform = dataRange.data();
205 
206  fWaveformTool->firstDerivative(waveform, rawDerivativeVec);
207  fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
208 
209  // Now we get the erosion/dilation vectors
210  Waveform erosionVec;
211  Waveform dilationVec;
212  Waveform averageVec;
213  Waveform differenceVec;
214 
215  reco_tool::HistogramMap histogramMap;
216 
217  // Compute the morphological filter vectors
218  fWaveformTool->getErosionDilationAverageDifference(waveform,
220  histogramMap,
221  erosionVec,
222  dilationVec,
223  averageVec,
224  differenceVec);
225 
226  // Now find the hits
227  findHitCandidates(derivativeVec.begin(),
228  derivativeVec.end(),
229  erosionVec.begin(),
230  erosionVec.end(),
231  dilationVec.begin(),
232  dilationVec.end(),
233  roiStartTick,
235  hitCandidateVec);
236 
237  // Limit start and stop tick to the neighborhood of the peak
238  if (fFitNSigmaFromCenter > 0) {
239  for (auto& hc : hitCandidateVec) {
240  auto startCand = hc.hitCenter - fFitNSigmaFromCenter * hc.hitSigma;
241  if (startCand >= 0) hc.startTick = std::max(hc.startTick, size_t(startCand));
242  hc.stopTick =
243  std::min(hc.stopTick, size_t(hc.hitCenter + fFitNSigmaFromCenter * hc.hitSigma));
244  }
245  }
246 
247  // Reset the hit height from the input waveform
248  for (auto& hitCandidate : hitCandidateVec) {
249  size_t centerIdx = hitCandidate.hitCenter;
250 
251  hitCandidate.hitHeight = waveform.at(centerIdx);
252  }
253 
254  // Keep track of histograms if requested
255  if (fOutputWaveforms) {
256  // Recover the details...
257  std::vector<geo::WireID> wids = fWireReadoutGeom->ChannelToWire(channel);
258  size_t plane = wids[0].Plane;
259  size_t cryo = wids[0].Cryostat;
260  size_t tpc = wids[0].TPC;
261  size_t wire = wids[0].Wire;
262 
263  if (channel != fLastChannel) fChannelCnt = 0;
264 
265  // Make a directory for these histograms
266  art::TFileDirectory dir = fHistDirectory->mkdir(
267  Form("Event%04zu/c%1zuT%1zuP%1zu/Wire_%05zu", eventCount, cryo, tpc, plane, wire));
268 
269  size_t waveformSize = waveform.size();
270  size_t waveStart = dataRange.begin_index();
271 
272  TProfile* waveHist =
273  dir.make<TProfile>(Form("HWfm_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
274  "Waveform",
275  waveformSize,
276  0,
277  waveformSize,
278  -500.,
279  500.);
280  TProfile* derivHist =
281  dir.make<TProfile>(Form("HDer_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
282  "Derivative",
283  waveformSize,
284  0,
285  waveformSize,
286  -500.,
287  500.);
288  TProfile* erosionHist =
289  dir.make<TProfile>(Form("HEro_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
290  "Erosion",
291  waveformSize,
292  0,
293  waveformSize,
294  -500.,
295  500.);
296  TProfile* dilationHist =
297  dir.make<TProfile>(Form("HDil_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
298  "Dilation",
299  waveformSize,
300  0,
301  waveformSize,
302  -500.,
303  500.);
304  TProfile* candHitHist =
305  dir.make<TProfile>(Form("HCan_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
306  "Cand Hits",
307  waveformSize,
308  0,
309  waveformSize,
310  -500.,
311  500.);
312  TProfile* maxDerivHist =
313  dir.make<TProfile>(Form("HMax_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
314  "Maxima",
315  waveformSize,
316  0,
317  waveformSize,
318  -500.,
319  500.);
320  TProfile* strtStopHist =
321  dir.make<TProfile>(Form("HSSS_%03zu_roiStart-%05zu", fChannelCnt, waveStart),
322  "Start/Stop",
323  waveformSize,
324  0,
325  waveformSize,
326  -500.,
327  500.);
328 
329  // Fill wave/derivative
330  for (size_t idx = 0; idx < waveform.size(); idx++) {
331  waveHist->Fill(roiStartTick + idx, waveform.at(idx));
332  derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
333  erosionHist->Fill(roiStartTick + idx, erosionVec.at(idx));
334  dilationHist->Fill(roiStartTick + idx, dilationVec.at(idx));
335  }
336 
337  // Fill hits
338  for (const auto& hitCandidate : hitCandidateVec) {
339  candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
340  maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
341  maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
342  strtStopHist->Fill(hitCandidate.startTick, waveform.at(hitCandidate.startTick));
343  strtStopHist->Fill(hitCandidate.stopTick, waveform.at(hitCandidate.stopTick));
344  }
345 
346  fLastChannel = channel;
347  fChannelCnt++;
348  }
349 
350  if (fOutputHistograms) {
351  // Fill hits
352  for (const auto& hitCandidate : hitCandidateVec) {
353  fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
354  fDMaxTickMinTickHist->Fill(hitCandidate.minTick - hitCandidate.maxTick, 1.);
355  fDMaxDerivMinDerivHist->Fill(hitCandidate.maxDerivative - hitCandidate.minDerivative, 1.);
356  }
357 
358  // Get the max dilation/erosion
359  Waveform::const_iterator maxDilationItr =
360  std::max_element(dilationVec.begin(), dilationVec.end());
361  Waveform::const_iterator maxErosionItr =
362  std::max_element(erosionVec.begin(), erosionVec.end());
363 
364  float dilEroRat(1.);
365 
366  if (std::abs(*maxDilationItr) > 0.) dilEroRat = *maxErosionItr / *maxDilationItr;
367 
368  fMaxErosionHist->Fill(*maxErosionItr, 1.);
369  fMaxDilationHist->Fill(*maxDilationItr, 1.);
370  fMaxDilEroRatHist->Fill(dilEroRat, 1.);
371  }
372 
373  return;
374  }
375 
377  Waveform::const_iterator derivStopItr,
378  Waveform::const_iterator erosionStartItr,
379  Waveform::const_iterator erosionStopItr,
380  Waveform::const_iterator dilationStartItr,
381  Waveform::const_iterator dilationStopItr,
382  const size_t roiStartTick,
383  float dilationThreshold,
384  HitCandidateVec& hitCandidateVec) const
385  {
386  // This function aims to use the erosion/dilation vectors to find candidate hit regions
387  // Once armed with a region then the "standard" differential approach is used to return the candidate peaks
388 
389  // Don't do anything if not enough ticks
390  int ticksInInputWaveform = std::distance(derivStartItr, derivStopItr);
391 
392  if (ticksInInputWaveform < fMinDeltaTicks) return;
393 
394  // This function is recursive, we start by finding the largest element of the dilation vector
395  Waveform::const_iterator maxItr = std::max_element(dilationStartItr, dilationStopItr);
396  float maxVal = *maxItr;
397 
398  // Check that the peak is of reasonable height...
399  if (maxVal < dilationThreshold) return;
400 
401  int maxBin = std::distance(dilationStartItr, maxItr);
402 
403  // The candidate hit region we want lies between the two nearest minima to the maximum we just found
404  // subject to the condition that the erosion vector has return to less than zero
405  Waveform::const_iterator firstDerItr = derivStartItr + maxBin;
406  Waveform::const_iterator erosionItr = erosionStartItr + maxBin;
407 
408  float firstDerivValue = -1.;
409  float erosionCutValue = fErosionFraction * maxVal;
410 
411  // Search for starting point
412  while (firstDerItr != derivStartItr) {
413  // Look for the turnover point
414  if (*erosionItr-- < erosionCutValue) {
415  // We are looking for the zero crossing signifying a minimum value in the waveform
416  // (so the previous derivative < 0 while current is > 0)
417  // We are moving "backwards" so the current value <= 0, the previous value > 0
418  if (*firstDerItr * firstDerivValue <= 0. && firstDerivValue > 0.) break;
419  }
420 
421  firstDerivValue = *firstDerItr--;
422  }
423 
424  // Set the start bin
425  int hitRegionStart = std::distance(derivStartItr, firstDerItr);
426 
427  // Now go the other way
428  Waveform::const_iterator lastDerItr = derivStartItr + maxBin;
429 
430  // Reset the local variables
431  float lastDerivValue = 1.;
432 
433  erosionItr = erosionStartItr + maxBin;
434 
435  // Search for starting point
436  while (lastDerItr != derivStopItr) {
437  if (*erosionItr++ <= erosionCutValue) {
438  // We are again looking for the zero crossing signifying a minimum value in the waveform
439  // This time we are moving forward, so test is that previous value < 0, new value >= 0
440  if (*lastDerItr * lastDerivValue <= 0. && lastDerivValue < 0.) break;
441  }
442 
443  lastDerivValue = *lastDerItr++;
444  }
445 
446  // Set the stop bin
447  int hitRegionStop = std::distance(derivStartItr, lastDerItr);
448 
449  // Recursive call to find any hits in front of where we are now
450  if (hitRegionStart > fMinDeltaTicks)
451  findHitCandidates(derivStartItr,
452  derivStartItr + hitRegionStart,
453  erosionStartItr,
454  erosionStartItr + hitRegionStart,
455  dilationStartItr,
456  dilationStartItr + hitRegionStart,
457  roiStartTick,
459  hitCandidateVec);
460 
461  // Call the differential hit finding to get the actual hits within the region
462  findHitCandidates(derivStartItr + hitRegionStart,
463  derivStartItr + hitRegionStop,
464  roiStartTick + hitRegionStart,
467  hitCandidateVec);
468 
469  // Now call ourselves again to find any hits trailing the region we just identified
470  if (std::distance(lastDerItr, derivStopItr) > fMinDeltaTicks)
471  findHitCandidates(derivStartItr + hitRegionStop,
472  derivStopItr,
473  erosionStartItr + hitRegionStop,
474  erosionStopItr,
475  dilationStartItr + hitRegionStop,
476  dilationStopItr,
477  roiStartTick + hitRegionStop,
479  hitCandidateVec);
480 
481  return;
482  }
483 
485  Waveform::const_iterator stopItr,
486  const size_t roiStartTick,
487  int dTicksThreshold,
488  float dPeakThreshold,
489  HitCandidateVec& hitCandidateVec) const
490  {
491  // Search for candidate hits...
492  // Strategy is to get the list of all possible max/min pairs of the input derivative vector and then
493  // look for candidate hits in that list
494  CandHitParamsVec candHitParamsVec;
495 
497  startItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec)) {
498  // We've been given a list of candidate hits so now convert to hits
499  // Version one... simply convert all the candidates
500  for (const auto& tuple : candHitParamsVec) {
501  // Create a new hit candidate and store away
502  HitCandidate hitCandidate;
503 
504  Waveform::const_iterator candStartItr = std::get<0>(tuple);
505  Waveform::const_iterator maxItr = std::get<1>(tuple);
506  Waveform::const_iterator minItr = std::get<2>(tuple);
507  Waveform::const_iterator candStopItr = std::get<3>(tuple);
508 
509  Waveform::const_iterator peakItr =
510  std::min_element(maxItr, minItr, [](const auto& left, const auto& right) {
511  return std::fabs(left) < std::fabs(right);
512  });
513 
514  // Check balance
515  if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
516  peakItr--;
517  else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
518  peakItr++;
519 
520  // Special handling of the start tick for multiple hits
521  size_t hitCandidateStartTick = roiStartTick + std::distance(startItr, candStartItr);
522 
523  if (!hitCandidateVec.empty()) {
524  int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
525 
526  if (deltaTicks > 0) {
527  hitCandidateStartTick -= deltaTicks / 2;
528  hitCandidateVec.back().stopTick += deltaTicks / 2;
529  }
530  }
531 
532  hitCandidate.startTick = hitCandidateStartTick;
533  hitCandidate.stopTick = roiStartTick + std::distance(startItr, candStopItr);
534  hitCandidate.maxTick = roiStartTick + std::distance(startItr, maxItr);
535  hitCandidate.minTick = roiStartTick + std::distance(startItr, minItr);
536  hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
537  hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
538  hitCandidate.hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
539  hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
540  hitCandidate.hitHeight = hitCandidate.hitSigma *
541  (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
542 
543  hitCandidateVec.push_back(hitCandidate);
544  }
545  }
546 
547  // // The idea will be to find the largest deviation in the input derivative waveform as the starting point. Depending
548  // // on if a maximum or minimum, we search forward or backward to find the minimum or maximum that our extremum
549  // // corresponds to.
550  // std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair = std::minmax_element(startItr, stopItr);
551  //
552  // Waveform::const_iterator maxItr = minMaxPair.second;
553  // Waveform::const_iterator minItr = minMaxPair.first;
554  //
555  // // Use the larger of the two as the starting point and recover the nearest max or min
556  // if (std::fabs(*maxItr) > std::fabs(*minItr)) minItr = findNearestMin(maxItr, stopItr);
557  // else maxItr = findNearestMax(minItr,startItr);
558  //
559  // int deltaTicks = std::distance(maxItr,minItr);
560  // float range = *maxItr - *minItr;
561  //
562  // // At some point small rolling oscillations on the waveform need to be ignored...
563  // if (deltaTicks >= dTicksThreshold && range > dPeakThreshold)
564  // {
565  // // Need to back up to find zero crossing, this will be the starting point of our
566  // // candidate hit but also the endpoint of the pre sub-waveform we'll search next
567  // Waveform::const_iterator newEndItr = findStartTick(maxItr, startItr);
568  //
569  // int startTick = std::distance(startItr,newEndItr);
570  //
571  // // Now need to go forward to again get close to zero, this will then be the end point
572  // // of our candidate hit and the starting point for the post sub-waveform to search
573  // Waveform::const_iterator newStartItr = findStopTick(minItr, stopItr);
574  //
575  // int stopTick = std::distance(startItr,newStartItr);
576  //
577  // // Find hits in the section of the waveform leading up to this candidate hit
578  // if (startTick > 2)
579  // {
580  // // Special handling for merged hits
581  // if (*(newEndItr-1) > 0.) {dTicksThreshold = 2; dPeakThreshold = 0.; }
582  // else {dTicksThreshold = fMinDeltaTicks; dPeakThreshold = fMinDeltaPeaks;}
583  //
584  // findHitCandidates(startItr,newEndItr+1,roiStartTick,dTicksThreshold,dPeakThreshold,hitCandidateVec);
585  // }
586  //
587  // // Create a new hit candidate and store away
588  // HitCandidate_t hitCandidate;
589  //
590  // Waveform::const_iterator peakItr = std::min_element(maxItr,minItr,[](const auto& left, const auto& right){return std::fabs(left) < std::fabs(right);});
591  //
592  // // Check balance
593  // if (2 * std::distance(peakItr,minItr) < std::distance(maxItr,peakItr)) peakItr--;
594  // else if (2 * std::distance(maxItr,peakItr) < std::distance(peakItr,minItr)) peakItr++;
595  //
596  // // Special handling of the start tick for multiple hits
597  // size_t hitCandidateStartTick = roiStartTick + startTick;
598  //
599  // if (!hitCandidateVec.empty())
600  // {
601  // int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
602  //
603  // if (deltaTicks > 0)
604  // {
605  // hitCandidateStartTick -= deltaTicks / 2;
606  // hitCandidateVec.back().stopTick += deltaTicks / 2;
607  // }
608  // }
609  //
610  // hitCandidate.startTick = hitCandidateStartTick;
611  // hitCandidate.stopTick = roiStartTick + stopTick;
612  // hitCandidate.maxTick = roiStartTick + std::distance(startItr,maxItr);
613  // hitCandidate.minTick = roiStartTick + std::distance(startItr,minItr);
614  // hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
615  // hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
616  // hitCandidate.hitCenter = roiStartTick + std::distance(startItr,peakItr) + 0.5;
617  // hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
618  // hitCandidate.hitHeight = hitCandidate.hitSigma * (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
619  //
620  // hitCandidateVec.push_back(hitCandidate);
621  //
622  // // Finally, search the section of the waveform following this candidate for more hits
623  // if (std::distance(newStartItr,stopItr) > 2)
624  // {
625  // // Special handling for merged hits
626  // if (*(newStartItr+1) < 0.) {dTicksThreshold = 2; dPeakThreshold = 0.; }
627  // else {dTicksThreshold = fMinDeltaTicks; dPeakThreshold = fMinDeltaPeaks;}
628  //
629  // findHitCandidates(newStartItr,stopItr,roiStartTick + stopTick,dTicksThreshold,dPeakThreshold,hitCandidateVec);
630  // }
631  // }
632 
633  return;
634  }
635 
637  Waveform::const_iterator stopItr,
638  int dTicksThreshold,
639  float dPeakThreshold,
640  CandHitParamsVec& candHitParamsVec) const
641  {
642  // We'll check if any of our candidates meet the requirements so declare the result here
643  bool foundCandidate(false);
644 
645  int dTicks = std::distance(startItr, stopItr);
646 
647  // Search for candidate hits...
648  // But only if enough ticks
649  if (dTicks < fMinDeltaTicks) return foundCandidate;
650 
651  // Generally, the mission is simple... the goal is to find all possible combinations of maximum/minimum pairs in
652  // the input (presumed) derivative waveform. We can do this with a divice and conquer approach where we start by
653  // finding the largerst max or min and start from there
654  MaxMinPair minMaxPair = std::minmax_element(startItr, stopItr);
655 
656  Waveform::const_iterator maxItr = minMaxPair.second;
657  Waveform::const_iterator minItr = minMaxPair.first;
658 
659  // Use the larger of the two as the starting point and recover the nearest max or min
660  if (std::fabs(*maxItr) > std::fabs(*minItr))
661  minItr = findNearestMin(maxItr, stopItr);
662  else
663  maxItr = findNearestMax(minItr, startItr);
664 
665  int deltaTicks = std::distance(maxItr, minItr);
666  float range = *maxItr - *minItr;
667 
668  if (deltaTicks < 2) return foundCandidate;
669 
670  // Check if this particular max/min pair would meet the requirements...
671  if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) foundCandidate = true;
672 
673  // Need to back up to find zero crossing, this will be the starting point of our
674  // candidate hit but also the endpoint of the pre sub-waveform we'll search next
675  Waveform::const_iterator candStartItr = findStartTick(maxItr, startItr);
676 
677  // Now need to go forward to again get close to zero, this will then be the end point
678  // of our candidate hit and the starting point for the post sub-waveform to search
679  Waveform::const_iterator candStopItr = findStopTick(minItr, stopItr);
680 
681  // Call ourself to find hit candidates preceding this one
682  bool prevTicks = getListOfHitCandidates(
683  startItr, candStartItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
684 
685  // The above call will have populated the list of candidate max/min pairs preceding this one, so now add our contribution
686  candHitParamsVec.emplace_back(candStartItr, maxItr, minItr, candStopItr);
687 
688  // Now catch any that might follow this one
689  bool postTicks = getListOfHitCandidates(
690  candStopItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
691 
692  return foundCandidate || prevTicks || postTicks;
693  }
694 
696  const recob::Wire::RegionsOfInterest_t::datarange_t&,
697  const HitCandidateVec& hitCandidateVec,
698  MergeHitCandidateVec& mergedHitsVec) const
699  {
700  // If nothing on the input end then nothing to do
701  if (hitCandidateVec.empty()) return;
702 
703  // The idea is to group hits that "touch" so they can be part of common fit, those that
704  // don't "touch" are fit independently. So here we build the output vector to achieve that
705  // Get a container for the hits...
706  HitCandidateVec groupedHitVec;
707 
708  // Initialize the end of the last hit which we'll set to the first input hit's stop
709  size_t lastStopTick = hitCandidateVec.front().stopTick;
710 
711  // Step through the input hit candidates and group them by proximity
712  for (const auto& hitCandidate : hitCandidateVec) {
713  // Small pulse height hits should not be considered?
714  if (hitCandidate.hitHeight > fMinHitHeight) {
715  // Check condition that we have a new grouping
716  if (hitCandidate.startTick > lastStopTick + fNumInterveningTicks &&
717  !groupedHitVec.empty()) {
718  mergedHitsVec.emplace_back(groupedHitVec);
719 
720  groupedHitVec.clear();
721  }
722 
723  // Add the current hit to the current group
724  groupedHitVec.emplace_back(hitCandidate);
725 
726  lastStopTick = hitCandidate.stopTick;
727  }
728  }
729 
730  // Check end condition
731  if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
732 
733  return;
734  }
735 
738  Waveform::const_iterator stopItr) const
739  {
740  // reset the min iterator and search forward to find the nearest minimum
741  Waveform::const_iterator lastItr = maxItr;
742 
743  // The strategy is simple...
744  // We are at a maximum so we search forward until we find the lowest negative point
745  while ((lastItr + 1) != stopItr) {
746  if (*(lastItr + 1) > *lastItr) break;
747 
748  lastItr++;
749  }
750 
751  // The minimum will be the last iterator value...
752  return lastItr;
753  }
754 
757  Waveform::const_iterator startItr) const
758  {
759  // Set the internal loop variable...
760  Waveform::const_iterator lastItr = minItr;
761 
762  // One extra condition to watch for here, make sure we can actually back up!
763  if (std::distance(startItr, minItr) > 0) {
764  // Similar to searching for a maximum, we loop backward over ticks looking for the waveform to start decreasing
765  while ((lastItr - 1) != startItr) {
766  if (*(lastItr - 1) < *lastItr) break;
767 
768  lastItr--;
769  }
770  }
771 
772  return lastItr;
773  }
774 
777  Waveform::const_iterator startItr) const
778  {
779  Waveform::const_iterator lastItr = maxItr;
780 
781  // If we can't back up then there is nothing to do
782  if (std::distance(startItr, lastItr) > 0) {
783  // In theory, we are starting at a maximum and want to find the "start" of the candidate peak
784  // Ideally we would look to search backward to the point where the (derivative) waveform crosses zero again.
785  // However, the complication is that we need to watch for the case where two peaks are merged together and
786  // we might run through another peak before crossing zero...
787  // So... loop until we hit the startItr...
788  Waveform::const_iterator loopItr = lastItr - 1;
789 
790  while (loopItr != startItr) {
791  // Ideal world case, we cross zero... but we might encounter a minimum... or an inflection point
792  if (*loopItr < 0. || !(*loopItr < *lastItr)) break;
793 
794  lastItr = loopItr--;
795  }
796  }
797  else
798  lastItr = startItr;
799 
800  return lastItr;
801  }
802 
805  Waveform::const_iterator stopItr) const
806  {
807  Waveform::const_iterator lastItr = minItr;
808 
809  // If we can't go forward then there is really nothing to do
810  if (std::distance(minItr, stopItr) > 1) {
811  // Pretty much the same strategy as for finding the start tick...
812  Waveform::const_iterator loopItr = lastItr + 1;
813 
814  while (loopItr != stopItr) {
815  // Ideal case that we have crossed zero coming from a minimum... but watch for a maximum as well
816  if (*loopItr > 0. || !(*loopItr > *lastItr)) break;
817 
818  lastItr = loopItr++;
819  }
820  }
821 
822  return lastItr;
823  }
824 
826 }
Waveform::const_iterator findNearestMin(Waveform::const_iterator, Waveform::const_iterator) const
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
std::vector< CandHitParams > CandHitParamsVec
T * get() const
Definition: ServiceHandle.h:69
This is the interface class for tools/algorithms that perform various operations on waveforms...
constexpr auto abs(T v)
Returns the absolute value of the argument.
Waveform::const_iterator findStopTick(Waveform::const_iterator, Waveform::const_iterator) const
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
intermediate_table::const_iterator const_iterator
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
Waveform::const_iterator findStartTick(Waveform::const_iterator, Waveform::const_iterator) const
std::tuple< Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator, Waveform::const_iterator > CandHitParams
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
Interface for a class providing readout channel mapping to geometry.
bool getListOfHitCandidates(Waveform::const_iterator, Waveform::const_iterator, int, float, CandHitParamsVec &) const
T get(std::string const &key) const
Definition: ParameterSet.h:314
const geo::WireReadoutGeom * fWireReadoutGeom
std::map< int, TProfile * > HistogramMap
Definition: IWaveformTool.h:42
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
TDirectory * dir
Definition: macro.C:5
std::pair< Waveform::const_iterator, Waveform::const_iterator > MaxMinPair
This provides an interface for tools which are tasked with finding candidate hits on input waveforms...
Waveform::const_iterator findNearestMax(Waveform::const_iterator, Waveform::const_iterator) const
void MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const HitCandidateVec &, MergeHitCandidateVec &) const override
static Globals * instance()
Definition: Globals.cc:17
std::vector< HitCandidateVec > MergeHitCandidateVec
art framework interface to geometry description
std::vector< HitCandidate > HitCandidateVec