LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Protected Member Functions | Private Attributes | List of all members
pmtana::PedAlgoRollingMean Class Referenceabstract

#include "PedAlgoRollingMean.h"

Inheritance diagram for pmtana::PedAlgoRollingMean:
pmtana::PMTPedestalBase

Public Member Functions

 PedAlgoRollingMean (const std::string name="PedRollingMean")
 Default constructor. More...
 
 PedAlgoRollingMean (const fhicl::ParameterSet &pset, const std::string name="PedRollingMean")
 Alternative ctor. More...
 
const std::string & Name () const
 Name getter. More...
 
bool Evaluate (const pmtana::Waveform_t &wf)
 Method to compute a pedestal. More...
 
double Mean (size_t i) const
 Getter of the pedestal mean value. More...
 
const pmtana::PedestalMean_t & Mean () const
 Getter of the pedestal mean value. More...
 
double Sigma (size_t i) const
 Getter of the pedestal standard deviation. More...
 
const pmtana::PedestalSigma_t & Sigma () const
 Getter of the pedestal standard deviation. More...
 

Protected Member Functions

bool ComputePedestal (const pmtana::Waveform_t &wf, pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)
 Method to compute a pedestal of the input waveform using "nsample" ADC samples from "start" index. More...
 
virtual bool ComputePedestal (const ::pmtana::Waveform_t &wf, pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)=0
 

Private Attributes

size_t _sample_size
 
float _max_sigma
 
float _ped_range_max
 
float _ped_range_min
 
double _threshold
 
double _diff_threshold
 
double _diff_adc_count
 
int _n_presamples
 

Detailed Description

A class that calculates pedestal mean & standard deviation (here and elsewhere called as "RMS").

Definition at line 33 of file PedAlgoRollingMean.h.

Constructor & Destructor Documentation

pmtana::PedAlgoRollingMean::PedAlgoRollingMean ( const std::string  name = "PedRollingMean")

Default constructor.

Definition at line 16 of file PedAlgoRollingMean.cxx.

16  : PMTPedestalBase(name)
17  //*****************************************************************
18  {
19  srand(static_cast<unsigned int>(time(0)));
20  }
pmtana::PMTPedestalBase::PMTPedestalBase
PMTPedestalBase(std::string name="noname")
Default constructor.
Definition: PMTPedestalBase.cxx:13
pmtana::PedAlgoRollingMean::PedAlgoRollingMean ( const fhicl::ParameterSet &  pset,
const std::string  name = "PedRollingMean" 
)

Alternative ctor.

Definition at line 23 of file PedAlgoRollingMean.cxx.

References _diff_adc_count, _diff_threshold, _max_sigma, _n_presamples, _ped_range_max, _ped_range_min, _sample_size, _threshold, and fhicl::ParameterSet::get().

27  : PMTPedestalBase(name)
28  //############################################################
29  {
30 
31  _sample_size = pset.get<size_t>("SampleSize");
32  _max_sigma = pset.get<float>("MaxSigma");
33  _ped_range_max = pset.get<float>("PedRangeMax");
34  _ped_range_min = pset.get<float>("PedRangeMin");
35 
36  // _range = pset.get<int> ("RandomRange");
37  // _divisions = pset.get<double>("RandomRangeDivisions");
38  _threshold = pset.get<double>("Threshold");
39  _diff_threshold = pset.get<double>("DiffBetweenGapsThreshold");
40  _diff_adc_count = pset.get<double>("DiffADCCounts");
41 
42  _n_presamples = pset.get<int>("NPrePostSamples");
43 
44  //_random_shift = pset.get<double>("RandomRangeShift");
45  // Random seed number generator
46  //srand(static_cast<unsigned int>(time(0)));
47  }
pmtana::PMTPedestalBase::PMTPedestalBase
PMTPedestalBase(std::string name="noname")
Default constructor.
Definition: PMTPedestalBase.cxx:13
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:314

Member Function Documentation

