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