LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
AlgoSlidingWindow.cxx
Go to the documentation of this file.
1 //
3 // AlgoSlidingWindow source
4 //
6 
8 
9 #include "AlgoSlidingWindow.h"
10 
11 namespace pmtana {
12 
13  //*********************************************************************
15  //*********************************************************************
16  {}
17 
18  //*********************************************************************
20  const fhicl::ParameterSet& pset,
21  std::unique_ptr<pmtana::RiseTimeCalculatorBase> risetimecalculator,
22  //AlgoSlidingWindow::AlgoSlidingWindow(const ::fcllite::PSet &pset,
23  const std::string name)
24  : PMTPulseRecoBase(name)
25  //*********************************************************************
26  {
27  _positive = pset.get<bool>("PositivePolarity", true);
28 
29  _adc_thres = pset.get<float>("ADCThreshold");
30 
31  _tail_adc_thres = pset.get<float>("TailADCThreshold", _adc_thres);
32 
33  _end_adc_thres = pset.get<float>("EndADCThreshold");
34 
35  _nsigma = pset.get<float>("NSigmaThreshold");
36 
37  _tail_nsigma = pset.get<float>("TailNSigma", _nsigma);
38 
39  _end_nsigma = pset.get<float>("EndNSigmaThreshold");
40 
41  _verbose = pset.get<bool>("Verbosity");
42 
43  _num_presample = pset.get<size_t>("NumPreSample");
44 
45  _num_postsample = pset.get<size_t>("NumPostSample", 0);
46 
47  _min_width = pset.get<size_t>("MinPulseWidth", 0);
48 
49  _risetime_calc_ptr = std::move(risetimecalculator);
50 
51  Reset();
52  }
53 
54  //***************************************************************
56  //***************************************************************
57  {
59  }
60 
61  //***************************************************************
63  const pmtana::PedestalMean_t& mean_v,
64  const pmtana::PedestalSigma_t& sigma_v)
65  //***************************************************************
66  {
67 
68  bool fire = false;
69 
70  bool in_tail = false;
71 
72  bool in_post = false;
73 
74  double pulse_tail_threshold = 0;
75 
76  double pulse_end_threshold = 0;
77 
78  double pulse_start_baseline = 0;
79 
80  int post_integration = 0;
81 
82  assert(wf.size() == mean_v.size() && wf.size() == sigma_v.size());
83 
84  //double threshold = ( _adc_thres > (_nsigma * _ped_rms) ? _adc_thres : (_nsigma * _ped_rms) );
85 
86  //threshold += _ped_mean;
87 
88  Reset();
89 
90  for (size_t i = 0; i < wf.size(); ++i) {
91 
92  double value = 0.;
93  if (_positive)
94  value = ((double)(wf[i])) - mean_v[i];
95  else
96  value = mean_v[i] - ((double)(wf[i]));
97 
98  float start_threshold = 0.;
99  float tail_threshold = 0.;
100  if (sigma_v[i] * _nsigma < _adc_thres)
101  start_threshold = _adc_thres;
102  else
103  start_threshold = sigma_v[i] * _nsigma;
104 
105  if (sigma_v[i] * _tail_nsigma < _tail_adc_thres)
106  tail_threshold = _tail_adc_thres;
107  else
108  tail_threshold = sigma_v[i] * _tail_nsigma;
109 
110  // End pulse if significantly high peak found (new pulse)
111  if ((!fire || in_tail || in_post) && ((double)value > start_threshold)) {
112 
113  // If there's a pulse, end it
114  if (in_tail) {
115  _pulse.t_end = i - 1;
116 
117  // Register if width is acceptable
118  if ((_pulse.t_end - _pulse.t_start) >= _min_width) {
119  if (_risetime_calc_ptr)
120  _pulse.t_rise = _risetime_calc_ptr->RiseTime(
121  {wf.begin() + _pulse.t_start, wf.begin() + _pulse.t_end},
122  {mean_v.begin() + _pulse.t_start, mean_v.begin() + _pulse.t_end},
123  _positive);
124 
125  _pulse_v.push_back(_pulse);
126  }
127 
129 
130  if (_verbose)
131  std::cout << "\033[93mPulse End\033[00m: "
132  << "baseline: " << mean_v[i] << " ... "
133  << " ... adc above: " << value << " T=" << i << std::endl;
134  }
135 
136  //
137  // Found a new pulse ... try to get a few samples prior to this
138  //
139 
140  pulse_tail_threshold = tail_threshold;
141  pulse_start_baseline = mean_v[i];
142 
143  pulse_end_threshold = 0.;
144  if (sigma_v[i] * _end_nsigma < _end_adc_thres)
145  pulse_end_threshold = _end_adc_thres;
146  else
147  pulse_end_threshold = sigma_v[i] * _end_nsigma;
148 
149  int buffer_num_index = 0;
150  if (_pulse_v.size())
151  buffer_num_index = (int)i - _pulse_v.back().t_end - 1;
152  else
153  buffer_num_index = std::min(_num_presample, i);
154 
155  if (buffer_num_index > (int)_num_presample) buffer_num_index = _num_presample;
156 
157  if (buffer_num_index < 0) {
158  std::cerr << "\033[95m[ERROR]\033[00m Logic error! Negative buffer_num_index..."
159  << std::endl;
160  throw std::exception();
161  }
162 
163  // If there's a pulse, end we where in in_post, end the previous pulse first
164  if (in_post) {
165  // Find were
166  _pulse.t_end = static_cast<int>(i) - buffer_num_index;
167  if (_pulse.t_end > 0) --_pulse.t_end; // leave a gap, if we can
168 
169  // Register if width is acceptable
170  if ((_pulse.t_end - _pulse.t_start) >= _min_width) {
171  if (_risetime_calc_ptr)
172  _pulse.t_rise = _risetime_calc_ptr->RiseTime(
173  {wf.begin() + _pulse.t_start, wf.begin() + _pulse.t_end},
174  {mean_v.begin() + _pulse.t_start, mean_v.begin() + _pulse.t_end},
175  _positive);
176 
177  _pulse_v.push_back(_pulse);
178  }
179 
181 
182  if (_verbose)
183  std::cout << "\033[93mPulse End\033[00m: new pulse starts during in_post: "
184  << "baseline: " << mean_v[i] << " ... "
185  << " ... adc above: " << value << " T=" << i << std::endl;
186  }
187 
188  _pulse.t_start = i - buffer_num_index;
189  _pulse.ped_mean = pulse_start_baseline;
190  _pulse.ped_sigma = sigma_v[i];
191 
192  for (size_t pre_index = _pulse.t_start; pre_index < i; ++pre_index) {
193 
194  double pre_adc = wf[pre_index];
195  if (_positive)
196  pre_adc -= pulse_start_baseline;
197  else
198  pre_adc = pulse_start_baseline - pre_adc;
199 
200  if (pre_adc > 0.) _pulse.area += pre_adc;
201  }
202 
203  if (_verbose)
204  std::cout << "\033[93mPulse Start\033[00m: "
205  << "baseline: " << mean_v[i] << " ... threshold: " << start_threshold
206  << " ... adc above baseline: " << value << " ... pre-adc sum: " << _pulse.area
207  << " T=" << i << std::endl;
208 
209  fire = true;
210  in_tail = false;
211  in_post = false;
212  }
213 
214  if (fire && value < pulse_tail_threshold) {
215  fire = false;
216  in_tail = true;
217  in_post = false;
218  }
219 
220  if ((fire || in_tail || in_post) && _verbose) {
221  std::cout << (fire ? "\033[93mPulsing\033[00m: " : "\033[93mIn-tail\033[00m: ")
222  << "baseline: " << mean_v[i] << " std: " << sigma_v[i]
223  << " ... adc above baseline " << value << " T=" << i << std::endl;
224  }
225 
226  if ((fire || in_tail) && value < pulse_end_threshold) {
227  in_post = true;
228  fire = in_tail = false;
229  post_integration = _num_postsample;
230  }
231 
232  if (in_post && post_integration < 1) {
233  // Found the end of a pulse
234  _pulse.t_end = i - 1;
235 
236  // Register if width is acceptable
237  if ((_pulse.t_end - _pulse.t_start) >= _min_width) {
238  if (_risetime_calc_ptr)
239  _pulse.t_rise = _risetime_calc_ptr->RiseTime(
240  {wf.begin() + _pulse.t_start, wf.begin() + _pulse.t_end},
241  {mean_v.begin() + _pulse.t_start, mean_v.begin() + _pulse.t_end},
242  _positive);
243 
244  _pulse_v.push_back(_pulse);
245  }
246 
247  if (_verbose)
248  std::cout << "\033[93mPulse End\033[00m: "
249  << "baseline: " << mean_v[i] << " ... adc: " << value << " T=" << i
250  << " ... area sum " << _pulse.area << std::endl;
251 
253 
254  fire = false;
255  in_tail = false;
256  in_post = false;
257  }
258 
259  if (fire || in_tail || in_post) {
260 
261  //_pulse.area += ((double)value - (double)mean_v[i]);
262  _pulse.area += value;
263 
264  if (_pulse.peak < value) {
265 
266  // Found a new maximum
267  _pulse.peak = value;
268 
269  _pulse.t_max = i;
270  }
271 
272  if (in_post) --post_integration;
273  }
274  }
275 
276  if (fire || in_tail || in_post) {
277 
278  // Take care of a pulse that did not finish within the readout window.
279 
280  fire = false;
281  in_tail = false;
282 
283  _pulse.t_end = wf.size() - 1;
284 
285  // Register if width is acceptable
286  if ((_pulse.t_end - _pulse.t_start) >= _min_width) {
287  if (_risetime_calc_ptr)
288  _pulse.t_rise = _risetime_calc_ptr->RiseTime(
289  {wf.begin() + _pulse.t_start, wf.begin() + _pulse.t_end},
290  {mean_v.begin() + _pulse.t_start, mean_v.begin() + _pulse.t_end},
291  _positive);
292  _pulse_v.push_back(_pulse);
293  }
294 
296  }
297 
298  return true;
299  }
300 
301 }
std::vector< double > PedestalSigma_t
virtual void Reset()
A method to be called event-wise to reset parameters.
float _nsigma
A variable holder for a multiplicative factor for the pedestal standard deviation to define the thres...
AlgoSlidingWindow(const std::string name="SlidingWindow")
Default constructor.
bool RecoPulse(const pmtana::Waveform_t &, const pmtana::PedestalMean_t &, const pmtana::PedestalSigma_t &)
Implementation of AlgoSlidingWindow::reco() method.
pulse_param _pulse
A subject pulse_param object to be filled with the last reconstructed pulse parameters.
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::vector< short > Waveform_t
double value
Definition: spectrum.C:18
std::unique_ptr< pmtana::RiseTimeCalculatorBase > _risetime_calc_ptr
Tool for rise time calculation.
bool _positive
A boolean to set waveform positive/negative polarity.
Class definition file of AlgoSlidingWindow.
void Reset()
Implementation of AlgoSlidingWindow::reset() method.
float _adc_thres
A variable holder for a user-defined absolute ADC threshold value.
std::vector< double > PedestalMean_t
size_t _min_width
A variable holder to ensure the minimum pulse width.
pulse_param_array _pulse_v
A container array of pulse_param struct objects to store (possibly multiple) reconstructed pulse(s)...
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33