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

Public Types

using PeakTuple = std::tuple< size_t, size_t, size_t >
 
using PeakTupleVec = std::vector< PeakTuple >
 

Public Member Functions

 WaveformTools (const fhicl::ParameterSet &pset)
 
 ~WaveformTools ()
 
void configure (const fhicl::ParameterSet &pset) override
 
void triangleSmooth (const std::vector< float > &, std::vector< float > &, size_t=0) const override
 
void triangleSmooth (const std::vector< double > &, std::vector< double > &, size_t=0) const override
 
void medianSmooth (const std::vector< float > &, std::vector< float > &, size_t=3) const override
 
void medianSmooth (const std::vector< double > &, std::vector< double > &, size_t=3) const override
 
void getTruncatedMeanRMS (const std::vector< double > &, double &, double &, double &, int &) const override
 
void getTruncatedMeanRMS (const std::vector< float > &, float &, float &, float &, int &) const override
 
void firstDerivative (const std::vector< float > &, std::vector< float > &) const override
 
void firstDerivative (const std::vector< double > &, std::vector< double > &) const override
 
void findPeaks (std::vector< float >::iterator, std::vector< float >::iterator, PeakTupleVec &, float, size_t) const override
 
void findPeaks (std::vector< double >::iterator, std::vector< double >::iterator, PeakTupleVec &, double, size_t) const override
 
void getFFTPower (const std::vector< float > &inputVec, std::vector< float > &outputPowerVec) const override
 
void getFFTPower (const std::vector< double > &inputVec, std::vector< double > &outputPowerVec) const override
 
void getErosionDilationAverageDifference (const Waveform< short > &, int, HistogramMap &, Waveform< short > &, Waveform< short > &, Waveform< short > &, Waveform< short > &) const override
 
void getErosionDilationAverageDifference (const Waveform< float > &, int, HistogramMap &, Waveform< float > &, Waveform< float > &, Waveform< float > &, Waveform< float > &) const override
 
void getErosionDilationAverageDifference (const Waveform< double > &, int, HistogramMap &, Waveform< double > &, Waveform< double > &, Waveform< double > &, Waveform< double > &) const override
 
void getOpeningAndClosing (const Waveform< short > &, const Waveform< short > &, int, HistogramMap &, Waveform< short > &, Waveform< short > &) const override
 
void getOpeningAndClosing (const Waveform< float > &, const Waveform< float > &, int, HistogramMap &, Waveform< float > &, Waveform< float > &) const override
 
void getOpeningAndClosing (const Waveform< double > &, const Waveform< double > &, int, HistogramMap &, Waveform< double > &, Waveform< double > &) const override
 

Private Member Functions

template<typename T >
void triangleSmooth (const std::vector< T > &, std::vector< T > &, size_t=0) const
 
template<typename T >
void medianSmooth (const std::vector< T > &, std::vector< T > &, size_t=3) const
 
template<typename T >
void getTruncatedMeanRMS (const std::vector< T > &, T &, T &, T &, int &) const
 
template<typename T >
void firstDerivative (const std::vector< T > &, std::vector< T > &) const
 
template<typename T >
void findPeaks (typename std::vector< T >::iterator, typename std::vector< T >::iterator, PeakTupleVec &, T, size_t) const
 
template<typename T >
void getErosionDilationAverageDifference (const Waveform< T > &, int, HistogramMap &, Waveform< T > &, Waveform< T > &, Waveform< T > &, Waveform< T > &) const
 
template<typename T >
void getOpeningAndClosing (const Waveform< T > &, const Waveform< T > &, int, HistogramMap &, Waveform< T > &, Waveform< T > &) const
 

Detailed Description

Definition at line 16 of file WaveformTools_tool.cc.

Member Typedef Documentation

using reco_tool::WaveformTools::PeakTuple = std::tuple<size_t, size_t, size_t>

Definition at line 24 of file WaveformTools_tool.cc.

Definition at line 25 of file WaveformTools_tool.cc.

