LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PedAlgoRmsSlider.cxx
Go to the documentation of this file.
1 //
3 // PedAlgoRmsSlider source
4 //
6 
7 #include "PedAlgoRmsSlider.h"
8 #include "OpticalRecoException.h"
9 #include "UtilFunc.h"
10 #include "fhiclcpp/ParameterSet.h"
11 
12 #include <fstream>
13 #include <iostream>
14 #include <numeric>
15 
16 namespace pmtana {
17 
18  //*****************************************************************
19  PedAlgoRmsSlider::PedAlgoRmsSlider(const std::string name) : PMTPedestalBase(name)
20  //*****************************************************************
21  {
22  srand(static_cast<unsigned int>(time(0)));
23  }
24 
25  //**************************************************************************
26  PedAlgoRmsSlider::PedAlgoRmsSlider(const fhicl::ParameterSet& pset, const std::string name)
27  : PMTPedestalBase(name)
28  //############################################################
29  {
30 
31  _sample_size = pset.get<size_t>("SampleSize", 7);
32  _threshold = pset.get<double>("Threshold", 0.6);
33  _max_sigma = pset.get<float>("MaxSigma", 0.5);
34  _ped_range_max = pset.get<float>("PedRangeMax", 2150);
35  _ped_range_min = pset.get<float>("PedRangeMin", 100);
36  _num_presample = pset.get<int>("NumPreSample", 0);
37  _num_postsample = pset.get<int>("NumPostSample", 0);
38  _verbose = pset.get<bool>("Verbose", true);
39  _n_wf_to_csvfile = pset.get<int>("NWaveformsToFile", 12);
40 
41  if (_n_wf_to_csvfile > 0) {
42  _csvfile.open("wf_pedalgormsslider.csv", std::ofstream::out | std::ofstream::trunc);
43  _csvfile << "n,time,wf,wf_ped_mean,wf_ped_rms" << std::endl;
44  }
45  }
46 
47  //*******************************************
49  //*******************************************
50  {
51  std::cout << "PedAlgoRmsSlider setting:"
52  << "\n\t SampleSize: " << _sample_size
53  << "\n\t Threshold: " << _threshold << "\n\t Verbose: " << _verbose
54  << "\n\t NWaveformsToFile: " << _n_wf_to_csvfile << std::endl;
55  }
56 
57  //****************************************************************************
58  double PedAlgoRmsSlider::CalcMean(const std::vector<double>& wf, size_t start, size_t nsample)
59  //****************************************************************************
60  {
61  if (!nsample) nsample = wf.size();
62  if (start > wf.size() || (start + nsample) > wf.size())
63  throw OpticalRecoException("Invalid start/end index!");
64 
65  double sum = std::accumulate(wf.begin() + start, wf.begin() + start + nsample, 0.0);
66 
67  sum /= ((double)nsample);
68 
69  return sum;
70  }
71 
72  //****************************************************************************
73  double PedAlgoRmsSlider::CalcStd(const std::vector<double>& wf,
74  const double ped_mean,
75  size_t start,
76  size_t nsample)
77  //****************************************************************************
78  {
79  if (!nsample) nsample = wf.size();
80  if (start > wf.size() || (start + nsample) > wf.size())
81  throw OpticalRecoException("Invalid start/end index!");
82 
83  double sigma = 0;
84 
85  for (size_t index = start; index < (start + nsample); ++index) {
86  sigma += pow((wf[index] - ped_mean), 2);
87  }
88 
89  sigma = sqrt(sigma / ((double)(nsample)));
90 
91  return sigma;
92  }
93 
94  //****************************************************************************
96  pmtana::PedestalMean_t& mean_v,
97  pmtana::PedestalSigma_t& sigma_v)
98  //****************************************************************************
99  {
100 
101  if (_verbose) this->PrintInfo();
102 
103  if (wf.size() <= (_sample_size * 2)) return false;
104 
105  // Prepare output
106  mean_v.resize(wf.size(), 0);
107  sigma_v.resize(wf.size(), 0);
108 
109  // **********
110  // To start, set the pedestal equal to
111  // the wf itself
112  // **********
113 
114  pmtana::PedestalMean_t mean_temp_v;
115  mean_temp_v.resize(wf.size(), 0);
116 
117  for (size_t i = 0; i < wf.size(); ++i) {
118  mean_temp_v[i] = wf[i];
119  sigma_v[i] = 0;
120  }
121 
122  // **********
123  // Now look for rms variations
124  // and change the mean and rms accordingly
125  // **********
126  int last_good_index = -1;
127  double local_mean, local_rms;
128  std::vector<double> local_mean_v(wf.size(), -1.);
129  std::vector<double> local_sigma_v(wf.size(), -1.);
130 
131  for (size_t i = 0; i < wf.size() - _sample_size; i++) {
132 
133  local_mean = mean(wf, i, _sample_size);
134  local_rms = std(wf, local_mean, i, _sample_size);
135 
136  if (_verbose)
137  std::cout << "\033[93mPedAlgoRmsSlider\033[00m: i " << i << " local_mean: " << local_mean
138  << " local_rms: " << local_rms << std::endl;
139 
140  if (local_rms < _threshold) {
141 
142  local_mean_v[i] = local_mean;
143  local_sigma_v[i] = local_rms;
144 
145  if (_verbose)
146  std::cout << "\033[93mBelow threshold\033[00m: "
147  << "at i " << i << " last good index was: " << last_good_index << std::endl;
148  }
149  }
150 
151  // find the gaps (regions to be interpolated
152  last_good_index = -1;
153  std::vector<bool> ped_interapolated(wf.size(), false);
154  for (size_t i = 0; i < wf.size() - _sample_size; i++) {
155 
156  if (local_mean_v[i] > -0.1) {
157  // good pedestal!
158 
159  if ((last_good_index + 1) < (int)i) {
160  // finished the gap. try interpolation
161  // 0) find where to start/end interpolation
162  int start_tick = last_good_index;
163  int end_tick = i;
164  int start_bound = std::max(last_good_index - _num_presample, 0);
165  int end_bound = std::min(i + _num_postsample, (int)(wf.size()) - _sample_size);
166  for (int j = start_tick; j >= start_bound; --j) {
167  if (local_mean_v[j] < 0) continue;
168  start_tick = j;
169  }
170  for (int j = end_tick; j <= end_bound; ++j) {
171  if (local_mean_v[j] < 0) continue;
172  end_tick = j;
173  }
174 
175  //this should become generic interpolation function, for now lets leave.
176  float slope =
177  (local_mean_v[end_tick] - local_mean_v[start_tick]) / (float(end_tick - start_tick));
178 
179  for (int j = start_tick + 1; j < end_tick; ++j) {
180  mean_temp_v[j] = slope * (float(j - start_tick)) + local_mean_v[start_tick];
181  // for sigma, put either the sigma in the region before the pulse or
182  // after the pulse, depending on which one if != 0. If both are !=0 put the one after
183  // the pulse (just default), if both are zero then put zero
184  sigma_v[j] = (local_sigma_v[end_tick] != 0 ?
185  local_sigma_v[end_tick] :
186  local_sigma_v[start_tick]); // todo: fluctuate baseline
187  ped_interapolated[j] = true;
188  }
189  }
190 
191  last_good_index = i;
192  }
193  }
194 
195  /*
196 
197  for (size_t i = 0; i < wf.size() - _sample_size; i++) {
198 
199  local_mean = mean(wf, i, _sample_size);
200  local_rms = std(wf, local_mean, i, _sample_size);
201 
202  if(_verbose) std::cout << "\033[93mPedAlgoRmsSlider\033[00m: i " << i << " local_mean: " << local_mean << " local_rms: " << local_rms << std::endl;
203 
204  if (local_rms < _threshold) {
205 
206  local_mean_v[i] = local_mean;
207  local_sigma_v[i] = local_sigma;
208 
209  if(_verbose)
210  std::cout << "\033[93mBelow threshold\033[00m: "
211  << "at i " << i
212  << " last good index was: " << last_good_index
213  << std::endl;
214 
215  if(last_good_index<0) {
216  last_good_index = (int)i;
217  last_local_mean = local_mean;
218  last_local_rms = local_rms;
219  continue;
220  }
221 
222 
223  if( ( last_good_index + 1 ) < (int)i ) {
224 
225  //this should become generic interpolation function, for now lets leave.
226  float slope = (local_mean - last_local_mean) / (float(i - last_good_index));
227 
228  for(size_t j = last_good_index + 1; j < i && j < wf.size(); ++j) {
229  mean_temp_v.at(j) = slope * ( float(j - last_good_index) ) + mean_temp_v.at(last_good_index);
230  // for sigma, put either the sigma in the region before the pulse or
231  // after the pulse, depending on which one if != 0. If both are !=0 put the one after
232  // the pulse (just default), if both are zero then put zero
233  sigma_v.at(j) = (local_rms != 0 ? local_rms : last_local_rms); // todo: fluctuate baseline
234  ped_interapolated.at(j) = true;
235  }
236  }
237 
238  // record this mean & rms as good mean value
239  last_good_index = i;
240  last_local_mean = local_mean;
241  last_local_rms = local_rms;
242  // if _num_postsample is specified, go back in time to look for it
243  if(_num_postsample >0 && (i>_num_postsample)) {
244  int loop_start = std::max(((int)i) - _num_postsample, 0);
245  for(int j=loop_start; j>=0; --j) {
246  if(local_mean_v[j] <0) continue;
247  last_good_index = j;
248  last_local_mean = local_mean_v[j];
249  last_local_rms = local_sigma_v[j];
250  break;
251  }
252  }
253  }
254 
255 
256  }
257  */
258 
259  // **********
260  // Now look at special cases, if wf starts or
261  // ends with a pulse
262  // **********
263 
264  // At start
265 
266  bool end_found = false;
267 
268  local_mean = mean(wf, 0, _sample_size);
269  local_rms = std(wf, local_mean, 0, _sample_size);
270 
271  if (local_rms >= _threshold) {
272 
273  for (size_t i = 1; i < wf.size() - _sample_size; i++) {
274 
275  local_mean = mean(wf, i, _sample_size);
276  local_rms = std(wf, local_mean, i, _sample_size);
277 
278  if (local_rms < _threshold) {
279 
280  end_found = true;
281 
282  for (size_t j = 0; j < i; j++) {
283  mean_temp_v[j] = local_mean;
284  sigma_v[j] = local_rms;
285  ped_interapolated[j] = true;
286  }
287  break;
288  }
289  }
290 
291  if (!end_found) {
292  std::cerr << "\033[93m<<" << __FUNCTION__
293  << ">>\033[00m Could not find good pedestal for CDF"
294  << "There is pulse on first sample and baseline never went back down. Returning "
295  "false here.";
296  return false;
297  }
298  }
299 
300  // At end
301 
302  bool start_found = false;
303 
304  local_mean = mean(wf, wf.size() - 1 - _sample_size, _sample_size);
305  local_rms = std(wf, local_mean, wf.size() - 1 - _sample_size, _sample_size);
306 
307  if (local_rms >= _threshold) {
308 
309  size_t i = wf.size() - 1 - _sample_size;
310  while (i-- > 0) {
311  local_mean = mean(wf, i, _sample_size);
312  local_rms = std(wf, local_mean, i, _sample_size);
313 
314  if (local_rms < _threshold) {
315 
316  start_found = true;
317 
318  for (size_t j = wf.size() - 1; j > i; j--) {
319  mean_temp_v[j] = local_mean;
320  sigma_v[j] = local_rms;
321  ped_interapolated[j] = true;
322  }
323  break;
324  }
325  }
326 
327  if (!start_found) {
328  std::cerr << "\033[93m<<" << __FUNCTION__
329  << ">>\033[00m Could not find good pedestal for CDF"
330  << "There is pulse on last sample and baseline never went back down. Returning "
331  "false here.";
332  return false;
333  }
334  }
335 
336  // **********
337  // Now smooth it to estimate the final pedestal
338  // **********
339 
340  const size_t window_size = _sample_size * 2;
341 
342  // middle mean
343  for (size_t i = 0; i < mean_temp_v.size(); ++i) {
344 
345  if (i < _sample_size || i >= (wf.size() - _sample_size)) continue;
346 
347  mean_v[i] = this->CalcMean(mean_temp_v, i - _sample_size, window_size);
348  if (!ped_interapolated[i]) {
349  sigma_v[i] = this->CalcStd(mean_temp_v, mean_v[i], i - _sample_size, window_size);
350  }
351  }
352 
353  // front mean
354  for (size_t i = 0; i < _sample_size; ++i) {
355 
356  mean_v[i] = mean_v[_sample_size];
357  if (!ped_interapolated[i]) { sigma_v[i] = sigma_v[_sample_size]; }
358  }
359 
360  // tail mean
361  for (size_t i = (mean_temp_v.size() - _sample_size); i < mean_temp_v.size(); ++i) {
362 
363  mean_v[i] = mean_v[wf.size() - _sample_size - 1];
364  if (!ped_interapolated[i]) { sigma_v[i] = sigma_v[wf.size() - _sample_size - 1]; }
365  }
366 
367  // Save to file
368  if (_wf_saved + 1 <= _n_wf_to_csvfile) {
369  _wf_saved++;
370  for (size_t i = 0; i < wf.size(); i++) {
371  _csvfile << _wf_saved - 1 << "," << i << "," << wf[i] << "," << mean_v[i] << ","
372  << sigma_v[i] << std::endl;
373  }
374  }
375 
376  bool is_sane = this->CheckSanity(mean_v, sigma_v);
377 
378  return is_sane;
379  }
380 
381  //*******************************************
383  pmtana::PedestalSigma_t& sigma_v)
384  //*******************************************
385  {
386 
387  float best_sigma = 1.1e9;
388  size_t best_sigma_index = 0;
389  size_t num_good_adc = 0;
390 
391  for (size_t i = 0; i < sigma_v.size(); ++i) {
392  // Only consider adcs which mean is in the allowed range
393  auto const& mean = mean_v[i];
394 
395  if (mean < _ped_range_min || mean > _ped_range_max) continue;
396 
397  auto const& sigma = sigma_v[i];
398  if (sigma < best_sigma) {
399  best_sigma = sigma;
400  best_sigma_index = i;
401  }
402 
403  if (sigma < _max_sigma) num_good_adc += 1;
404  }
405 
406  if (num_good_adc < 1) {
407  std::cerr << "\033[93m<<" << __FUNCTION__
408  << ">>\033[00m Could not find good pedestal at all..." << std::endl;
409  return false;
410  }
411 
412  // If not enough # of good mean indices, use the best guess within this waveform
413  if (best_sigma > _max_sigma || num_good_adc < 3) {
414 
415  if (_verbose) {
416  std::cout << "\033[93mPedAlgoRmsSlider\033[00m: Not enough number of good mean indices."
417  << "Using the best guess within this waveform." << std::endl;
418  }
419 
420  for (size_t i = 0; i < mean_v.size(); ++i) {
421  mean_v[i] = mean_v[best_sigma_index];
422  sigma_v[i] = sigma_v[best_sigma_index];
423  }
424  }
425 
426  return true;
427  }
428 }
Class def header for exception classes in OpticalDetector package.
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:43
void PrintInfo()
Print settings.
bool CheckSanity(pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)
Checks the sanity of the estimated pedestal, returns false if not sane.
std::vector< double > PedestalSigma_t
float _ped_range_min
Min value of adc to consider adc as &#39;sane&#39;.
int _num_postsample
number of ADCs to sample after the gap
float _max_sigma
Max sigma to consider adc as &#39;sane&#39;.
PedAlgoRmsSlider(const std::string name="PedRmsSlider")
Default constructor.
int _num_presample
number of ADCs to sample before the gap
double CalcStd(const std::vector< double > &wf, const double ped_mean, size_t start, size_t nsample)
Returns the std of the elements of the vector from start to start+nsample.
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...
double _threshold
Threshold applied to local rms to claim a pulse.
float _ped_range_max
Max value of adc to consider adc as &#39;sane&#39;.
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::vector< short > Waveform_t
int _n_wf_to_csvfile
If greater than zero saves firsts waveforms with pedestal to csv file.
size_t _sample_size
How many samples are used to calculate local rms and mean.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
Class definition file of PedAlgoRmsSlider.
double CalcMean(const std::vector< double > &wf, size_t start, size_t nsample)
Returns the mean of the elements of the vector from start to start+nsample.
Double_t sum
Definition: plot.C:31
std::vector< double > PedestalMean_t
bool _verbose
For debugging.