bool pmtana::PedAlgoRollingMean::ComputePedestal ( const pmtana::Waveform_t &  wf,
pmtana::PedestalMean_t &  mean_v,
pmtana::PedestalSigma_t &  sigma_v 
)
protected

Method to compute a pedestal of the input waveform using "nsample" ADC samples from "start" index.

Definition at line 50 of file PedAlgoRollingMean.cxx.

References _diff_adc_count, _diff_threshold, _max_sigma, _n_presamples, _ped_range_max, _sample_size, _threshold, pmtana::BinnedMaxOccurrence(), pmtana::edge_aware_mean(), pmtana::mean(), and pmtana::std().

54  {
55 
56  // parameters
57  if (wf.size() <= (_sample_size * 2)) return false;
58 
59  const size_t window_size = _sample_size * 2;
60 
61  // middle mean
62  for (size_t i = 0; i < wf.size(); ++i) {
63 
64  mean_v[i] = 0;
65  sigma_v[i] = 0;
66 
67  if (i < _sample_size || i >= (wf.size() - _sample_size)) continue;
68 
69  mean_v[i] = mean(wf, i - _sample_size, window_size);
70  sigma_v[i] = std(wf, mean_v[i], i - _sample_size, window_size);
71  }
72 
73  // front mean
74  for (size_t i = 0; i < _sample_size; ++i) {
75 
76  mean_v[i] = mean_v[_sample_size];
77  sigma_v[i] = sigma_v[_sample_size];
78  }
79  // tail mean
80  for (size_t i = (wf.size() - _sample_size); i < wf.size(); ++i) {
81 
82  mean_v[i] = mean_v[wf.size() - _sample_size - 1];
83  sigma_v[i] = sigma_v[wf.size() - _sample_size - 1];
84  }
85 
86  float best_sigma = 1.1e9;
87  size_t best_sigma_index = 0;
88  size_t num_good_adc = 0;
89 
90  for (size_t i = 0; i < sigma_v.size(); ++i) {
91  // Only consider adcs which mean is in the allowed range
92  auto const& mean = mean_v[i];
93 
94  if (mean < _ped_range_min || mean > _ped_range_max) continue;
95 
96  auto const& sigma = sigma_v[i];
97  if (sigma < best_sigma) {
98  best_sigma = sigma;
99  best_sigma_index = i;
100  }
101 
102  if (sigma < _max_sigma) num_good_adc += 1;
103  }
104 
105  if (num_good_adc < 1) {
106  std::cerr << "\033[93m<<" << __FUNCTION__
107  << ">>\033[00m Could not find good pedestal at all..." << std::endl;
108  return false;
109  }
110 
111  // If not enough # of good mean indices, use the best guess within this waveform
112  if (best_sigma > _max_sigma || num_good_adc < 3) {
113  for (size_t i = 0; i < mean_v.size(); ++i) {
114  mean_v[i] = mean_v.at(best_sigma_index);
115  sigma_v[i] = sigma_v.at(best_sigma_index);
116  }
117 
118  return true;
119  }
120 
121  // Else do extrapolation, or random seed depending on what we find...
122 
123  unsigned nbins = 1000;
124 
126  const auto mode_mean = BinnedMaxOccurrence(mean_v, nbins);
127  const auto mode_sigma = BinnedMaxOccurrence(sigma_v, nbins);
128 
129  //auto mode_mean = BinnedMaxTH1D(mean_v ,nbins);
130  //auto mode_sigma = BinnedMaxTH1D(sigma_v,nbins);
131 
132  //std::cout<<mode_mean<<" +/- "<<mode_sigma<<std::endl;
133 
134  double const diff_cutoff = std::max(_diff_threshold * mode_sigma, _diff_adc_count);
135 
136  int last_good_index = -1;
137 
138  for (size_t i = 0; i < wf.size(); ++i) {
139 
140  auto const mean = mean_v[i];
141  auto const sigma = sigma_v[i];
142 
143  // if(sigma <= _max_sigma && mean < _ped_range_max && mean > _ped_range_min) {
144  // not sure if this works well for basline that is "linear" seen by David K
145 
146  if (sigma <= _threshold * mode_sigma && fabs(mean - mode_mean) <= _threshold * mode_sigma) {
147 
148  if (last_good_index < 0) {
149  last_good_index = (int)i;
150  continue;
151  }
152 
153  if ((last_good_index + 1) < (int)i) {
154 
155  auto diff = fabs(mean_v.at(last_good_index) - mean);
156 
157  if (diff > diff_cutoff) {
158  //this should become generic interpolation function, for now lets leave.
159  float slope = (mean - mean_v.at(last_good_index)) / (float(i - last_good_index));
160 
161  for (size_t j = last_good_index + 1; j < i; ++j) {
162  mean_v.at(j) = slope * (float(j - last_good_index)) + mean_v.at(last_good_index);
163  sigma_v.at(j) = mode_sigma;
164  }
165  }
166  else {
167  //difference roughly the same lets fill with constant baseline
168 
169  auto presample_mean =
170  edge_aware_mean(wf, last_good_index - _n_presamples, last_good_index);
171  auto postsample_mean = edge_aware_mean(wf, i, i + _n_presamples);
172 
173  auto diff_pre = fabs(presample_mean - mode_mean);
174  auto diff_post = fabs(postsample_mean - mode_mean);
175 
176  auto constant_mean = diff_pre <= diff_post ? presample_mean : postsample_mean;
177 
178  for (size_t j = last_good_index + 1; j < i; ++j) {
179  //mean_v.at(j) = floor( mean_v.at(last_good_index) ) + _random_shift + (double) ( rand() % _range) / _divisions;
180  mean_v.at(j) = constant_mean;
181  sigma_v.at(j) = mode_sigma;
182  }
183  }
184  }
185  last_good_index = (int)i;
186  }
187  }
188 
189  // Next do extrapolation to the first and end (if needed)
190  // vic: for now we leave this i'm not sure if this really needs
191  // to be tuned until I can make unit test
192  // update: yes this needs work...
193 
194  if (sigma_v.front() > mode_sigma) {
195 
196  int first_index = -1;
197  int second_index = -1;
198 
199  for (size_t i = 0; i < wf.size(); ++i) {
200  if (sigma_v.at(i) < mode_sigma) {
201  if (first_index < 0)
202  first_index = (int)i;
203  else if (second_index < 0) {
204  second_index = (int)i;
205  break;
206  }
207  }
208  }
209 
210  if (first_index < 0 || second_index < 0) {
211  std::cerr << "\033[93m<<" << __FUNCTION__
212  << ">>\033[00m Could not find good pedestal for CDF"
213  << "\n"
214  << "first_index: " << first_index << "\n"
215  << "second_index: " << second_index << "\n"
216  << "If one of these is less than zero, this means there is pulse\n"
217  << "on first sample and baseline never went back down. Returning false here.";
218  return false;
219  }
220 
221  auto diff = fabs(mean_v.at(first_index) - mean_v.at(second_index));
222 
223  if (diff > diff_cutoff) {
224 
225  float slope =
226  (mean_v.at(second_index) - mean_v.at(first_index)) / (float(second_index - first_index));
227 
228  for (int j = 0; j < first_index; ++j) {
229  mean_v.at(j) = mean_v.at(first_index) - slope * (first_index - j);
230  sigma_v.at(j) = _max_sigma;
231  }
232  }
233  else {
234 
235  auto postsample_mean = edge_aware_mean(wf, first_index, first_index + _n_presamples);
236 
237  for (int j = 0; j < first_index; ++j) {
238  //mean_v.at(j) = floor( mean_v.at(second_index) ) + _random_shift + (double) ( rand() % _range) / _divisions;
239  mean_v.at(j) = postsample_mean;
240  sigma_v.at(j) = mode_sigma;
241  }
242  }
243  }
244 
245  if (sigma_v.back() > mode_sigma) {
246 
247  int first_index = -1;
248  int second_index = -1;
249 
250  for (int i = wf.size() - 1; i >= 0; --i) {
251  if (sigma_v.at(i) < mode_sigma) {
252  if (second_index < 0)
253  second_index = (int)i;
254  else if (first_index < 0) {
255  first_index = (int)i;
256  break;
257  }
258  }
259  }
260 
261  if (first_index < 0 || second_index < 0) {
262  std::cerr << "\033[93m<<" << __FUNCTION__
263  << ">>\033[00m Could not find good pedestal for CDF"
264  << "\n"
265  << "first_index: " << first_index << "\n"
266  << "second_index: " << second_index << "\n"
267  << "If one of these is less than zero, this means there is pulse\n"
268  << "on the last sample and baseline never went back down. Returning false here.";
269  return false;
270  }
271 
272  auto diff = fabs(mean_v.at(second_index) - mean_v.at(first_index));
273 
274  if (diff > diff_cutoff) {
275 
276  float slope =
277  (mean_v.at(second_index) - mean_v.at(first_index)) / (float(second_index - first_index));
278  for (int j = second_index + 1; j < int(wf.size()); ++j) {
279  mean_v.at(j) = mean_v.at(second_index) + slope * (j - second_index);
280  sigma_v.at(j) = _max_sigma;
281  }
282  }
283  else {
284 
285  auto presample_mean = edge_aware_mean(wf, second_index - _n_presamples, second_index);
286 
287  for (int j = second_index + 1; j < int(wf.size()); ++j) {
288  //mean_v.at(j) = floor( mean_v.at(first_index) ) + _random_shift + (double) ( rand() % _range) / _divisions;
289  mean_v.at(j) = presample_mean;
290  sigma_v.at(j) = mode_sigma;
291  }
292  }
293  }
294 
295  return true;
296  }
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:43
pmtana::BinnedMaxOccurrence
double BinnedMaxOccurrence(const PedestalMean_t &mean_v, const size_t nbins)
Definition: UtilFunc.cxx:60
pmtana::edge_aware_mean
double edge_aware_mean(const std::vector< short > &wf, int start, int end)
Definition: UtilFunc.cxx:25
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
virtual bool pmtana::PMTPedestalBase::ComputePedestal ( const ::pmtana::Waveform_t &  wf,
pmtana::PedestalMean_t &  mean_v,
pmtana::PedestalSigma_t &  sigma_v 
)
protectedpure virtualinherited

