CFDHitFinderAlg.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <map>
4 #include <iostream>
5 #include <vector>
6 #include <cstdio>
7 #include <cstddef> // for size_t
8 #include <cmath> // for fabs
9 
10 #include <TH1D.h>
11 #include <TF1.h>
12 #include <TMath.h>
13 #include <TH1D.h>
14 #include <TGraph.h>
15 #include <TFitResult.h>
16 #include <TFitResultPtr.h>
17 #include <TVirtualFitter.h>
18 
19 #include "KalmanFilterAlg.h"
20 
21 namespace beamlinereco {
22  template <class T>
23  struct hit_t {
24  unsigned int DigitizerChannel;
25  uint32_t Timestamp;
36  };
37 
38  enum CFDParams {
45  kDiscriminationThreshold, // percentage for CFD, number of noise bands for LE
59  };
60 
61  std::vector<double> sg_smooth(const std::vector<double> &v, const int w, const int deg);
62  std::vector<double> sg_derivative(const std::vector<double> &v, const int w, const int deg, const double h = 1.0);
63 
64  class SGSmoothing {
65  public:
66  static void Smooth(size_t nsamples, double* in, double* out, const int w, const int deg);
67  static void Derivative(size_t nsamples, double* in, double* out, const int w, const int deg, const double h = 1.0);
68  };
69 
70  template <class T>
71  class CFDHitFinder {
72  public:
74  _HitCollection.clear();
75  _RawHitLogicNanosec.clear();
76  _PedestalInADC = 0;
77  _NoiseSigmaInADC = 0;
78  _DiscriminationThresholdInADC = 0;
79  _WaveformADCNanosec.clear();
80  _ChannelNumber = 993;
81  _Timestamp = 0;
82  _NonfilterWaveformADCNanosec.clear();
83  _CFDParamSet.clear();
84  }
85  virtual ~CFDHitFinder() = default;
86 
87 
88  void SetWaveform(std::vector<uint16_t>& waveform, unsigned int channelNo, uint32_t timestamp);
89  void SetParams(const std::map<CFDParams, double>& paramSet);
90  void SetDiscriminationType(const std::string type);
91  void SetInterpolationType(const std::string type);
92  void SetFilterType(const std::string type);
93 
94  inline const std::string GetDiscriminationType() const { return _DiscriminationType; }
95  inline const std::string GetInterpolationType() const { return _InterpolationType; }
96  inline const std::string GetFilterType() const { return _FilterType; }
97  inline const std::map<T, hit_t<T> >& GetHitCollection() const { return _HitCollection; }
98  inline const T& GetPedestal() const { return _PedestalInADC; }
99  inline const T& GetNoiseSigma() const { return _NoiseSigmaInADC; }
100  inline const T& GetDiscriminationThreshold() const { return _DiscriminationThresholdInADC; }
101 
102  virtual void Go();
103 
104  private:
105 
106  void SetChannel(unsigned int channelNo);
107  void SetTimestamp(uint32_t timestamp);
108  void SetParam(CFDParams param, T value);
109 
110  private:
111  unsigned int _ChannelNumber;
112  uint32_t _Timestamp;
116  std::vector<T> _WaveformADCNanosec;
118  std::vector<bool> _RawHitLogicNanosec;
122 
123  std::map<CFDParams, T> _CFDParamSet;
124  std::map<T, hit_t<T> > _HitCollection;
125 
126  private:
127  virtual void FindPedestal();
128  virtual void FindRawHitLogic();
129  virtual void FindCFDHits();
130  virtual bool BackwardFindingOfHitStart(size_t hitPeakTimeIndex, T hitPeakValue, T& hitStartTimeIndex, T& hitRiseTimeInIndex);
131  virtual bool ForwardFindingOfHitFallTime(size_t hitPeakTimeIndex, T& hitFallTimeInIndex);
132  virtual T IntegrateWaveformInADC(size_t hitStartTimeIndex, size_t hitEndTimeIndex);
133  virtual T FindHitWidth(size_t hitPeakTimeIndex, T hitPeakValueADC);
134  virtual T FindPeakValue(size_t hitPeakTimeIndex, T &hitPeakValue);
135 
136  virtual inline void Reset() {
137  _HitCollection.clear();
138  _RawHitLogicNanosec.clear();
139  _PedestalInADC = 0;
140  _NoiseSigmaInADC = 0;
141  _DiscriminationThresholdInADC = 0;
142  _WaveformADCNanosec.clear();
143  _ChannelNumber = 993;
144  _Timestamp = 0;
145  _NonfilterWaveformADCNanosec.clear();
146  }
147  };
148 }
149 
150 template struct beamlinereco::hit_t<double>;
152 
153 template<class T>
154 void beamlinereco::CFDHitFinder<T>::SetWaveform(std::vector<uint16_t>& waveform, unsigned int channelNo, uint32_t timestamp) {
155  Reset();
156 
157  SetChannel(channelNo);
158  SetTimestamp(timestamp);
159 
160  double polaritySignFactor = _CFDParamSet[kIsWaveformNegativePolarity] ? 1. : -1.;
161 
162  T* nonfilterwaveform = (T*) malloc(_CFDParamSet[kNSamplingPoints] * sizeof(T));
163  T* filterwaveform = (T*) malloc(_CFDParamSet[kNSamplingPoints] * sizeof(T));
164 
165  for (auto itr = waveform.begin(); itr != waveform.end(); itr++) {
166  _NonfilterWaveformADCNanosec.push_back((T)(*itr) * polaritySignFactor);
167  nonfilterwaveform[itr - waveform.begin()] = (T)(*itr) * polaritySignFactor;
168  }
169 
170  if (_FilterType == "GS" || _FilterType == "gs") { // Savitzky-Golay smoothing
172  nonfilterwaveform,
173  filterwaveform,
174  _CFDParamSet[kGSFilterWindow],
175  _CFDParamSet[kGSFilterDegree]);
176  } else if (_FilterType == "Kalman" || _FilterType == "kalman" || _FilterType == "KALMAN") { // Kalman filter
177  float q = _CFDParamSet[kKalmanFilterProcessNoiseCovariance];
178  float r = _CFDParamSet[kKalmanFilterMeasurementNoiseCovariance];
179  float P = .05;
180  float K = _CFDParamSet[kKalmanFilterGain];
181  Kalman filter = Kalman(q,r,P,K);
182  for (unsigned int i=0; i<_CFDParamSet[kNSamplingPoints]; i++) {
183  filterwaveform[i] = filter.GetFilteredValue(nonfilterwaveform[i]);
184  }
185  } else if (_FilterType == "" || _FilterType == "None" || _FilterType == "none" || _FilterType == "NONE") {
186  // no smoothing
187  for (auto itr = waveform.begin(); itr != waveform.end(); itr++) {
188  filterwaveform[itr - waveform.begin()] = (T)(*itr) * polaritySignFactor;
189  }
190  } else {
191  // unknown -- throw exception
193  << "Unknown raw hit finder filter type: " << _FilterType << "." << std::endl;
194  }
195 
196  for (size_t i = 0; i < _CFDParamSet[kNSamplingPoints]; i++) {
197  _WaveformADCNanosec.push_back(*(filterwaveform + i));
198  }
199 
200  free(nonfilterwaveform);
201  free(filterwaveform);
202 }
203 
204 template<class T>
206  _CFDParamSet[param] = value;
207 }
208 
209 template<class T>
210 void beamlinereco::CFDHitFinder<T>::SetParams(const std::map<CFDParams, double>& paramSet) {
211  for (auto & itr : paramSet) {
212  SetParam(itr.first, itr.second);
213  }
214 }
215 
216 template<class T>
217 void beamlinereco::CFDHitFinder<T>::SetChannel(unsigned int channelNo) {
218  _ChannelNumber = channelNo;
219 }
220 
221 template<class T>
223  _Timestamp = ts;
224 }
225 
226 template<class T>
228  _InterpolationType = type;
229 }
230 
231 template<class T>
233  _DiscriminationType = type;
234 }
235 
236 template<class T>
238  _FilterType = type;
239 }
240 
241 template<class T>
243 /* 2nd way
244  TH1D* pedHist = new TH1D("", "", 10000, 0, 5000);
245  for (size_t i = 700; i < 1024; i++) {
246  pedHist->Fill(_NonfilterWaveformADCNanosec.at(i));
247  }
248 
249  _PedestalInADC = (T) pedHist->GetMean();
250  _NoiseSigmaInADC = (T) pedHist->GetStdDev();
251 
252  delete (pedHist);
253 */
254 
255 /* 3rd way*/
256  T NoiseBand = 1E9;
257  T OldPedestal = 1E9;
258  T Pedestal = 0;
259 
260  while (TMath::Abs(Pedestal - OldPedestal) >= 2) {
261  T average = 0;
262  T stdDev = 0;
263  size_t countSamples = 0;
264  for (typename std::vector<T>::iterator itr = _NonfilterWaveformADCNanosec.begin(); itr != _NonfilterWaveformADCNanosec.end(); itr++) {
265  if ((*itr <= Pedestal + NoiseBand) && (*itr >= Pedestal - NoiseBand)) {
266  average += *itr;
267  countSamples++;
268  }
269  }
270  average = average / (T)countSamples;
271  for (typename std::vector<T>::iterator itr = _NonfilterWaveformADCNanosec.begin(); itr != _NonfilterWaveformADCNanosec.end(); itr++) {
272  if ((*itr <= Pedestal + NoiseBand) && (*itr >= Pedestal - NoiseBand)) {
273  stdDev += TMath::Power((*itr) - average, 2);
274  }
275  }
276 
277  stdDev = TMath::Sqrt(stdDev / (T)countSamples);
278  OldPedestal = Pedestal;
279  Pedestal = average;
280  NoiseBand = stdDev;
281  }
282 
283  _PedestalInADC = Pedestal;
284  _NoiseSigmaInADC = NoiseBand;
285 /**/
286 }
287 
288 template<class T>
290  T rawHitThreshold = (_PedestalInADC) - ((_NoiseSigmaInADC) * _CFDParamSet[kRawHitFinderThresholdInNoiseSigma]);
291  for (size_t i = 0; i < _WaveformADCNanosec.size(); i++) {
292  _RawHitLogicNanosec.push_back(_WaveformADCNanosec.at(i) < rawHitThreshold);
293  }
294 
295  for (size_t i = 1; i < _RawHitLogicNanosec.size() - 1; i++) {
296  if ((not _RawHitLogicNanosec.at(i - 1)) and _RawHitLogicNanosec.at(i)) {
297  size_t countDuration = 0;
298  for (size_t j = i; j < _RawHitLogicNanosec.size() - 1; j++) {
299  if (_RawHitLogicNanosec.at(j)) countDuration++;
300  if (not(_RawHitLogicNanosec.at(j + 1))) break;
301  }
302 
303  if (countDuration < _CFDParamSet[kShortRawHitIgnoringDurationInTicks]) {
304  for (size_t j = i; j < i + countDuration; j++) {
305  _RawHitLogicNanosec.at(j) = false;
306  }
307  }
308  }
309  }
310 
311  for (size_t i = 1; i < _RawHitLogicNanosec.size() - 1; i++) {
312  if (_RawHitLogicNanosec.at(i - 1) && not(_RawHitLogicNanosec.at(i))) {
313  size_t countDuration = 0;
314  for (size_t j = i; j < _RawHitLogicNanosec.size() - 1; j++) {
315  if (not(_RawHitLogicNanosec.at(j))) countDuration++;
316  if (_RawHitLogicNanosec.at(j + 1)) break;
317  }
318 
319  if (countDuration < _CFDParamSet[kConsecutiveHitSeperationDurationInTicks]) {
320  for (size_t j = i; j < i + countDuration; j++) {
321  _RawHitLogicNanosec.at(j) = true;
322  }
323  }
324  }
325  }
326 
327 }
328 
329 template<class T>
331  std::vector<T> hitPeakValueVector; hitPeakValueVector.clear();
332  std::vector<size_t> hitPeakTimeIndexVector; hitPeakTimeIndexVector.clear();
333  std::vector<T> hitPeakTimeInNanoSecVector; hitPeakTimeInNanoSecVector.clear();
334  std::vector<T> hitStartTimeIndexVector; hitStartTimeIndexVector.clear();
335  std::vector<T> hitRiseTimeInIndexVector; hitRiseTimeInIndexVector.clear();
336  std::vector<T> hitFallTimeInIndexVector; hitFallTimeInIndexVector.clear();
337  std::vector<T> hitWidthValueVector; hitWidthValueVector.clear();
338  std::vector<bool> isHitContainedVector; isHitContainedVector.clear();
339 
340  bool logicFallingPart = true;
341  for (size_t i = 1; i < _RawHitLogicNanosec.size() - 1; i++) {
342  T hitPeakValue = 1E10;
343  size_t hitPeakTimeIndex = 0;
344  if (_RawHitLogicNanosec.at(i) && !(_RawHitLogicNanosec.at(i - 1)) && logicFallingPart) {
345  logicFallingPart = false;
346  for (size_t j = i; j < _RawHitLogicNanosec.size() - 1; j++) {
347  hitPeakValue = hitPeakValue > _WaveformADCNanosec.at(j) ? (T) _WaveformADCNanosec.at(j) : hitPeakValue;
348  hitPeakTimeIndex = hitPeakValue == _WaveformADCNanosec.at(j) ? j : hitPeakTimeIndex;
349  if (((_RawHitLogicNanosec.at(j) && !(_RawHitLogicNanosec.at(j + 1))) || ((_RawHitLogicNanosec.size() - j) < _CFDParamSet[kRawHitFinderTicksFromEnd])) && !logicFallingPart) {
350  logicFallingPart = true;
351  T hitStartTimeIndex = (T)hitPeakTimeIndex;
352  T hitRiseTimeInIndex = (T)hitPeakTimeIndex;
353  T hitFallTimeInIndex = (T)hitPeakTimeIndex;
354  T hitPeakTimeInNanoSec = FindPeakValue(hitPeakTimeIndex, hitPeakValue);
355  T hitWidthInTicks = FindHitWidth(hitPeakTimeIndex, hitPeakValue);
356  bool containedFromRisingEdge = BackwardFindingOfHitStart(hitPeakTimeIndex, hitPeakValue, hitStartTimeIndex, hitRiseTimeInIndex);
357  bool containedFromFallingEdge = ForwardFindingOfHitFallTime(hitPeakTimeIndex, hitFallTimeInIndex);
358 
359  hitPeakValueVector.push_back(hitPeakValue);
360  hitPeakTimeIndexVector.push_back(hitPeakTimeIndex);
361  hitPeakTimeInNanoSecVector.push_back(hitPeakTimeInNanoSec);
362  hitStartTimeIndexVector.push_back(hitStartTimeIndex);
363  hitRiseTimeInIndexVector.push_back(hitRiseTimeInIndex);
364  hitFallTimeInIndexVector.push_back(hitFallTimeInIndex);
365  hitWidthValueVector.push_back(hitWidthInTicks);
366  isHitContainedVector.push_back(containedFromRisingEdge && containedFromFallingEdge);
367  break;
368  }
369  }
370  }
371  }
372 
373  for (size_t i = 0; i < hitStartTimeIndexVector.size(); i++) {
374  hit_t<double> aHit;
375  aHit.Timestamp = _Timestamp;
376  aHit.DigitizerChannel = _ChannelNumber;
377  aHit.TStartInNanoSec = hitStartTimeIndexVector.at(i) * _CFDParamSet[kTimeSamplingInterval];
378  aHit.TPeakInNanoSec = hitPeakTimeInNanoSecVector.at(i) * _CFDParamSet[kTimeSamplingInterval];
379  aHit.AmplitudeInMiliVolt = hitPeakValueVector.at(i) * (_CFDParamSet[kADCDynamicRange] / (TMath::Power(2, _CFDParamSet[kADCNBits]) - 1));
380  aHit.AmplitudeInADC = hitPeakValueVector.at(i);
381  aHit.WidthInNanoSec = hitWidthValueVector.at(i) * _CFDParamSet[kTimeSamplingInterval];
382  aHit.RiseTimeInNanoSec = hitRiseTimeInIndexVector.at(i) * _CFDParamSet[kTimeSamplingInterval];
383  aHit.FallTimeInNanoSec = hitFallTimeInIndexVector.at(i) * _CFDParamSet[kTimeSamplingInterval];
384  aHit.IsContained = isHitContainedVector.at(i);
385  aHit.IntegratedChargeInADCTimeTicks = IntegrateWaveformInADC(hitPeakTimeIndexVector.at(i) - hitRiseTimeInIndexVector.at(i),
386  hitPeakTimeIndexVector.at(i) + hitFallTimeInIndexVector.at(i));
388 
389  _HitCollection[hitStartTimeIndexVector.at(i)] = aHit;
390  }
391 }
392 
393 template<class T>
394 T beamlinereco::CFDHitFinder<T>::FindPeakValue(size_t hitPeakTimeIndex, T &hitPeakValue) {
395  // Fit a third-order polynomial
396  // Translate x to near zero to help out the fitter
397  T x[7] = { -3, -2, -1, 0, 1, 2, 3 };
398  T y[7] = {
399  _WaveformADCNanosec.at(hitPeakTimeIndex - 3),
400  _WaveformADCNanosec.at(hitPeakTimeIndex - 2),
401  _WaveformADCNanosec.at(hitPeakTimeIndex - 1),
402  _WaveformADCNanosec.at(hitPeakTimeIndex),
403  _WaveformADCNanosec.at(hitPeakTimeIndex + 1),
404  _WaveformADCNanosec.at(hitPeakTimeIndex + 2),
405  _WaveformADCNanosec.at(hitPeakTimeIndex + 3)
406  };
407 
408  TGraph *gfit = new TGraph(7, x, y);
409 
410  // Make sure fitter is Minuit2 -- works better than Minut here
411  TVirtualFitter::SetDefaultFitter("Minuit2");
412  TFitResultPtr frp = gfit->Fit("pol3", "SFNQ");
413  T p0 = frp->Parameter(0);
414  T p1 = frp->Parameter(1);
415  T p2 = frp->Parameter(2);
416  T p3 = frp->Parameter(3);
417 
418  T disc = p2*p2 - 3.*p1*p3;
419  if (disc < 0) { // Discriminant negative => no extrema
420  hitPeakValue = _WaveformADCNanosec.at(hitPeakTimeIndex);
421  return hitPeakTimeIndex;
422  }
423 
424  // Positive root corresponds to the minimum for a third-order polynomial
425  T hitPeakTime = (-p2 + sqrt(disc))/(3.*p3);
426  hitPeakValue = p0 + p1*hitPeakTime + p2*hitPeakTime*hitPeakTime + p3*hitPeakTime*hitPeakTime*hitPeakTime;
427  hitPeakTime += (T)hitPeakTimeIndex; // translate back to actual time
428  return hitPeakTime;
429 }
430 
431 template<class T>
432 T beamlinereco::CFDHitFinder<T>::IntegrateWaveformInADC(size_t hitStartTimeIndex, size_t hitEndTimeIndex) {
433 
434  T integration = 0;
435 
436  if (!_CFDParamSet[kIntergratedWindowFixed]) {
437  hitStartTimeIndex = hitStartTimeIndex < 0 ? 0 : hitStartTimeIndex;
438  hitEndTimeIndex = hitEndTimeIndex > (_CFDParamSet[kNSamplingPoints] - 1) ? _CFDParamSet[kNSamplingPoints] - 1 : hitEndTimeIndex;
439  for (size_t i = hitStartTimeIndex; i <= hitEndTimeIndex; i++) {
440  integration = integration + (_WaveformADCNanosec.at(i) - _PedestalInADC);
441  }
442  } else {
443  for (size_t i = _CFDParamSet[kIntergratedWindowLowerLimitIndex]; i <= _CFDParamSet[kIntergratedWindowUpperLimitIndex]; i++) {
444  integration = integration + (_WaveformADCNanosec.at(i) - _PedestalInADC);
445  }
446  }
447 
448  return integration;
449 }
450 
451 template<class T>
452 bool beamlinereco::CFDHitFinder<T>::BackwardFindingOfHitStart(size_t hitPeakTimeIndex, T hitPeakValue, T& hitStartTimeIndex, T& hitRiseTimeInIndex) {
453  bool isThisHitContainedFromTheRisingEdge = true;
454 
455  if (_DiscriminationType == "CFD") {
456  _DiscriminationThresholdInADC = _PedestalInADC - _CFDParamSet[kDiscriminationThreshold] * (_PedestalInADC - hitPeakValue);
457  } else if (_DiscriminationType == "LE") {
458  _DiscriminationThresholdInADC = _PedestalInADC - _CFDParamSet[kDiscriminationThreshold] * _NoiseSigmaInADC;
459  } else {
461  << "Unknown DiscriminationType: " << _DiscriminationType << "." << std::endl;
462  }
463 
464  size_t tmp_index = hitPeakTimeIndex;
465  while (_WaveformADCNanosec.at(tmp_index) < _DiscriminationThresholdInADC) {
466  if (tmp_index == 2 && _WaveformADCNanosec.at(tmp_index - 1) < _DiscriminationThresholdInADC) {
467  isThisHitContainedFromTheRisingEdge = false;
468  break;
469  }
470  tmp_index = tmp_index - 1;
471  }
472 
473  if ((_InterpolationType == "LinearNearest") || (tmp_index < 2)) {
474  // Draw a line between two points
475  T V1 = _WaveformADCNanosec.at(tmp_index);
476  T V2 = _WaveformADCNanosec.at(tmp_index - 1);
477  size_t t1 = tmp_index;
478  size_t t2 = tmp_index - 1;
479  T slope = (V2 - V1) / ((T)t2 - (T)t1);
480  hitStartTimeIndex = (_DiscriminationThresholdInADC - V1) / slope + (T)t1;
481  }
482  if ((_InterpolationType == "LinearFit") && (tmp_index >= 2)) {
483  // Fit the surrounding 4
484  T x[4] = {
485  (T)(tmp_index - 2),
486  (T)(tmp_index - 1),
487  (T)(tmp_index),
488  (T)(tmp_index + 1)
489  };
490  T y[4] = {
491  _WaveformADCNanosec.at(tmp_index - 2),
492  _WaveformADCNanosec.at(tmp_index - 1),
493  _WaveformADCNanosec.at(tmp_index),
494  _WaveformADCNanosec.at(tmp_index + 1)
495  };
496  TGraph *gfit = new TGraph(4, x, y);
497  TFitResultPtr frp = gfit->Fit("pol1", "SQNF");
498  double p0 = frp->Parameter(0);
499  double p1 = frp->Parameter(1);
500  T V = p0 + p1 * (T)tmp_index;
501  T t = (T)tmp_index;
502  hitStartTimeIndex = (T)t + (_DiscriminationThresholdInADC - V) / p1;
503  }
504 
505  if (!isThisHitContainedFromTheRisingEdge) {
506  hitRiseTimeInIndex = (T) hitPeakTimeIndex - (T) tmp_index;
507  } else {
508  tmp_index = hitPeakTimeIndex;
509  while (_WaveformADCNanosec.at(tmp_index) <= _PedestalInADC) {
510  if (tmp_index == 2 && _WaveformADCNanosec.at(tmp_index - 1) < _PedestalInADC) {
511  isThisHitContainedFromTheRisingEdge = false;
512  break;
513  }
514  tmp_index = tmp_index - 1;
515  }
516  hitRiseTimeInIndex = (T) hitPeakTimeIndex - (T) tmp_index;
517  }
518 
519  return isThisHitContainedFromTheRisingEdge;
520 }
521 
522 template<class T>
523 bool beamlinereco::CFDHitFinder<T>::ForwardFindingOfHitFallTime(size_t hitPeakTimeIndex, T& hitFallTimeInIndex) {
524  bool isThisHitContainedFromTheFallingEdge = true;
525  size_t tmp_index = hitPeakTimeIndex;
526  while (_WaveformADCNanosec.at(tmp_index) < _PedestalInADC) {
527  if (tmp_index == _CFDParamSet[kNSamplingPoints] - 2 && _WaveformADCNanosec.at(tmp_index + 1) < _PedestalInADC) {
528  isThisHitContainedFromTheFallingEdge = false;
529  break;
530  }
531  tmp_index = tmp_index + 1;
532  }
533 
534  hitFallTimeInIndex = (T)tmp_index - (T)hitPeakTimeIndex;
535  return isThisHitContainedFromTheFallingEdge;
536 }
537 
538 template<class T>
539 T beamlinereco::CFDHitFinder<T>::FindHitWidth(size_t hitPeakTimeIndex, T hitPeakValueADC) {
540  // Calculate FWHM
541 
542  // Find half-max
543  T halfMax = _PedestalInADC - (_PedestalInADC - hitPeakValueADC)/2.;
544 
545  T halfMaxLeft = 9993; // left point of half-max in ticks
546  T halfMaxRight = 9993; // right point of half-max in ticks
547 
548  // Scan backwards
549  for (size_t i=hitPeakTimeIndex; i>0; i--) {
550  if (_WaveformADCNanosec.at(i) == halfMax) { // exactly equal
551  halfMaxLeft = (T)i;
552  break;
553  } else if (_WaveformADCNanosec.at(i) > halfMax) {
554  // Simple two-point linear interpolation
555  T y1 = _WaveformADCNanosec.at(i);
556  T y2 = _WaveformADCNanosec.at(i + 1);
557  size_t t1 = i;
558  size_t t2 = i + 1;
559  T slope = (y2 - y1) / ((T)t2 - (T)t1);
560  halfMaxLeft = (halfMax - y1) / slope + (T)t1;
561  break;
562  }
563  }
564 
565  // Scan forwards
566  for (size_t i=hitPeakTimeIndex; i<_CFDParamSet[kNSamplingPoints]-1; i++) {
567  if (_WaveformADCNanosec.at(i) == halfMax) { // exactly equal
568  halfMaxRight = (T)i;
569  break;
570  } else if (_WaveformADCNanosec.at(i) > halfMax) {
571  // Simple two-point linear interpolation
572  T y1 = _WaveformADCNanosec.at(i - 1);
573  T y2 = _WaveformADCNanosec.at(i);
574  size_t t1 = i - 1;
575  size_t t2 = i;
576  T slope = (y2 - y1) / ((T)t2 - (T)t1);
577  halfMaxRight = (halfMax - y1) / slope + (T)t1;
578  break;
579  }
580  }
581 
582  T widthInTicks = halfMaxRight - halfMaxLeft;
583  return widthInTicks;
584 }
585 
586 template<class T>
588  if (_WaveformADCNanosec.size() == 0) return;
589 
590  FindPedestal(); // Get pedestal and noise band sigma
591  FindRawHitLogic(); // From pedestal and noise band, get rawHitThreshold and RawHitLogic vector
592  FindCFDHits(); // From RawHitLogic vector, get hitPeakValue and hitPeakTimeIndex. From hitPeakValue and pedestal, get cfdHitThresholdValue. Go back from hitPeakTimeIndex to the point when
593 }
594 
595 
596 
597 
598 
599 
600 
601 
602 
603 
604 
605 
606 
607 
608 
609 
610 
611 
612 
613 
614 
615 
616 
617 
618 
619 
620 
621 
622 
623 
624 // GOLAY_SAVITZKY ALGORITHM
625 
626 //! default convergence
627 static const double TINY_FLOAT = 1.0e-300;
628 
629 //! comfortable array of doubles
630 using float_vect = std::vector<double>;
631 //! comfortable array of ints;
632 using int_vect = std::vector<int>;
633 
634 /*! matrix class.
635  *
636  * This is a matrix class derived from a vector of float_vects. Note that
637  * the matrix elements indexed [row][column] with indices starting at 0 (c
638  * style). Also note that because of its design looping through rows should
639  * be faster than looping through columns.
640  *
641  * \brief two dimensional floating point array
642  */
643 class float_mat : public std::vector<float_vect> {
644 private:
645  //! disable the default constructor
646  explicit float_mat() {};
647  //! disable assignment operator until it is implemented.
648  float_mat &operator =(const float_mat &) { return *this; };
649 public:
650  //! constructor with sizes
651  float_mat(const size_t rows, const size_t cols, const double def=0.0);
652  //! copy constructor for matrix
653  float_mat(const float_mat &m);
654  //! copy constructor for vector
655  float_mat(const float_vect &v);
656 
657  //! use default destructor
658  // ~float_mat() {};
659 
660  //! get size
661  size_t nr_rows(void) const { return size(); };
662  //! get size
663  size_t nr_cols(void) const { return front().size(); };
664 };
665 
666 
667 
668 // constructor with sizes
669 float_mat::float_mat(const size_t rows,const size_t cols,const double defval)
670  : std::vector<float_vect>(rows) {
671 
672  for (unsigned int i = 0; i < rows; ++i) {
673  (*this)[i].resize(cols, defval);
674  }
675  if ((rows < 1) || (cols < 1)) {
676  char buffer[1024];
677 
678  sprintf(buffer, "cannot build matrix with %d rows and %d columns\n",
679  (int)rows, (int)cols);
680  //sgs_error(buffer);
681  }
682 }
683 
684 // copy constructor for matrix
686 
687  float_mat::iterator inew = begin();
688  float_mat::const_iterator iold = m.begin();
689  for (/* empty */; iold < m.end(); ++inew, ++iold) {
690  const size_t oldsz = iold->size();
691  inew->resize(oldsz);
692  const float_vect oldvec(*iold);
693  *inew = oldvec;
694  }
695 }
696 
697 // copy constructor for vector
699  : std::vector<float_vect>(1) {
700 
701  const size_t oldsz = v.size();
702  front().resize(oldsz);
703  front() = v;
704 }
705 
706 //////////////////////
707 // Helper functions //
708 //////////////////////
709 
710 //! permute() orders the rows of A to match the integers in the index array.
712 {
713  int_vect i(idx.size());
714  unsigned int j,k;
715 
716  for (j = 0; j < A.nr_rows(); ++j) {
717  i[j] = j;
718  }
719 
720  // loop over permuted indices
721  for (j = 0; j < A.nr_rows(); ++j) {
722  if (i[j] != idx[j]) {
723 
724  // search only the remaining indices
725  for (k = j+1; k < A.nr_rows(); ++k) {
726  if (i[k] ==idx[j]) {
727  std::swap(A[j],A[k]); // swap the rows and
728  i[k] = i[j]; // the elements of
729  i[j] = idx[j]; // the ordered index.
730  break; // next j
731  }
732  }
733  }
734  }
735 }
736 
737 /*! \brief Implicit partial pivoting.
738  *
739  * The function looks for pivot element only in rows below the current
740  * element, A[idx[row]][column], then swaps that row with the current one in
741  * the index map. The algorithm is for implicit pivoting (i.e., the pivot is
742  * chosen as if the max coefficient in each row is set to 1) based on the
743  * scaling information in the vector scale. The map of swapped indices is
744  * recorded in swp. The return value is +1 or -1 depending on whether the
745  * number of row swaps was even or odd respectively. */
746 static int partial_pivot(float_mat &A, const size_t row, const size_t col,
747  float_vect &scale, int_vect &idx, double tol)
748 {
749  if (tol <= 0.0)
750  tol = TINY_FLOAT;
751 
752  int swapNum = 1;
753 
754  // default pivot is the current position, [row,col]
755  unsigned int pivot = row;
756  double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]];
757 
758  // loop over possible pivots below current
759  unsigned int j;
760  for (j = row + 1; j < A.nr_rows(); ++j) {
761 
762  const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]];
763 
764  // if this elem is larger, then it becomes the pivot
765  if (tmp > piv_elem) {
766  pivot = j;
767  piv_elem = tmp;
768  }
769  }
770 
771 #if 0
772  if(piv_elem < tol) {
773  //sgs_error("partial_pivot(): Zero pivot encountered.\n")
774 #endif
775 
776  if(pivot > row) { // bring the pivot to the diagonal
777  j = idx[row]; // reorder swap array
778  idx[row] = idx[pivot];
779  idx[pivot] = j;
780  swapNum = -swapNum; // keeping track of odd or even swap
781  }
782  return swapNum;
783 }
784 
785 /*! \brief Perform backward substitution.
786  *
787  * Solves the system of equations A*b=a, ASSUMING that A is upper
788  * triangular. If diag==1, then the diagonal elements are additionally
789  * assumed to be 1. Note that the lower triangular elements are never
790  * checked, so this function is valid to use after a LU-decomposition in
791  * place. A is not modified, and the solution, b, is returned in a. */
792 static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false)
793 {
794  int r,c;
795  unsigned int k;
796 
797  for (r = (A.nr_rows() - 1); r >= 0; --r) {
798  for (c = (A.nr_cols() - 1); c > r; --c) {
799  for (k = 0; k < A.nr_cols(); ++k) {
800  a[r][k] -= A[r][c] * a[c][k];
801  }
802  }
803  if(!diag) {
804  for (k = 0; k < A.nr_cols(); ++k) {
805  a[r][k] /= A[r][r];
806  }
807  }
808  }
809 }
810 
811 /*! \brief Perform forward substitution.
812  *
813  * Solves the system of equations A*b=a, ASSUMING that A is lower
814  * triangular. If diag==1, then the diagonal elements are additionally
815  * assumed to be 1. Note that the upper triangular elements are never
816  * checked, so this function is valid to use after a LU-decomposition in
817  * place. A is not modified, and the solution, b, is returned in a. */
818 static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true)
819 {
820  unsigned int r, k, c;
821  for (r = 0;r < A.nr_rows(); ++r) {
822  for(c = 0; c < r; ++c) {
823  for (k = 0; k < A.nr_cols(); ++k) {
824  a[r][k] -= A[r][c] * a[c][k];
825  }
826  }
827  if(!diag) {
828  for (k = 0; k < A.nr_cols(); ++k) {
829  a[r][k] /= A[r][r];
830  }
831  }
832  }
833 }
834 
835 /*! \brief Performs LU factorization in place.
836  *
837  * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of
838  * swapped indeces is recorded in idx. The return value is +1 or -1
839  * depending on whether the number of row swaps was even or odd
840  * respectively. idx must be preinitialized to a valid set of indices
841  * (e.g., {1,2, ... ,A.nr_rows()}). */
842 static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
843 {
844  if ( tol <= 0.0)
845  tol = TINY_FLOAT;
846 
847  if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) {
848  //sgs_error("lu_factorize(): cannot handle empty "
849  // "or nonsquare matrices.\n");
850 
851  return 0;
852  }
853 
854  float_vect scale(A.nr_rows()); // implicit pivot scaling
855  unsigned int i;
856  int j;
857  for (i = 0; i < A.nr_rows(); ++i) {
858  double maxval = 0.0;
859  for (j = 0; j < (int)A.nr_cols(); ++j) {
860  if (fabs(A[i][j]) > maxval)
861  maxval = fabs(A[i][j]);
862  }
863  if (maxval == 0.0) {
864  //sgs_error("lu_factorize(): zero pivot found.\n");
865  return 0;
866  }
867  scale[i] = 1.0 / maxval;
868  }
869 
870  int swapNum = 1;
871  unsigned int c,r;
872  for (c = 0; c < A.nr_cols() ; ++c) { // loop over columns
873  swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal
874  for(r = 0; r < A.nr_rows(); ++r) { // loop over rows
875  int lim = (r < c) ? r : c;
876  for (j = 0; j < lim; ++j) {
877  A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c];
878  }
879  if (r > c)
880  A[idx[r]][c] /= A[idx[c]][c];
881  }
882  }
883  permute(A,idx);
884  return swapNum;
885 }
886 
887 /*! \brief Solve a system of linear equations.
888  * Solves the inhomogeneous matrix problem with lu-decomposition. Note that
889  * inversion may be accomplished by setting a to the identity_matrix. */
890 static float_mat lin_solve(const float_mat &A, const float_mat &a,
891  double tol=TINY_FLOAT)
892 {
893  float_mat B(A);
894  float_mat b(a);
895  int_vect idx(B.nr_rows());
896  unsigned int j;
897 
898  for (j = 0; j < B.nr_rows(); ++j) {
899  idx[j] = j; // init row swap label array
900  }
901  lu_factorize(B,idx,tol); // get the lu-decomp.
902  permute(b,idx); // sort the inhomogeneity to match the lu-decomp
903  lu_forwsubst(B,b); // solve the forward problem
904  lu_backsubst(B,b); // solve the backward problem
905  return b;
906 }
907 
908 ///////////////////////
909 // related functions //
910 ///////////////////////
911 
912 //! Returns the inverse of a matrix using LU-decomposition.
913 static float_mat invert(const float_mat &A)
914 {
915  const int n = A.size();
916  float_mat E(n, n, 0.0);
917  float_mat B(A);
918  int i;
919 
920  for (i = 0; i < n; ++i) {
921  E[i][i] = 1.0;
922  }
923 
924  return lin_solve(B, E);
925 }
926 
927 //! returns the transposed matrix.
929 {
930  float_mat res(a.nr_cols(), a.nr_rows());
931  unsigned int i,j;
932 
933  for (i = 0; i < a.nr_rows(); ++i) {
934  for (j = 0; j < a.nr_cols(); ++j) {
935  res[j][i] = a[i][j];
936  }
937  }
938  return res;
939 }
940 
941 //! matrix multiplication.
943 {
944  float_mat res(a.nr_rows(), b.nr_cols());
945  if (a.nr_cols() != b.nr_rows()) {
946  //sgs_error("incompatible matrices in multiplication\n");
947  return res;
948  }
949 
950  unsigned int i,j,k;
951 
952  for (i = 0; i < a.nr_rows(); ++i) {
953  for (j = 0; j < b.nr_cols(); ++j) {
954  double sum(0.0);
955  for (k = 0; k < a.nr_cols(); ++k) {
956  sum += a[i][k] * b[k][j];
957  }
958  res[i][j] = sum;
959  }
960  }
961  return res;
962 }
963 
964 
965 //! calculate savitzky golay coefficients.
966 static float_vect sg_coeff(const float_vect &b, const size_t deg)
967 {
968  const size_t rows(b.size());
969  const size_t cols(deg + 1);
970  float_mat A(rows, cols);
971  float_vect res(rows);
972 
973  // generate input matrix for least squares fit
974  unsigned int i,j;
975  for (i = 0; i < rows; ++i) {
976  for (j = 0; j < cols; ++j) {
977  A[i][j] = pow(double(i), double(j));
978  }
979  }
980 
981  float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));
982 
983  for (i = 0; i < b.size(); ++i) {
984  res[i] = c[0][0];
985  for (j = 1; j <= deg; ++j) {
986  res[i] += c[j][0] * pow(double(i), double(j));
987  }
988  }
989  return res;
990 }
991 
992 /*! \brief savitzky golay smoothing.
993  *
994  * This method means fitting a polynome of degree 'deg' to a sliding window
995  * of width 2w+1 throughout the data. The needed coefficients are
996  * generated dynamically by doing a least squares fit on a "symmetric" unit
997  * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome
998  * yields the sg-coefficients. at the border non symmectric vectors b are
999  * used. */
1000 float_vect beamlinereco::sg_smooth(const float_vect &v, const int width, const int deg)
1001 {
1002  float_vect res(v.size(), 0.0);
1003  if ((width < 1) || (deg < 0) || ((int)v.size() < (2 * width + 2))) {
1004  //sgs_error("sgsmooth: parameter error.\n");
1005  return res;
1006  }
1007 
1008  const int window = 2 * width + 1;
1009  const int endidx = v.size() - 1;
1010 
1011  // do a regular sliding window average
1012  int i,j;
1013  if (deg == 0) {
1014  // handle border cases first because we need different coefficients
1015 #if defined(_OPENMP)
1016 #pragma omp parallel for private(i,j) schedule(static)
1017 #endif
1018  for (i = 0; i < (int)width; ++i) {
1019  const double scale = 1.0/double(i+1);
1020  const float_vect c1(width, scale);
1021  for (j = 0; j <= i; ++j) {
1022  res[i] += c1[j] * v[j];
1023  res[endidx - i] += c1[j] * v[endidx - j];
1024  }
1025  }
1026 
1027  // now loop over rest of data. reusing the "symmetric" coefficients.
1028  const double scale = 1.0/double(window);
1029  const float_vect c2(window, scale);
1030 #if defined(_OPENMP)
1031 #pragma omp parallel for private(i,j) schedule(static)
1032 #endif
1033  for (i = 0; i <= (int) (v.size() - window); ++i) {
1034  for (j = 0; j < (int)window; ++j) {
1035  res[i + width] += c2[j] * v[i + j];
1036  }
1037  }
1038  return res;
1039  }
1040 
1041  // handle border cases first because we need different coefficients
1042 #if defined(_OPENMP)
1043 #pragma omp parallel for private(i,j) schedule(static)
1044 #endif
1045  for (i = 0; i < (int) width; ++i) {
1046  float_vect b1(window, 0.0);
1047  b1[i] = 1.0;
1048 
1049  const float_vect c1(sg_coeff(b1, deg));
1050  for (j = 0; j < (int) window; ++j) {
1051  res[i] += c1[j] * v[j];
1052  res[endidx - i] += c1[j] * v[endidx - j];
1053  }
1054  }
1055 
1056  // now loop over rest of data. reusing the "symmetric" coefficients.
1057  float_vect b2(window, 0.0);
1058  b2[width] = 1.0;
1059  const float_vect c2(sg_coeff(b2, deg));
1060 
1061 #if defined(_OPENMP)
1062 #pragma omp parallel for private(i,j) schedule(static)
1063 #endif
1064  for (i = 0; i <= (int) (v.size() - window); ++i) {
1065  for (j = 0; j < (int) window; ++j) {
1066  res[i + width] += c2[j] * v[i + j];
1067  }
1068  }
1069  return res;
1070 }
1071 
1072 /*! least squares fit a polynome of degree 'deg' to data in 'b'.
1073  * then calculate the first derivative and return it. */
1074 static float_vect lsqr_fprime(const float_vect &b, const int deg)
1075 {
1076  const int rows(b.size());
1077  const int cols(deg + 1);
1078  float_mat A(rows, cols);
1079  float_vect res(rows);
1080 
1081  // generate input matrix for least squares fit
1082  int i,j;
1083  for (i = 0; i < (int) rows; ++i) {
1084  for (j = 0; j < (int) cols; ++j) {
1085  A[i][j] = pow(double(i), double(j));
1086  }
1087  }
1088 
1089  float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));
1090 
1091  for (i = 0; i < (int) b.size(); ++i) {
1092  res[i] = c[1][0];
1093  for (j = 1; j < (int) deg; ++j) {
1094  res[i] += c[j + 1][0] * double(j+1)
1095  * pow(double(i), double(j));
1096  }
1097  }
1098  return res;
1099 }
1100 
1101 /*! \brief savitzky golay smoothed numerical derivative.
1102  *
1103  * This method means fitting a polynome of degree 'deg' to a sliding window
1104  * of width 2w+1 throughout the data.
1105  *
1106  * In contrast to the sg_smooth function we do a brute force attempt by
1107  * always fitting the data to a polynome of degree 'deg' and using the
1108  * result. */
1109 float_vect beamlinereco::sg_derivative(const float_vect &v, const int width,
1110  const int deg, const double h)
1111 {
1112  float_vect res(v.size(), 0.0);
1113  if ((width < 1) || (deg < 1) || ((int) v.size() < (2 * width + 2))) {
1114  //sgs_error("sgsderiv: parameter error.\n");
1115  return res;
1116  }
1117 
1118  const int window = 2 * width + 1;
1119 
1120  // handle border cases first because we do not repeat the fit
1121  // lower part
1122  float_vect b(window, 0.0);
1123  int i,j;
1124 
1125  for (i = 0; i < (int) window; ++i) {
1126  b[i] = v[i] / h;
1127  }
1128  const float_vect c(lsqr_fprime(b, deg));
1129  for (j = 0; j <= (int) width; ++j) {
1130  res[j] = c[j];
1131  }
1132  // upper part. direction of fit is reversed
1133  for (i = 0; i < (int) window; ++i) {
1134  b[i] = v[v.size() - 1 - i] / h;
1135  }
1136  const float_vect d(lsqr_fprime(b, deg));
1137  for (i = 0; i <= (int) width; ++i) {
1138  res[v.size() - 1 - i] = -d[i];
1139  }
1140 
1141  // now loop over rest of data. wasting a lot of least squares calcs
1142  // since we only use the middle value.
1143 #if defined(_OPENMP)
1144 #pragma omp parallel for private(i,j) schedule(static)
1145 #endif
1146  for (i = 1; i < (int) (v.size() - window); ++i) {
1147  for (j = 0; j < (int) window; ++j) {
1148  b[j] = v[i + j] / h;
1149  }
1150  res[i + width] = lsqr_fprime(b, deg)[width];
1151  }
1152  return res;
1153 }
1154 
1155 void beamlinereco::SGSmoothing::Smooth(size_t nsamples, double *in, double *out, const int w, const int deg) {
1156  std::vector<double> in_vec;
1157  for (size_t idx = 0; idx < nsamples; idx++) {
1158  in_vec.push_back(*(in + idx));
1159  }
1160 
1161  std::vector<double> out_vec = sg_smooth(in_vec, w, deg);
1162 
1163  for (size_t idx = 0; idx < nsamples; idx++) {
1164  *(out + idx) = out_vec.at(idx);
1165  }
1166 }
1167 
1168 void beamlinereco::SGSmoothing::Derivative(size_t nsamples, double *in, double *out, const int w, const int deg, const double h) {
1169  std::vector<double> in_vec;
1170  for (size_t idx = 0; idx < nsamples; idx++) {
1171  in_vec.push_back(*(in + idx));
1172  }
1173 
1174  std::vector<double> out_vec = sg_derivative(in_vec, w, deg, h);
1175 
1176  for (size_t idx = 0; idx < nsamples; idx++) {
1177  *(out + idx) = out_vec.at(idx);
1178  }
1179 }
void SetFilterType(const std::string type)
static const double TINY_FLOAT
default convergence
std::map< T, hit_t< T > > _HitCollection
static float_vect lsqr_fprime(const float_vect &b, const int deg)
const std::string GetInterpolationType() const
static int partial_pivot(float_mat &A, const size_t row, const size_t col, float_vect &scale, int_vect &idx, double tol)
Implicit partial pivoting.
static constexpr Double_t deg
Definition: Munits.h:165
const std::map< T, hit_t< T > > & GetHitCollection() const
static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false)
Perform backward substitution.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< T > _WaveformADCNanosec
virtual bool BackwardFindingOfHitStart(size_t hitPeakTimeIndex, T hitPeakValue, T &hitStartTimeIndex, T &hitRiseTimeInIndex)
void permute(float_mat &A, int_vect &idx)
permute() orders the rows of A to match the integers in the index array.
Float_t y1[n_points_granero]
Definition: compare.C:5
static void Smooth(size_t nsamples, double *in, double *out, const int w, const int deg)
static float_mat lin_solve(const float_mat &A, const float_mat &a, double tol=TINY_FLOAT)
Solve a system of linear equations. Solves the inhomogeneous matrix problem with lu-decomposition. Note that inversion may be accomplished by setting a to the identity_matrix.
T sqrt(T number)
Definition: d0nt_math.hpp:156
two dimensional floating point array
constexpr T pow(T x)
Definition: pow.h:75
std::map< CFDParams, T > _CFDParamSet
const T & GetDiscriminationThreshold() const
const std::string GetFilterType() const
Float_t tmp
Definition: plot.C:36
std::vector< double > float_vect
comfortable array of doubles
static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
Performs LU factorization in place.
::xsd::cxx::tree::buffer< char > buffer
Definition: Database.h:179
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
void SetParams(const std::map< CFDParams, double > &paramSet)
Double_t K
static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true)
Perform forward substitution.
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
c2
Definition: demo5.py:33
#define P(a, b, c, d, e, x)
void SetInterpolationType(const std::string type)
virtual T FindPeakValue(size_t hitPeakTimeIndex, T &hitPeakValue)
void SetChannel(unsigned int channelNo)
std::vector< double > sg_derivative(const std::vector< double > &v, const int w, const int deg, const double h=1.0)
float_mat operator*(const float_mat &a, const float_mat &b)
matrix multiplication.
Double_t scale
Definition: plot.C:25
double GetFilteredValue(double measurement)
size_t nr_rows(void) const
use default destructor
const XML_Char int const XML_Char * value
Definition: expat.h:331
const int cols[3]
const std::string GetDiscriminationType() const
Float_t E
Definition: plot.C:20
Int_t col[ntarg]
Definition: Style.C:29
const double a
void Reset()
virtual T FindHitWidth(size_t hitPeakTimeIndex, T hitPeakValueADC)
Float_t d
Definition: plot.C:236
virtual bool ForwardFindingOfHitFallTime(size_t hitPeakTimeIndex, T &hitFallTimeInIndex)
const double j
Definition: BetheBloch.cxx:29
virtual T IntegrateWaveformInADC(size_t hitStartTimeIndex, size_t hitEndTimeIndex)
double t2
void SetWaveform(std::vector< uint16_t > &waveform, unsigned int channelNo, uint32_t timestamp)
static float_vect sg_coeff(const float_vect &b, const size_t deg)
calculate savitzky golay coefficients.
std::vector< double > sg_smooth(const std::vector< double > &v, const int w, const int deg)
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double A
Definition: Units.h:82
ifstream in
Definition: comparison.C:7
size_t nr_cols(void) const
get size
void SetParam(CFDParams param, T value)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
std::vector< T > _NonfilterWaveformADCNanosec
static void Derivative(size_t nsamples, double *in, double *out, const int w, const int deg, const double h=1.0)
const T & GetNoiseSigma() const
const hit & b
Definition: hits.cxx:21
std::vector< int > int_vect
comfortable array of ints;
TRandom3 r(0)
unsigned int DigitizerChannel
c1
Definition: demo5.py:24
void SetDiscriminationType(const std::string type)
const T & GetPedestal() const
double T
Definition: Xdiff_gwt.C:5
Double_t sum
Definition: plot.C:31
static float_mat transpose(const float_mat &a)
returns the transposed matrix.
void SetTimestamp(uint32_t timestamp)
Float_t w
Definition: plot.C:20
float_mat()
disable the default constructor
static float_mat invert(const float_mat &A)
Returns the inverse of a matrix using LU-decomposition.
std::vector< bool > _RawHitLogicNanosec