Constructor & Destructor Documentation

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

Definition at line 137 of file WaveformTools_tool.cc.

References configure().

138  {
139  configure(pset);
140  }
void configure(const fhicl::ParameterSet &pset) override
reco_tool::WaveformTools::~WaveformTools ( )
inline

Definition at line 20 of file WaveformTools_tool.cc.

References configure().

20 {}

Member Function Documentation

void reco_tool::WaveformTools::configure ( const fhicl::ParameterSet pset)
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 142 of file WaveformTools_tool.cc.

Referenced by WaveformTools(), and ~WaveformTools().

143  {
144  // Start by recovering the parameters
145  // fThisPlane = pset.get<size_t>("Plane");
146 
147  return;
148  }
void reco_tool::WaveformTools::findPeaks ( std::vector< float >::iterator  startItr,
std::vector< float >::iterator  stopItr,
PeakTupleVec peakTupleVec,
float  threshold,
size_t  firstTick 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 383 of file WaveformTools_tool.cc.

Referenced by findPeaks().

388  {
389  findPeaks<float>(startItr, stopItr, peakTupleVec, threshold, firstTick);
390 
391  return;
392  }
void reco_tool::WaveformTools::findPeaks ( std::vector< double >::iterator  startItr,
std::vector< double >::iterator  stopItr,
PeakTupleVec peakTupleVec,
double  threshold,
size_t  firstTick 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 372 of file WaveformTools_tool.cc.

377  {
378  findPeaks<double>(startItr, stopItr, peakTupleVec, threshold, firstTick);
379 
380  return;
381  }
template<typename T >
void reco_tool::WaveformTools::findPeaks ( typename std::vector< T >::iterator  startItr,
typename std::vector< T >::iterator  stopItr,
PeakTupleVec peakTupleVec,
threshold,
size_t  firstTick 
) const
private

Definition at line 395 of file WaveformTools_tool.cc.

References findPeaks(), art::left(), and art::right().

400  {
401  // Need a minimum distance or else nothing to do
402  if (std::distance(startItr, stopItr) > 4) {
403  // This is a divide and conquer algorithm, start by finding the maximum element.
404  typename std::vector<T>::iterator firstItr =
405  std::max_element(startItr, stopItr, [](float left, float right) {
406  return std::fabs(left) < std::fabs(right);
407  });
408 
409  // Are we over threshold?
410  if (std::fabs(*firstItr) > threshold) {
411  // What am I thinking?
412  // First task is to find the "other" lobe max point
413  // Set one to the "first", the other to the "second"
414  // Search backward from first to find start point, forward from second to find end point
415  // Set mid point between first and second as "peak"?
416  typename std::vector<T>::iterator secondItr = firstItr;
417 
418  // Assume if max bin is positive then second lobe is later
419  if (*firstItr > 0) {
420  typename std::vector<T>::iterator tempItr = secondItr;
421 
422  while (tempItr != stopItr) {
423  if (*++tempItr < -threshold) {
424  if (*tempItr < *secondItr) secondItr = tempItr;
425  }
426  else if (secondItr != firstItr)
427  break;
428  }
429  }
430  // Otherwise it goes the other way
431  else {
432  typename std::vector<T>::iterator tempItr = secondItr;
433 
434  while (tempItr != startItr) {
435  if (*--tempItr > threshold) {
436  if (*tempItr > *secondItr) secondItr = tempItr;
437  }
438  else if (secondItr != firstItr)
439  break;
440  }
441 
442  std::swap(firstItr, secondItr);
443  }
444 
445  // It might that no real pulse was found
446  if (firstItr != secondItr) {
447  // Get the "peak" position
448  size_t peakBin =
449  std::distance(startItr, firstItr) + std::distance(firstItr, secondItr) / 2;
450 
451  // Advance (forward or backward) the first and second iterators to get back to zero crossing
452  while (firstItr != startItr)
453  if (*--firstItr < 0.) break;
454  while (secondItr != stopItr)
455  if (*++secondItr > 0.) break;
456 
457  size_t firstBin = std::distance(startItr, firstItr);
458  size_t lastBin = std::distance(startItr, secondItr);
459 
460  // Find leading peaks
461  findPeaks(startItr, firstItr, peakTupleVec, threshold, firstTick);
462 
463  // Save this peak
464  peakTupleVec.push_back(
465  PeakTuple(firstBin + firstTick, peakBin + firstTick, lastBin + firstTick));
466 
467  // Find downstream peaks
468  findPeaks(secondItr,
469  stopItr,
470  peakTupleVec,
471  threshold,
472  firstTick + std::distance(startItr, secondItr));
473  }
474  }
475  }
476 
477  return;
478  }
intermediate_table::iterator iterator
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
void findPeaks(std::vector< float >::iterator, std::vector< float >::iterator, PeakTupleVec &, float, size_t) const override
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
std::tuple< size_t, size_t, size_t > PeakTuple
void reco_tool::WaveformTools::firstDerivative ( const std::vector< float > &  inputVec,
std::vector< float > &  derivVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 352 of file WaveformTools_tool.cc.

354  {
355  firstDerivative<float>(inputVec, derivVec);
356 
357  return;
358  }
void reco_tool::WaveformTools::firstDerivative ( const std::vector< double > &  inputVec,
std::vector< double > &  derivVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 344 of file WaveformTools_tool.cc.

346  {
347  firstDerivative<double>(inputVec, derivVec);
348 
349  return;
350  }
template<typename T >
void reco_tool::WaveformTools::firstDerivative ( const std::vector< T > &  inputVec,
std::vector< T > &  derivVec 
) const
private

Definition at line 361 of file WaveformTools_tool.cc.

363  {
364  derivVec.resize(inputVec.size(), 0.);
365 
366  for (size_t idx = 1; idx < derivVec.size() - 1; idx++)
367  derivVec.at(idx) = 0.5 * (inputVec.at(idx + 1) - inputVec.at(idx - 1));
368 
369  return;
370  }
void reco_tool::WaveformTools::getErosionDilationAverageDifference ( const Waveform< short > &  waveform,
int  structuringElement,
HistogramMap histogramMap,
Waveform< short > &  erosionVec,
Waveform< short > &  dilationVec,
Waveform< short > &  averageVec,
Waveform< short > &  differenceVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 530 of file WaveformTools_tool.cc.

537  {
538  getErosionDilationAverageDifference<short>(waveform,
539  structuringElement,
540  histogramMap,
541  erosionVec,
542  dilationVec,
543  averageVec,
544  differenceVec);
545 
546  return;
547  }
void reco_tool::WaveformTools::getErosionDilationAverageDifference ( const Waveform< float > &  waveform,
int  structuringElement,
HistogramMap histogramMap,
Waveform< float > &  erosionVec,
Waveform< float > &  dilationVec,
Waveform< float > &  averageVec,
Waveform< float > &  differenceVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 549 of file WaveformTools_tool.cc.

556  {
557  getErosionDilationAverageDifference<float>(waveform,
558  structuringElement,
559  histogramMap,
560  erosionVec,
561  dilationVec,
562  averageVec,
563  differenceVec);
564 
565  return;
566  }
void reco_tool::WaveformTools::getErosionDilationAverageDifference ( const Waveform< double > &  waveform,
int  structuringElement,
HistogramMap histogramMap,
Waveform< double > &  erosionVec,
Waveform< double > &  dilationVec,
Waveform< double > &  averageVec,
Waveform< double > &  differenceVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 568 of file WaveformTools_tool.cc.

575  {
576  getErosionDilationAverageDifference<double>(waveform,
577  structuringElement,
578  histogramMap,
579  erosionVec,
580  dilationVec,
581  averageVec,
582  differenceVec);
583 
584  return;
585  }
template<typename T >
void reco_tool::WaveformTools::getErosionDilationAverageDifference ( const Waveform< T > &  inputWaveform,
int  structuringElement,
HistogramMap histogramMap,
Waveform< T > &  erosionVec,
Waveform< T > &  dilationVec,
Waveform< T > &  averageVec,
Waveform< T > &  differenceVec 
) const
private

Definition at line 588 of file WaveformTools_tool.cc.

References reco_tool::AVERAGE, reco_tool::DIFFERENCE, reco_tool::DILATION, reco_tool::EROSION, and reco_tool::WAVEFORM.

595  {
596  // Set the window size
597  int halfWindowSize(structuringElement / 2);
598 
599  // Initialize min and max elements
601  minMaxItr =
602  std::minmax_element(inputWaveform.begin(), inputWaveform.begin() + halfWindowSize);
603 
604  typename Waveform<T>::const_iterator minElementItr = minMaxItr.first;
605  typename Waveform<T>::const_iterator maxElementItr = minMaxItr.second;
606 
607  // Initialize the erosion and dilation vectors
608  erosionVec.resize(inputWaveform.size());
609  dilationVec.resize(inputWaveform.size());
610  averageVec.resize(inputWaveform.size());
611  differenceVec.resize(inputWaveform.size());
612 
613  // Now loop through remaining elements and complete the vectors
614  typename Waveform<T>::iterator minItr = erosionVec.begin();
615  typename Waveform<T>::iterator maxItr = dilationVec.begin();
616  typename Waveform<T>::iterator aveItr = averageVec.begin();
617  typename Waveform<T>::iterator difItr = differenceVec.begin();
618 
619  for (typename Waveform<T>::const_iterator inputItr = inputWaveform.begin();
620  inputItr != inputWaveform.end();
621  inputItr++) {
622  // There are two conditions to check:
623  // 1) is the current min/max element outside the current window?
624  // 2) is the new element smaller/larger than the current min/max?
625 
626  // Make sure we are not running off the end of the vector
627  if (std::distance(inputItr, inputWaveform.end()) > halfWindowSize) {
628  if (std::distance(minElementItr, inputItr) >= halfWindowSize)
629  minElementItr =
630  std::min_element(inputItr - halfWindowSize + 1, inputItr + halfWindowSize + 1);
631  else if (*(inputItr + halfWindowSize) < *minElementItr)
632  minElementItr = inputItr + halfWindowSize;
633 
634  if (std::distance(maxElementItr, inputItr) >= halfWindowSize)
635  maxElementItr =
636  std::max_element(inputItr - halfWindowSize + 1, inputItr + halfWindowSize + 1);
637  else if (*(inputItr + halfWindowSize) > *maxElementItr)
638  maxElementItr = inputItr + halfWindowSize;
639  }
640 
641  // Update the vectors
642  *minItr++ = *minElementItr;
643  *maxItr++ = *maxElementItr;
644  *aveItr++ = 0.5 * (*maxElementItr + *minElementItr);
645  *difItr++ = *maxElementItr - *minElementItr;
646 
647  if (!histogramMap.empty()) {
648  int curBin = std::distance(inputWaveform.begin(), inputItr);
649 
650  histogramMap.at(WAVEFORM)->Fill(curBin, *inputItr);
651  histogramMap.at(EROSION)->Fill(curBin, *minElementItr);
652  histogramMap.at(DILATION)->Fill(curBin, *maxElementItr);
653  histogramMap.at(AVERAGE)->Fill(curBin, 0.5 * (*maxElementItr + *minElementItr));
654  histogramMap.at(DIFFERENCE)->Fill(curBin, *maxElementItr - *minElementItr);
655  }
656  }
657 
658  return;
659  }
intermediate_table::iterator iterator
intermediate_table::const_iterator const_iterator
void reco_tool::WaveformTools::getFFTPower ( const std::vector< float > &  inputVec,
std::vector< float > &  outputPowerVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 480 of file WaveformTools_tool.cc.

482  {
483  std::vector<double> inputDoubleVec(inputVec.size());
484  std::vector<double> outputDoubleVec(inputVec.size() / 2);
485 
486  std::copy(inputVec.begin(), inputVec.end(), inputDoubleVec.begin());
487 
488  getFFTPower(inputDoubleVec, outputDoubleVec);
489 
490  if (outputDoubleVec.size() != outputPowerVec.size())
491  outputPowerVec.resize(outputDoubleVec.size());
492 
493  std::copy(outputDoubleVec.begin(), outputDoubleVec.end(), outputPowerVec.begin());
494 
495  return;
496  }
void getFFTPower(const std::vector< float > &inputVec, std::vector< float > &outputPowerVec) const override
void reco_tool::WaveformTools::getFFTPower ( const std::vector< double > &  inputVec,
std::vector< double > &  outputPowerVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 498 of file WaveformTools_tool.cc.

500  {
501  // Get the FFT of the response
502  int fftDataSize = inputVec.size();
503 
504  TVirtualFFT* fftr2c = TVirtualFFT::FFT(1, &fftDataSize, "R2C");
505 
506  fftr2c->SetPoints(inputVec.data());
507  fftr2c->Transform();
508 
509  // Recover the results so we can compute the power spectrum
510  size_t halfFFTDataSize(fftDataSize / 2 + 1);
511 
512  std::vector<double> realVals(halfFFTDataSize);
513  std::vector<double> imaginaryVals(halfFFTDataSize);
514 
515  fftr2c->GetPointsComplex(realVals.data(), imaginaryVals.data());
516 
517  if (outputPowerVec.size() != halfFFTDataSize) outputPowerVec.resize(halfFFTDataSize, 0.);
518 
519  std::transform(realVals.begin(),
520  realVals.begin() + halfFFTDataSize,
521  imaginaryVals.begin(),
522  outputPowerVec.begin(),
523  [](const double& real, const double& imaginary) {
524  return std::sqrt(real * real + imaginary * imaginary);
525  });
526 
527  return;
528  }
void reco_tool::WaveformTools::getOpeningAndClosing ( const Waveform< short > &  erosionVec,
const Waveform< short > &  dilationVec,
int  structuringElement,
HistogramMap histogramMap,
Waveform< short > &  openingVec,
Waveform< short > &  closingVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 661 of file WaveformTools_tool.cc.

667  {
668  getOpeningAndClosing<short>(
669  erosionVec, dilationVec, structuringElement, histogramMap, openingVec, closingVec);
670 
671  return;
672  }
void reco_tool::WaveformTools::getOpeningAndClosing ( const Waveform< float > &  erosionVec,
const Waveform< float > &  dilationVec,
int  structuringElement,
HistogramMap histogramMap,
Waveform< float > &  openingVec,
Waveform< float > &  closingVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 674 of file WaveformTools_tool.cc.

680  {
681  getOpeningAndClosing<float>(
682  erosionVec, dilationVec, structuringElement, histogramMap, openingVec, closingVec);
683 
684  return;
685  }
void reco_tool::WaveformTools::getOpeningAndClosing ( const Waveform< double > &  erosionVec,
const Waveform< double > &  dilationVec,
int  structuringElement,
HistogramMap histogramMap,
Waveform< double > &  openingVec,
Waveform< double > &  closingVec 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 687 of file WaveformTools_tool.cc.

693  {
694  getOpeningAndClosing<double>(
695  erosionVec, dilationVec, structuringElement, histogramMap, openingVec, closingVec);
696 
697  return;
698  }
template<typename T >
void reco_tool::WaveformTools::getOpeningAndClosing ( const Waveform< T > &  erosionVec,
const Waveform< T > &  dilationVec,
int  structuringElement,
HistogramMap histogramMap,
Waveform< T > &  openingVec,
Waveform< T > &  closingVec 
) const
private

Definition at line 701 of file WaveformTools_tool.cc.

References reco_tool::CLOSING, DEFINE_ART_CLASS_TOOL, reco_tool::DOPENCLOSING, and reco_tool::OPENING.

707  {
708  // Set the window size
709  int halfWindowSize(structuringElement / 2);
710 
711  // Start with the opening, here we get the max element in the input erosion vector
712  typename Waveform<T>::const_iterator maxElementItr =
713  std::max_element(erosionVec.begin(), erosionVec.begin() + halfWindowSize);
714 
715  // Initialize the opening vector
716  openingVec.resize(erosionVec.size());
717 
718  // Now loop through remaining elements and complete the vectors
719  typename Waveform<T>::iterator maxItr = openingVec.begin();
720 
721  for (typename Waveform<T>::const_iterator inputItr = erosionVec.begin();
722  inputItr != erosionVec.end();
723  inputItr++) {
724  // There are two conditions to check:
725  // 1) is the current min/max element outside the current window?
726  // 2) is the new element smaller/larger than the current min/max?
727 
728  // Make sure we are not running off the end of the vector
729  if (std::distance(inputItr, erosionVec.end()) > halfWindowSize) {
730  if (std::distance(maxElementItr, inputItr) >= halfWindowSize)
731  maxElementItr =
732  std::max_element(inputItr - halfWindowSize + 1, inputItr + halfWindowSize + 1);
733  else if (*(inputItr + halfWindowSize) > *maxElementItr)
734  maxElementItr = inputItr + halfWindowSize;
735  }
736 
737  // Update the vectors
738  *maxItr++ = *maxElementItr;
739 
740  if (!histogramMap.empty()) {
741  int curBin = std::distance(erosionVec.begin(), inputItr);
742 
743  histogramMap.at(OPENING)->Fill(curBin, *maxElementItr);
744  }
745  }
746 
747  // Now go with the closling, here we get the min element in the input dilation vector
748  typename Waveform<T>::const_iterator minElementItr =
749  std::min_element(dilationVec.begin(), dilationVec.begin() + halfWindowSize);
750 
751  // Initialize the opening and closing vectors
752  closingVec.resize(dilationVec.size());
753 
754  // Now loop through remaining elements and complete the vectors
755  typename Waveform<T>::iterator minItr = closingVec.begin();
756 
757  for (typename Waveform<T>::const_iterator inputItr = dilationVec.begin();
758  inputItr != dilationVec.end();
759  inputItr++) {
760  // There are two conditions to check:
761  // 1) is the current min/max element outside the current window?
762  // 2) is the new element smaller/larger than the current min/max?
763 
764  // Make sure we are not running off the end of the vector
765  if (std::distance(inputItr, dilationVec.end()) > halfWindowSize) {
766  if (std::distance(minElementItr, inputItr) >= halfWindowSize)
767  minElementItr =
768  std::min_element(inputItr - halfWindowSize + 1, inputItr + halfWindowSize + 1);
769  else if (*(inputItr + halfWindowSize) < *minElementItr)
770  minElementItr = inputItr + halfWindowSize;
771  }
772 
773  // Update the vectors
774  *minItr++ = *minElementItr;
775 
776  if (!histogramMap.empty()) {
777  int curBin = std::distance(dilationVec.begin(), inputItr);
778 
779  histogramMap.at(CLOSING)->Fill(curBin, *minElementItr);
780  histogramMap.at(DOPENCLOSING)->Fill(curBin, *minElementItr - openingVec.at(curBin));
781  }
782  }
783 
784  return;
785  }
intermediate_table::iterator iterator
intermediate_table::const_iterator const_iterator
void reco_tool::WaveformTools::getTruncatedMeanRMS ( const std::vector< double > &  waveform,
double &  mean,
double &  rmsFull,
double &  rmsTrunc,
int &  nTrunc 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 258 of file WaveformTools_tool.cc.

References pmtana::mean().

263  {
264  getTruncatedMeanRMS<double>(waveform, mean, rmsFull, rmsTrunc, nTrunc);
265  }
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
void reco_tool::WaveformTools::getTruncatedMeanRMS ( const std::vector< float > &  waveform,
float &  mean,
float &  rmsFull,
float &  rmsTrunc,
int &  nTrunc 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 267 of file WaveformTools_tool.cc.

References pmtana::mean().

272  {
273  getTruncatedMeanRMS<float>(waveform, mean, rmsFull, rmsTrunc, nTrunc);
274  }
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
template<typename T >
void reco_tool::WaveformTools::getTruncatedMeanRMS ( const std::vector< T > &  waveform,
T &  mean,
T &  rmsFull,
T &  rmsTrunc,
int &  nTrunc 
) const
private

Definition at line 277 of file WaveformTools_tool.cc.

References art::left(), and art::right().

282  {
283  // We need to get a reliable estimate of the mean and can't assume the input waveform will be ~zero mean...
284  // Basic idea is to find the most probable value in the ROI presented to us
285  // From that we can develop an average of the true baseline of the ROI.
286  // To do that we employ a map based scheme
287  std::map<int, int> frequencyMap;
288  int mpCount(0);
289  int mpVal(0);
290 
291  for (const auto& val : waveform) {
292  int intVal = std::round(4. * val);
293 
294  frequencyMap[intVal]++;
295 
296  if (frequencyMap.at(intVal) > mpCount) {
297  mpCount = frequencyMap.at(intVal);
298  mpVal = intVal;
299  }
300  }
301 
302  // take a weighted average of two neighbor bins
303  int meanCnt = 0;
304  int meanSum = 0;
305  int binRange = std::min(16, int(frequencyMap.size() / 2 + 1));
306 
307  for (int idx = -binRange; idx <= binRange; idx++) {
308  std::map<int, int>::iterator neighborItr = frequencyMap.find(mpVal + idx);
309 
310  if (neighborItr != frequencyMap.end() && 5 * neighborItr->second > mpCount) {
311  meanSum += neighborItr->first * neighborItr->second;
312  meanCnt += neighborItr->second;
313  }
314  }
315 
316  mean = 0.25 * T(meanSum) / T(meanCnt); // Note that bins were expanded by a factor of 4 above
317 
318  // do rms calculation - the old fashioned way and over all adc values
319  typename std::vector<T> locWaveform = waveform;
320 
321  std::transform(locWaveform.begin(),
322  locWaveform.end(),
323  locWaveform.begin(),
324  std::bind(std::minus<T>(), std::placeholders::_1, mean));
325 
326  // sort in ascending order so we can truncate the sume
327  std::sort(locWaveform.begin(), locWaveform.end(), [](const auto& left, const auto& right) {
328  return std::fabs(left) < std::fabs(right);
329  });
330 
331  // recalculate the rms for truncation
332  rmsFull = std::inner_product(locWaveform.begin(), locWaveform.end(), locWaveform.begin(), 0.);
333  rmsFull = std::sqrt(std::max(T(0.), rmsFull / T(locWaveform.size())));
334 
335  // recalculate the rms for truncation
336  rmsTrunc = std::inner_product(
337  locWaveform.begin(), locWaveform.begin() + meanCnt, locWaveform.begin(), 0.);
338  rmsTrunc = std::sqrt(std::max(T(0.), rmsTrunc / T(meanCnt)));
339  nTrunc = meanCnt;
340 
341  return;
342  }
intermediate_table::iterator iterator
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
void reco_tool::WaveformTools::medianSmooth ( const std::vector< float > &  inputVec,
std::vector< float > &  smoothVec,
size_t  nBins = 3 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 199 of file WaveformTools_tool.cc.

202  {
203  medianSmooth<float>(inputVec, smoothVec, nBins);
204 
205  return;
206  }
void reco_tool::WaveformTools::medianSmooth ( const std::vector< double > &  inputVec,
std::vector< double > &  smoothVec,
size_t  nBins = 3 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 208 of file WaveformTools_tool.cc.

211  {
212  medianSmooth<double>(inputVec, smoothVec, nBins);
213 
214  return;
215  }
template<typename T >
void reco_tool::WaveformTools::medianSmooth ( const std::vector< T > &  inputVec,
std::vector< T > &  smoothVec,
size_t  nBins = 3 
) const
private

Definition at line 218 of file WaveformTools_tool.cc.

221  {
222  // For our purposes, nBins must be odd
223  if (nBins % 2 == 0) nBins++;
224 
225  // Make sure the input vector is right sized
226  if (inputVec.size() != smoothVec.size()) smoothVec.resize(inputVec.size());
227 
228  // Basic set up
229  typename std::vector<T> medianVec(nBins);
230  typename std::vector<T>::const_iterator startItr = inputVec.begin();
231  typename std::vector<T>::const_iterator stopItr = startItr;
232 
233  std::advance(stopItr, inputVec.size() - nBins);
234 
235  size_t medianBin = nBins / 2;
236  size_t smoothBin = medianBin;
237 
238  // First bins are not smoothed
239  std::copy(startItr, startItr + medianBin, smoothVec.begin());
240 
241  while (std::distance(startItr, stopItr) > 0) {
242  std::copy(startItr, startItr + nBins, medianVec.begin());
243  std::sort(medianVec.begin(), medianVec.end());
244 
245  T medianVal = medianVec[medianBin];
246 
247  smoothVec[smoothBin++] = medianVal;
248 
249  startItr++;
250  }
251 
252  // Last bins are not smoothed
253  std::copy(startItr + medianBin, inputVec.end(), smoothVec.begin() + smoothBin);
254 
255  return;
256  }
intermediate_table::const_iterator const_iterator
void reco_tool::WaveformTools::triangleSmooth ( const std::vector< float > &  inputVec,
std::vector< float > &  smoothVec,
size_t  lowestBin = 0 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 159 of file WaveformTools_tool.cc.

162  {
163  triangleSmooth<float>(inputVec, smoothVec, lowestBin);
164 
165  return;
166  }
void reco_tool::WaveformTools::triangleSmooth ( const std::vector< double > &  inputVec,
std::vector< double > &  smoothVec,
size_t  lowestBin = 0 
) const
overridevirtual

Implements reco_tool::IWaveformTool.

Definition at line 150 of file WaveformTools_tool.cc.

153  {
154  triangleSmooth<double>(inputVec, smoothVec, lowestBin);
155 
156  return;
157  }
template<typename T >
void reco_tool::WaveformTools::triangleSmooth ( const std::vector< T > &  inputVec,
std::vector< T > &  smoothVec,
size_t  lowestBin = 0 
) const
private

Definition at line 169 of file WaveformTools_tool.cc.

172  {
173  if (inputVec.size() != smoothVec.size()) smoothVec.resize(inputVec.size());
174 
175  // Watch for edge condition
176  if (inputVec.size() > 4) {
177  std::copy(inputVec.begin(), inputVec.begin() + 2 + lowestBin, smoothVec.begin());
178  std::copy(inputVec.end() - 2, inputVec.end(), smoothVec.end() - 2);
179 
180  typename std::vector<T>::iterator curItr = smoothVec.begin() + 2 + lowestBin;
181  typename std::vector<T>::const_iterator curInItr = inputVec.begin() + 1 + lowestBin;
182  typename std::vector<T>::const_iterator stopInItr = inputVec.end() - 3;
183 
184  while (curInItr++ != stopInItr) {
185  // Take the weighted average of five consecutive points centered on current point
186  T newVal = (*(curInItr - 2) + 2. * *(curInItr - 1) + 3. * *curInItr + 2. * *(curInItr + 1) +
187  *(curInItr + 2)) /
188  9.;
189 
190  *curItr++ = newVal;
191  }
192  }
193  else
194  std::copy(inputVec.begin(), inputVec.end(), smoothVec.begin());
195 
196  return;
197  }
intermediate_table::iterator iterator
intermediate_table::const_iterator const_iterator

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