Method to compute pedestal: mean and sigma array should be filled per ADC. The length of each array is guaranteed to be same.

Referenced by pmtana::PMTPedestalBase::Evaluate().

bool pmtana::PMTPedestalBase::Evaluate ( const pmtana::Waveform_t &  wf)
inherited

Method to compute a pedestal.

Definition at line 30 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_mean_v, pmtana::PMTPedestalBase::_sigma_v, and pmtana::PMTPedestalBase::ComputePedestal().

Referenced by pmtana::PedAlgoUB::ComputePedestal(), and pmtana::PulseRecoManager::Reconstruct().

32  {
33  _mean_v.resize(wf.size(), 0);
34  _sigma_v.resize(wf.size(), 0);
35 
36  for (size_t i = 0; i < wf.size(); ++i)
37  _mean_v[i] = _sigma_v[i] = 0;
38 
39  const bool res = ComputePedestal(wf, _mean_v, _sigma_v);
40 
41  if (wf.size() != _mean_v.size())
42  throw OpticalRecoException("Internal error: computed pedestal mean array length changed!");
43  if (wf.size() != _sigma_v.size())
44  throw OpticalRecoException("Internal error: computed pedestal sigma array length changed!");
45 
46  return res;
47  }
pmtana::PMTPedestalBase::_mean_v
pmtana::PedestalMean_t _mean_v
A variable holder for pedestal mean value.
Definition: PMTPedestalBase.h:68
pmtana::PMTPedestalBase::ComputePedestal
virtual bool ComputePedestal(const ::pmtana::Waveform_t &wf, pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)=0
pmtana::PMTPedestalBase::_sigma_v
pmtana::PedestalSigma_t _sigma_v
A variable holder for pedestal standard deviation.
Definition: PMTPedestalBase.h:71
double pmtana::PMTPedestalBase::Mean ( size_t  i) const
inherited

Getter of the pedestal mean value.

Definition at line 50 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_mean_v.

Referenced by pmtana::PedAlgoUB::ComputePedestal(), and pmtana::PulseRecoManager::Reconstruct().

52  {
53  if (i > _mean_v.size()) {
54  std::stringstream ss;
55  ss << "Invalid index: no pedestal mean exist @ " << i;
56  throw OpticalRecoException(ss.str());
57  }
58  return _mean_v[i];
59  }
pmtana::PMTPedestalBase::_mean_v
pmtana::PedestalMean_t _mean_v
A variable holder for pedestal mean value.
Definition: PMTPedestalBase.h:68
const PedestalMean_t & pmtana::PMTPedestalBase::Mean ( ) const
inherited

Getter of the pedestal mean value.

Definition at line 74 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_mean_v.

76  {
77  return _mean_v;
78  }
pmtana::PMTPedestalBase::_mean_v
pmtana::PedestalMean_t _mean_v
A variable holder for pedestal mean value.
Definition: PMTPedestalBase.h:68
const std::string & pmtana::PMTPedestalBase::Name ( ) const
inherited

Name getter.

Definition at line 23 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_name.

25  {
26  return _name;
27  }
pmtana::PMTPedestalBase::_name
std::string _name
Name.
Definition: PMTPedestalBase.h:65
double pmtana::PMTPedestalBase::Sigma ( size_t  i) const
inherited

Getter of the pedestal standard deviation.

Definition at line 62 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_sigma_v.

Referenced by pmtana::PedAlgoUB::ComputePedestal(), and pmtana::PulseRecoManager::Reconstruct().

64  {
65  if (i > _sigma_v.size()) {
66  std::stringstream ss;
67  ss << "Invalid index: no pedestal sigma exist @ " << i;
68  throw OpticalRecoException(ss.str());
69  }
70  return _sigma_v[i];
71  }
pmtana::PMTPedestalBase::_sigma_v
pmtana::PedestalSigma_t _sigma_v
A variable holder for pedestal standard deviation.
Definition: PMTPedestalBase.h:71
const PedestalSigma_t & pmtana::PMTPedestalBase::Sigma ( ) const
inherited

Getter of the pedestal standard deviation.

Definition at line 81 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_sigma_v.

83  {
84  return _sigma_v;
85  }
pmtana::PMTPedestalBase::_sigma_v
pmtana::PedestalSigma_t _sigma_v
A variable holder for pedestal standard deviation.
Definition: PMTPedestalBase.h:71

Member Data Documentation

double pmtana::PedAlgoRollingMean::_diff_adc_count
private

Definition at line 59 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

double pmtana::PedAlgoRollingMean::_diff_threshold
private

Definition at line 58 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

float pmtana::PedAlgoRollingMean::_max_sigma
private

Definition at line 51 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

int pmtana::PedAlgoRollingMean::_n_presamples
private

Definition at line 61 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

float pmtana::PedAlgoRollingMean::_ped_range_max
private

Definition at line 52 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

float pmtana::PedAlgoRollingMean::_ped_range_min
private

Definition at line 53 of file PedAlgoRollingMean.h.

Referenced by PedAlgoRollingMean().

size_t pmtana::PedAlgoRollingMean::_sample_size
private

Definition at line 50 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

double pmtana::PedAlgoRollingMean::_threshold
private

Definition at line 57 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

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