LEHitFinderAlg.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 
15 
16 namespace beamlinereco {
17  template <class T>
18  struct hit_t {
19  unsigned int DigitizerChannel;
20  uint32_t Timestamp;
29  bool IsContained;
30  };
31 
32  enum LEParams {
33  kADCNBits,
35  kADCOffset,
43  kGSFilter,
49  };
50 
51  std::vector<double> sg_smooth(const std::vector<double> &v, const int w, const int deg);
52  std::vector<double> sg_derivative(const std::vector<double> &v, const int w, const int deg, const double h = 1.0);
53 
54  class SGSmoothing {
55  public:
56  static void Smooth(size_t nsamples, double* in, double* out, const int w, const int deg);
57  static void Derivative(size_t nsamples, double* in, double* out, const int w, const int deg, const double h = 1.0);
58  };
59 
60  template <class T>
61  class LEHitFinder {
62  public:
64  _HitCollection.clear();
65  _RawHitLogicNanosec.clear();
66  _PedestalInADC = 0;
67  _NoiseSigmaInADC = 0;
68  _LEThresholdInADC = 0;
69  _WaveformADCNanosec.clear();
70  _ChannelNumber = 993;
71  _Timestamp = 0;
72  _NonfilterWaveformADCNanosec.clear();
73  _LEParamSet.clear();
74  }
75  virtual ~LEHitFinder() = default;
76 
77 
78  void SetWaveform(std::vector<uint16_t>& waveform, unsigned int channelNo, uint32_t timestamp);
79  void SetParams(const std::map<LEParams, double>& paramSet);
80 
81  inline const std::map<T, hit_t<T> >& GetHitCollection() const { return _HitCollection; }
82  inline const T& GetPedestal() const { return _PedestalInADC; }
83  inline const T& GetNoiseSigma() const { return _NoiseSigmaInADC; }
84  inline const T& GetLEThreshold() const { return _LEThresholdInADC; }
85 
86  virtual void Go();
87 
88  private:
89 
90  void SetChannel(unsigned int channelNo);
91  void SetTimestamp(uint32_t timestamp);
92  void SetParam(LEParams param, T value);
93 
94  private:
95  unsigned int _ChannelNumber;
96  uint32_t _Timestamp;
97  std::vector<T> _WaveformADCNanosec;
99  std::vector<bool> _RawHitLogicNanosec;
103 
104  std::map<LEParams, T> _LEParamSet;
105  std::map<T, hit_t<T> > _HitCollection;
106 
107  private:
108  virtual void FindPedestal();
109  virtual void FindRawHitLogic();
110  virtual void FindLEHits();
111  virtual bool BackwardFindingOfHitStart(size_t hitPeakTimeIndex, T hitPeakValue, T& hitStartTimeIndex, T& hitRiseTimeInIndex);
112  virtual bool ForwardFindingOfHitFallTime(size_t hitPeakTimeIndex, T& hitFallTimeInIndex);
113  virtual T IntegrateWaveformInADC(size_t hitStartTimeIndex, size_t hitEndTimeIndex);
114 
115  virtual inline void Reset() {
116  _HitCollection.clear();
117  _RawHitLogicNanosec.clear();
118  _PedestalInADC = 0;
119  _NoiseSigmaInADC = 0;
120  _LEThresholdInADC = 0;
121  _WaveformADCNanosec.clear();
122  _ChannelNumber = 993;
123  _Timestamp = 0;
124  _NonfilterWaveformADCNanosec.clear();
125  }
126  };
127 }
128 
129 template struct beamlinereco::hit_t<double>;
130 template class beamlinereco::LEHitFinder<double>;
131 
132 template<class T>
133 void beamlinereco::LEHitFinder<T>::SetWaveform(std::vector<uint16_t>& waveform, unsigned int channelNo, uint32_t timestamp) {
134  Reset();
135 
136  SetChannel(channelNo);
137  SetTimestamp(timestamp);
138 
139  double polaritySignFactor = _LEParamSet[kIsWaveformNegativePolarity] ? 1. : -1.;
140 
141  T* nonfilterwaveform = (T*) malloc(_LEParamSet[kNSamplingPoints] * sizeof(T));
142  T* filterwaveform = (T*) malloc(_LEParamSet[kNSamplingPoints] * sizeof(T));
143 
144  for (auto itr = waveform.begin(); itr != waveform.end(); itr++) {
145  _NonfilterWaveformADCNanosec.push_back((T)(*itr) * polaritySignFactor);
146  nonfilterwaveform[itr - waveform.begin()] = (T)(*itr) * polaritySignFactor;
147  }
148 
149  if (_LEParamSet[kGSFilter]) {
151  nonfilterwaveform,
152  filterwaveform,
153  _LEParamSet[kGSFilterWindow],
154  _LEParamSet[kGSFilterDegree]);
155  } else {
156  for (auto itr = waveform.begin(); itr != waveform.end(); itr++) {
157  filterwaveform[itr - waveform.begin()] = (T)(*itr) * polaritySignFactor;
158  }
159  }
160 
161  for (size_t i = 0; i < _LEParamSet[kNSamplingPoints]; i++) {
162  _WaveformADCNanosec.push_back(*(filterwaveform + i));
163  }
164 
165  free(nonfilterwaveform);
166  free(filterwaveform);
167 }
168 
169 template<class T>
171  _LEParamSet[param] = value;
172 }
173 
174 template<class T>
175 void beamlinereco::LEHitFinder<T>::SetParams(const std::map<LEParams, double>& paramSet) {
176  for (auto & itr : paramSet) {
177  SetParam(itr.first, itr.second);
178  }
179 }
180 
181 template<class T>
182 void beamlinereco::LEHitFinder<T>::SetChannel(unsigned int channelNo) {
183  _ChannelNumber = channelNo;
184 }
185 
186 template<class T>
188  _Timestamp = ts;
189 }
190 
191 template<class T>
193 /* 2nd way
194  TH1D* pedHist = new TH1D("", "", 10000, 0, 5000);
195  for (size_t i = 700; i < 1024; i++) {
196  pedHist->Fill(_NonfilterWaveformADCNanosec.at(i));
197  }
198 
199  _PedestalInADC = (T) pedHist->GetMean();
200  _NoiseSigmaInADC = (T) pedHist->GetStdDev();
201 
202  delete (pedHist);
203 */
204 
205 /* 3rd way*/
206  T NoiseBand = 1E9;
207  T OldPedestal = 1E9;
208  T Pedestal = 0;
209 
210  while (TMath::Abs(Pedestal - OldPedestal) >= 2) {
211  T average = 0;
212  T stdDev = 0;
213  size_t countSamples = 0;
214  for (typename std::vector<T>::iterator itr = _NonfilterWaveformADCNanosec.begin(); itr != _NonfilterWaveformADCNanosec.end(); itr++) {
215  if ((*itr <= Pedestal + NoiseBand) && (*itr >= Pedestal - NoiseBand)) {
216  average += *itr;
217  countSamples++;
218  }
219  }
220  average = average / (T)countSamples;
221  for (typename std::vector<T>::iterator itr = _NonfilterWaveformADCNanosec.begin(); itr != _NonfilterWaveformADCNanosec.end(); itr++) {
222  if ((*itr <= Pedestal + NoiseBand) && (*itr >= Pedestal - NoiseBand)) {
223  stdDev += TMath::Power((*itr) - average, 2);
224  }
225  }
226 
227  stdDev = TMath::Sqrt(stdDev / (T)countSamples);
228  OldPedestal = Pedestal;
229  Pedestal = average;
230  NoiseBand = stdDev;
231  }
232 
233  _PedestalInADC = Pedestal;
234  _NoiseSigmaInADC = NoiseBand;
235 /**/
236 }
237 
238 template<class T>
240  T rawHitThreshold = (_PedestalInADC) - ((_NoiseSigmaInADC) * _LEParamSet[kRawHitFinderThresholdInNoiseSigma]);
241  //std::cout << "RawHitThreshold: " << rawHitThreshold << std::endl;
242  for (size_t i = 0; i < _WaveformADCNanosec.size(); i++) {
243  _RawHitLogicNanosec.push_back(_WaveformADCNanosec.at(i) < rawHitThreshold);
244  }
245 
246  for (size_t i = 1; i < _RawHitLogicNanosec.size() - 1; i++) {
247  if ((not _RawHitLogicNanosec.at(i - 1)) and _RawHitLogicNanosec.at(i)) {
248  size_t countDuration = 0;
249  for (size_t j = i; j < _RawHitLogicNanosec.size() - 1; j++) {
250  if (_RawHitLogicNanosec.at(j)) countDuration++;
251  if (not(_RawHitLogicNanosec.at(j + 1))) break;
252  }
253 
254  if (countDuration < _LEParamSet[kShortRawHitIgnoringDurationInTicks]) {
255  for (size_t j = i; j < i + countDuration; j++) {
256  _RawHitLogicNanosec.at(j) = false;
257  }
258  }
259  }
260  }
261 
262  for (size_t i = 1; i < _RawHitLogicNanosec.size() - 1; i++) {
263  if (_RawHitLogicNanosec.at(i - 1) && not(_RawHitLogicNanosec.at(i))) {
264  size_t countDuration = 0;
265  for (size_t j = i; j < _RawHitLogicNanosec.size() - 1; j++) {
266  if (not(_RawHitLogicNanosec.at(j))) countDuration++;
267  if (_RawHitLogicNanosec.at(j + 1)) break;
268  }
269 
270  if (countDuration < _LEParamSet[kConsecutiveHitSeperationDurationInTicks]) {
271  for (size_t j = i; j < i + countDuration; j++) {
272  _RawHitLogicNanosec.at(j) = true;
273  }
274  }
275  }
276  }
277 
278 }
279 
280 template<class T>
282  std::vector<T> hitPeakValueVector; hitPeakValueVector.clear();
283  std::vector<size_t> hitPeakTimeIndexVector; hitPeakTimeIndexVector.clear();
284  std::vector<T> hitStartTimeIndexVector; hitStartTimeIndexVector.clear();
285  std::vector<T> hitRiseTimeInIndexVector; hitRiseTimeInIndexVector.clear();
286  std::vector<T> hitFallTimeInIndexVector; hitFallTimeInIndexVector.clear();
287  std::vector<bool> isHitContainedVector; isHitContainedVector.clear();
288 
289  bool logicFallingPart = true;
290  for (size_t i = 1; i < _RawHitLogicNanosec.size() - 1; i++) {
291  T hitPeakValue = 1E10;
292  size_t hitPeakTimeIndex = 0;
293  if (_RawHitLogicNanosec.at(i) && !(_RawHitLogicNanosec.at(i - 1)) && logicFallingPart) {
294  logicFallingPart = false;
295  for (size_t j = i; j < _RawHitLogicNanosec.size() - 1; j++) {
296  hitPeakValue = hitPeakValue > _WaveformADCNanosec.at(j) ? (T) _WaveformADCNanosec.at(j) : hitPeakValue;
297  hitPeakTimeIndex = hitPeakValue == _WaveformADCNanosec.at(j) ? j : hitPeakTimeIndex;
298  if (_RawHitLogicNanosec.at(j) && !(_RawHitLogicNanosec.at(j + 1)) && !logicFallingPart) {
299  logicFallingPart = true;
300  T hitStartTimeIndex = (T)hitPeakTimeIndex;
301  T hitRiseTimeInIndex = (T)hitPeakTimeIndex;
302  T hitFallTimeInIndex = (T)hitPeakTimeIndex;
303  bool containedFromRisingEdge = BackwardFindingOfHitStart(hitPeakTimeIndex, hitPeakValue, hitStartTimeIndex, hitRiseTimeInIndex);
304  bool containedFromFallingEdge = ForwardFindingOfHitFallTime(hitPeakTimeIndex, hitFallTimeInIndex);
305 
306  hitPeakValueVector.push_back(hitPeakValue);
307  hitPeakTimeIndexVector.push_back(hitPeakTimeIndex);
308  hitStartTimeIndexVector.push_back(hitStartTimeIndex);
309  hitRiseTimeInIndexVector.push_back(hitRiseTimeInIndex);
310  hitFallTimeInIndexVector.push_back(hitFallTimeInIndex);
311  isHitContainedVector.push_back(containedFromRisingEdge && containedFromFallingEdge);
312  break;
313  }
314  }
315  }
316  }
317 
318  for (size_t i = 0; i < hitStartTimeIndexVector.size(); i++) {
319  hit_t<double> aHit;
320  aHit.Timestamp = _Timestamp;
321  aHit.DigitizerChannel = _ChannelNumber;
322  aHit.TStartInNanoSec = hitStartTimeIndexVector.at(i) * _LEParamSet[kTimeSamplingInterval];
323  aHit.TPeakInNanoSec = hitPeakTimeIndexVector.at(i) * _LEParamSet[kTimeSamplingInterval];
324  aHit.AmplitudeInMiliVolt = hitPeakValueVector.at(i) * (_LEParamSet[kADCDynamicRange] / (TMath::Power(2, _LEParamSet[kADCNBits]) - 1));
325  aHit.AmplitudeInADC = hitPeakValueVector.at(i);
326  aHit.RiseTimeInNanoSec = hitRiseTimeInIndexVector.at(i) * _LEParamSet[kTimeSamplingInterval];
327  aHit.FallTimeInNanoSec = hitFallTimeInIndexVector.at(i) * _LEParamSet[kTimeSamplingInterval];
328  aHit.IsContained = isHitContainedVector.at(i);
329  aHit.IntegratedChargeInADCTimeTicks = IntegrateWaveformInADC(hitPeakTimeIndexVector.at(i) - hitRiseTimeInIndexVector.at(i),
330  hitPeakTimeIndexVector.at(i) + hitFallTimeInIndexVector.at(i));
331  aHit.IntegratedChargeInADCNanoSec = aHit.IntegratedChargeInADCTimeTicks * _LEParamSet[kTimeSamplingInterval];
332 
333  _HitCollection[hitStartTimeIndexVector.at(i)] = aHit;
334  }
335 }
336 
337 template<class T>
338 T beamlinereco::LEHitFinder<T>::IntegrateWaveformInADC(size_t hitStartTimeIndex, size_t hitEndTimeIndex) {
339 
340  T integration = 0;
341 
342  if (!_LEParamSet[kIntergratedWindowFixed]) {
343  hitStartTimeIndex = hitStartTimeIndex < 0 ? 0 : hitStartTimeIndex;
344  hitEndTimeIndex = hitEndTimeIndex > (_LEParamSet[kNSamplingPoints] - 1) ? _LEParamSet[kNSamplingPoints] - 1 : hitEndTimeIndex;
345  for (size_t i = hitStartTimeIndex; i <= hitEndTimeIndex; i++) {
346  integration = integration + (_WaveformADCNanosec.at(i) - _PedestalInADC);
347  }
348  } else {
349  for (size_t i = _LEParamSet[kIntergratedWindowLowerLimitIndex]; i <= _LEParamSet[kIntergratedWindowUpperLimitIndex]; i++) {
350  integration = integration + (_WaveformADCNanosec.at(i) - _PedestalInADC);
351  }
352  }
353 
354  return integration;
355 }
356 
357 template<class T>
358 bool beamlinereco::LEHitFinder<T>::BackwardFindingOfHitStart(size_t hitPeakTimeIndex, T hitPeakValue, T& hitStartTimeIndex, T& hitRiseTimeInIndex) {
359  bool isThisHitContainedFromTheRisingEdge = true;
360  _LEThresholdInADC = _PedestalInADC - (_LEParamSet[kLEThresholdInNoiseBands] * (_NoiseSigmaInADC));
361 
362  size_t tmp_index = hitPeakTimeIndex;
363  while (_WaveformADCNanosec.at(tmp_index) < _LEThresholdInADC) {
364  if (tmp_index == 2 && _WaveformADCNanosec.at(tmp_index - 1) < _LEThresholdInADC) {
365  isThisHitContainedFromTheRisingEdge = false;
366  break;
367  }
368  tmp_index = tmp_index - 1;
369  }
370 
371  T V1 = _WaveformADCNanosec.at(tmp_index);
372  T V2 = _WaveformADCNanosec.at(tmp_index - 1);
373  size_t t1 = tmp_index;
374  size_t t2 = tmp_index - 1;
375  T slope = (V2 - V1) / ((T)t2 - (T)t1);
376  hitStartTimeIndex = (_LEThresholdInADC - V1) / slope + (T)t1;
377 
378  if (!isThisHitContainedFromTheRisingEdge) {
379  hitRiseTimeInIndex = (T) hitPeakTimeIndex - (T) tmp_index;
380  } else {
381  tmp_index = hitPeakTimeIndex;
382  while (_WaveformADCNanosec.at(tmp_index) <= _PedestalInADC) {
383  if (tmp_index == 2 && _WaveformADCNanosec.at(tmp_index - 1) < _PedestalInADC) {
384  isThisHitContainedFromTheRisingEdge = false;
385  break;
386  }
387  tmp_index = tmp_index - 1;
388  }
389  hitRiseTimeInIndex = (T) hitPeakTimeIndex - (T) tmp_index;
390  }
391 
392  return isThisHitContainedFromTheRisingEdge;
393 }
394 
395 template<class T>
396 bool beamlinereco::LEHitFinder<T>::ForwardFindingOfHitFallTime(size_t hitPeakTimeIndex, T& hitFallTimeInIndex) {
397  bool isThisHitContainedFromTheFallingEdge = true;
398  size_t tmp_index = hitPeakTimeIndex;
399  while (_WaveformADCNanosec.at(tmp_index) < _PedestalInADC) {
400  if (tmp_index == _LEParamSet[kNSamplingPoints] - 2 && _WaveformADCNanosec.at(tmp_index + 1) < _PedestalInADC) {
401  isThisHitContainedFromTheFallingEdge = false;
402  break;
403  }
404  tmp_index = tmp_index + 1;
405  }
406 
407  hitFallTimeInIndex = (T)tmp_index - (T)hitPeakTimeIndex;
408  return isThisHitContainedFromTheFallingEdge;
409 }
410 
411 template<class T>
413  if (_WaveformADCNanosec.size() == 0) return;
414 
415  FindPedestal(); // Get pedestal and noise band sigma
416  FindRawHitLogic(); // From pedestal and noise band, get rawHitThreshold and RawHitLogic vector
417  FindLEHits(); // From RawHitLogic vector, get hitPeakValue and hitPeakTimeIndex. From noise sigma and pedestal, get leHitThresholdValue. Go back from hitPeakTimeIndex to the point when
418 }
419 
420 
421 
422 
423 
424 
425 
426 
427 
428 
429 
430 
431 
432 
433 
434 
435 
436 
437 
438 
439 
440 
441 
442 
443 
444 
445 
446 
447 
448 
449 // GOLAY_SAVITZKY ALGORITHM
450 
451 //! default convergence
452 static const double TINY_FLOAT = 1.0e-300;
453 
454 //! comfortable array of doubles
455 using float_vect = std::vector<double>;
456 //! comfortable array of ints;
457 using int_vect = std::vector<int>;
458 
459 /*! matrix class.
460  *
461  * This is a matrix class derived from a vector of float_vects. Note that
462  * the matrix elements indexed [row][column] with indices starting at 0 (c
463  * style). Also note that because of its design looping through rows should
464  * be faster than looping through columns.
465  *
466  * \brief two dimensional floating point array
467  */
468 class float_mat : public std::vector<float_vect> {
469 private:
470  //! disable the default constructor
471  explicit float_mat() {};
472  //! disable assignment operator until it is implemented.
473  float_mat &operator =(const float_mat &) { return *this; };
474 public:
475  //! constructor with sizes
476  float_mat(const size_t rows, const size_t cols, const double def=0.0);
477  //! copy constructor for matrix
478  float_mat(const float_mat &m);
479  //! copy constructor for vector
480  float_mat(const float_vect &v);
481 
482  //! use default destructor
483  // ~float_mat() {};
484 
485  //! get size
486  size_t nr_rows(void) const { return size(); };
487  //! get size
488  size_t nr_cols(void) const { return front().size(); };
489 };
490 
491 
492 
493 // constructor with sizes
494 float_mat::float_mat(const size_t rows,const size_t cols,const double defval)
495  : std::vector<float_vect>(rows) {
496 
497  for (unsigned int i = 0; i < rows; ++i) {
498  (*this)[i].resize(cols, defval);
499  }
500  if ((rows < 1) || (cols < 1)) {
501  char buffer[1024];
502 
503  sprintf(buffer, "cannot build matrix with %d rows and %d columns\n",
504  (int)rows, (int)cols);
505  //sgs_error(buffer);
506  }
507 }
508 
509 // copy constructor for matrix
510 float_mat::float_mat(const float_mat &m) : std::vector<float_vect>(m.size()) {
511 
512  float_mat::iterator inew = begin();
513  float_mat::const_iterator iold = m.begin();
514  for (/* empty */; iold < m.end(); ++inew, ++iold) {
515  const size_t oldsz = iold->size();
516  inew->resize(oldsz);
517  const float_vect oldvec(*iold);
518  *inew = oldvec;
519  }
520 }
521 
522 // copy constructor for vector
524  : std::vector<float_vect>(1) {
525 
526  const size_t oldsz = v.size();
527  front().resize(oldsz);
528  front() = v;
529 }
530 
531 //////////////////////
532 // Helper functions //
533 //////////////////////
534 
535 //! permute() orders the rows of A to match the integers in the index array.
537 {
538  int_vect i(idx.size());
539  unsigned int j,k;
540 
541  for (j = 0; j < A.nr_rows(); ++j) {
542  i[j] = j;
543  }
544 
545  // loop over permuted indices
546  for (j = 0; j < A.nr_rows(); ++j) {
547  if (i[j] != idx[j]) {
548 
549  // search only the remaining indices
550  for (k = j+1; k < A.nr_rows(); ++k) {
551  if (i[k] ==idx[j]) {
552  std::swap(A[j],A[k]); // swap the rows and
553  i[k] = i[j]; // the elements of
554  i[j] = idx[j]; // the ordered index.
555  break; // next j
556  }
557  }
558  }
559  }
560 }
561 
562 /*! \brief Implicit partial pivoting.
563  *
564  * The function looks for pivot element only in rows below the current
565  * element, A[idx[row]][column], then swaps that row with the current one in
566  * the index map. The algorithm is for implicit pivoting (i.e., the pivot is
567  * chosen as if the max coefficient in each row is set to 1) based on the
568  * scaling information in the vector scale. The map of swapped indices is
569  * recorded in swp. The return value is +1 or -1 depending on whether the
570  * number of row swaps was even or odd respectively. */
571 static int partial_pivot(float_mat &A, const size_t row, const size_t col,
572  float_vect &scale, int_vect &idx, double tol)
573 {
574  if (tol <= 0.0)
575  tol = TINY_FLOAT;
576 
577  int swapNum = 1;
578 
579  // default pivot is the current position, [row,col]
580  unsigned int pivot = row;
581  double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]];
582 
583  // loop over possible pivots below current
584  unsigned int j;
585  for (j = row + 1; j < A.nr_rows(); ++j) {
586 
587  const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]];
588 
589  // if this elem is larger, then it becomes the pivot
590  if (tmp > piv_elem) {
591  pivot = j;
592  piv_elem = tmp;
593  }
594  }
595 
596 #if 0
597  if(piv_elem < tol) {
598  //sgs_error("partial_pivot(): Zero pivot encountered.\n")
599 #endif
600 
601  if(pivot > row) { // bring the pivot to the diagonal
602  j = idx[row]; // reorder swap array
603  idx[row] = idx[pivot];
604  idx[pivot] = j;
605  swapNum = -swapNum; // keeping track of odd or even swap
606  }
607  return swapNum;
608 }
609 
610 /*! \brief Perform backward substitution.
611  *
612  * Solves the system of equations A*b=a, ASSUMING that A is upper
613  * triangular. If diag==1, then the diagonal elements are additionally
614  * assumed to be 1. Note that the lower triangular elements are never
615  * checked, so this function is valid to use after a LU-decomposition in
616  * place. A is not modified, and the solution, b, is returned in a. */
617 static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false)
618 {
619  int r,c;
620  unsigned int k;
621 
622  for (r = (A.nr_rows() - 1); r >= 0; --r) {
623  for (c = (A.nr_cols() - 1); c > r; --c) {
624  for (k = 0; k < A.nr_cols(); ++k) {
625  a[r][k] -= A[r][c] * a[c][k];
626  }
627  }
628  if(!diag) {
629  for (k = 0; k < A.nr_cols(); ++k) {
630  a[r][k] /= A[r][r];
631  }
632  }
633  }
634 }
635 
636 /*! \brief Perform forward substitution.
637  *
638  * Solves the system of equations A*b=a, ASSUMING that A is lower
639  * triangular. If diag==1, then the diagonal elements are additionally
640  * assumed to be 1. Note that the upper triangular elements are never
641  * checked, so this function is valid to use after a LU-decomposition in
642  * place. A is not modified, and the solution, b, is returned in a. */
643 static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true)
644 {
645  unsigned int r, k, c;
646  for (r = 0;r < A.nr_rows(); ++r) {
647  for(c = 0; c < r; ++c) {
648  for (k = 0; k < A.nr_cols(); ++k) {
649  a[r][k] -= A[r][c] * a[c][k];
650  }
651  }
652  if(!diag) {
653  for (k = 0; k < A.nr_cols(); ++k) {
654  a[r][k] /= A[r][r];
655  }
656  }
657  }
658 }
659 
660 /*! \brief Performs LU factorization in place.
661  *
662  * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of
663  * swapped indeces is recorded in idx. The return value is +1 or -1
664  * depending on whether the number of row swaps was even or odd
665  * respectively. idx must be preinitialized to a valid set of indices
666  * (e.g., {1,2, ... ,A.nr_rows()}). */
667 static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
668 {
669  if ( tol <= 0.0)
670  tol = TINY_FLOAT;
671 
672  if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) {
673  //sgs_error("lu_factorize(): cannot handle empty "
674  // "or nonsquare matrices.\n");
675 
676  return 0;
677  }
678 
679  float_vect scale(A.nr_rows()); // implicit pivot scaling
680  unsigned int i;
681  int j;
682  for (i = 0; i < A.nr_rows(); ++i) {
683  double maxval = 0.0;
684  for (j = 0; j < (int)A.nr_cols(); ++j) {
685  if (fabs(A[i][j]) > maxval)
686  maxval = fabs(A[i][j]);
687  }
688  if (maxval == 0.0) {
689  //sgs_error("lu_factorize(): zero pivot found.\n");
690  return 0;
691  }
692  scale[i] = 1.0 / maxval;
693  }
694 
695  int swapNum = 1;
696  unsigned int c,r;
697  for (c = 0; c < A.nr_cols() ; ++c) { // loop over columns
698  swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal
699  for(r = 0; r < A.nr_rows(); ++r) { // loop over rows
700  int lim = (r < c) ? r : c;
701  for (j = 0; j < lim; ++j) {
702  A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c];
703  }
704  if (r > c)
705  A[idx[r]][c] /= A[idx[c]][c];
706  }
707  }
708  permute(A,idx);
709  return swapNum;
710 }
711 
712 /*! \brief Solve a system of linear equations.
713  * Solves the inhomogeneous matrix problem with lu-decomposition. Note that
714  * inversion may be accomplished by setting a to the identity_matrix. */
715 static float_mat lin_solve(const float_mat &A, const float_mat &a,
716  double tol=TINY_FLOAT)
717 {
718  float_mat B(A);
719  float_mat b(a);
720  int_vect idx(B.nr_rows());
721  unsigned int j;
722 
723  for (j = 0; j < B.nr_rows(); ++j) {
724  idx[j] = j; // init row swap label array
725  }
726  lu_factorize(B,idx,tol); // get the lu-decomp.
727  permute(b,idx); // sort the inhomogeneity to match the lu-decomp
728  lu_forwsubst(B,b); // solve the forward problem
729  lu_backsubst(B,b); // solve the backward problem
730  return b;
731 }
732 
733 ///////////////////////
734 // related functions //
735 ///////////////////////
736 
737 //! Returns the inverse of a matrix using LU-decomposition.
738 static float_mat invert(const float_mat &A)
739 {
740  const int n = A.size();
741  float_mat E(n, n, 0.0);
742  float_mat B(A);
743  int i;
744 
745  for (i = 0; i < n; ++i) {
746  E[i][i] = 1.0;
747  }
748 
749  return lin_solve(B, E);
750 }
751 
752 //! returns the transposed matrix.
754 {
755  float_mat res(a.nr_cols(), a.nr_rows());
756  unsigned int i,j;
757 
758  for (i = 0; i < a.nr_rows(); ++i) {
759  for (j = 0; j < a.nr_cols(); ++j) {
760  res[j][i] = a[i][j];
761  }
762  }
763  return res;
764 }
765 
766 //! matrix multiplication.
768 {
769  float_mat res(a.nr_rows(), b.nr_cols());
770  if (a.nr_cols() != b.nr_rows()) {
771  //sgs_error("incompatible matrices in multiplication\n");
772  return res;
773  }
774 
775  unsigned int i,j,k;
776 
777  for (i = 0; i < a.nr_rows(); ++i) {
778  for (j = 0; j < b.nr_cols(); ++j) {
779  double sum(0.0);
780  for (k = 0; k < a.nr_cols(); ++k) {
781  sum += a[i][k] * b[k][j];
782  }
783  res[i][j] = sum;
784  }
785  }
786  return res;
787 }
788 
789 
790 //! calculate savitzky golay coefficients.
791 static float_vect sg_coeff(const float_vect &b, const size_t deg)
792 {
793  const size_t rows(b.size());
794  const size_t cols(deg + 1);
795  float_mat A(rows, cols);
796  float_vect res(rows);
797 
798  // generate input matrix for least squares fit
799  unsigned int i,j;
800  for (i = 0; i < rows; ++i) {
801  for (j = 0; j < cols; ++j) {
802  A[i][j] = pow(double(i), double(j));
803  }
804  }
805 
806  float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));
807 
808  for (i = 0; i < b.size(); ++i) {
809  res[i] = c[0][0];
810  for (j = 1; j <= deg; ++j) {
811  res[i] += c[j][0] * pow(double(i), double(j));
812  }
813  }
814  return res;
815 }
816 
817 /*! \brief savitzky golay smoothing.
818  *
819  * This method means fitting a polynome of degree 'deg' to a sliding window
820  * of width 2w+1 throughout the data. The needed coefficients are
821  * generated dynamically by doing a least squares fit on a "symmetric" unit
822  * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome
823  * yields the sg-coefficients. at the border non symmectric vectors b are
824  * used. */
825 float_vect beamlinereco::sg_smooth(const float_vect &v, const int width, const int deg)
826 {
827  float_vect res(v.size(), 0.0);
828  if ((width < 1) || (deg < 0) || ((int)v.size() < (2 * width + 2))) {
829  //sgs_error("sgsmooth: parameter error.\n");
830  return res;
831  }
832 
833  const int window = 2 * width + 1;
834  const int endidx = v.size() - 1;
835 
836  // do a regular sliding window average
837  int i,j;
838  if (deg == 0) {
839  // handle border cases first because we need different coefficients
840 #if defined(_OPENMP)
841 #pragma omp parallel for private(i,j) schedule(static)
842 #endif
843  for (i = 0; i < (int)width; ++i) {
844  const double scale = 1.0/double(i+1);
845  const float_vect c1(width, scale);
846  for (j = 0; j <= i; ++j) {
847  res[i] += c1[j] * v[j];
848  res[endidx - i] += c1[j] * v[endidx - j];
849  }
850  }
851 
852  // now loop over rest of data. reusing the "symmetric" coefficients.
853  const double scale = 1.0/double(window);
854  const float_vect c2(window, scale);
855 #if defined(_OPENMP)
856 #pragma omp parallel for private(i,j) schedule(static)
857 #endif
858  for (i = 0; i <= (int) (v.size() - window); ++i) {
859  for (j = 0; j < (int)window; ++j) {
860  res[i + width] += c2[j] * v[i + j];
861  }
862  }
863  return res;
864  }
865 
866  // handle border cases first because we need different coefficients
867 #if defined(_OPENMP)
868 #pragma omp parallel for private(i,j) schedule(static)
869 #endif
870  for (i = 0; i < (int) width; ++i) {
871  float_vect b1(window, 0.0);
872  b1[i] = 1.0;
873 
874  const float_vect c1(sg_coeff(b1, deg));
875  for (j = 0; j < (int) window; ++j) {
876  res[i] += c1[j] * v[j];
877  res[endidx - i] += c1[j] * v[endidx - j];
878  }
879  }
880 
881  // now loop over rest of data. reusing the "symmetric" coefficients.
882  float_vect b2(window, 0.0);
883  b2[width] = 1.0;
884  const float_vect c2(sg_coeff(b2, deg));
885 
886 #if defined(_OPENMP)
887 #pragma omp parallel for private(i,j) schedule(static)
888 #endif
889  for (i = 0; i <= (int) (v.size() - window); ++i) {
890  for (j = 0; j < (int) window; ++j) {
891  res[i + width] += c2[j] * v[i + j];
892  }
893  }
894  return res;
895 }
896 
897 /*! least squares fit a polynome of degree 'deg' to data in 'b'.
898  * then calculate the first derivative and return it. */
899 static float_vect lsqr_fprime(const float_vect &b, const int deg)
900 {
901  const int rows(b.size());
902  const int cols(deg + 1);
903  float_mat A(rows, cols);
904  float_vect res(rows);
905 
906  // generate input matrix for least squares fit
907  int i,j;
908  for (i = 0; i < (int) rows; ++i) {
909  for (j = 0; j < (int) cols; ++j) {
910  A[i][j] = pow(double(i), double(j));
911  }
912  }
913 
914  float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));
915 
916  for (i = 0; i < (int) b.size(); ++i) {
917  res[i] = c[1][0];
918  for (j = 1; j < (int) deg; ++j) {
919  res[i] += c[j + 1][0] * double(j+1)
920  * pow(double(i), double(j));
921  }
922  }
923  return res;
924 }
925 
926 /*! \brief savitzky golay smoothed numerical derivative.
927  *
928  * This method means fitting a polynome of degree 'deg' to a sliding window
929  * of width 2w+1 throughout the data.
930  *
931  * In contrast to the sg_smooth function we do a brute force attempt by
932  * always fitting the data to a polynome of degree 'deg' and using the
933  * result. */
934 float_vect beamlinereco::sg_derivative(const float_vect &v, const int width,
935  const int deg, const double h)
936 {
937  float_vect res(v.size(), 0.0);
938  if ((width < 1) || (deg < 1) || ((int) v.size() < (2 * width + 2))) {
939  //sgs_error("sgsderiv: parameter error.\n");
940  return res;
941  }
942 
943  const int window = 2 * width + 1;
944 
945  // handle border cases first because we do not repeat the fit
946  // lower part
947  float_vect b(window, 0.0);
948  int i,j;
949 
950  for (i = 0; i < (int) window; ++i) {
951  b[i] = v[i] / h;
952  }
953  const float_vect c(lsqr_fprime(b, deg));
954  for (j = 0; j <= (int) width; ++j) {
955  res[j] = c[j];
956  }
957  // upper part. direction of fit is reversed
958  for (i = 0; i < (int) window; ++i) {
959  b[i] = v[v.size() - 1 - i] / h;
960  }
961  const float_vect d(lsqr_fprime(b, deg));
962  for (i = 0; i <= (int) width; ++i) {
963  res[v.size() - 1 - i] = -d[i];
964  }
965 
966  // now loop over rest of data. wasting a lot of least squares calcs
967  // since we only use the middle value.
968 #if defined(_OPENMP)
969 #pragma omp parallel for private(i,j) schedule(static)
970 #endif
971  for (i = 1; i < (int) (v.size() - window); ++i) {
972  for (j = 0; j < (int) window; ++j) {
973  b[j] = v[i + j] / h;
974  }
975  res[i + width] = lsqr_fprime(b, deg)[width];
976  }
977  return res;
978 }
979 
980 void beamlinereco::SGSmoothing::Smooth(size_t nsamples, double *in, double *out, const int w, const int deg) {
981  std::vector<double> in_vec;
982  for (size_t idx = 0; idx < nsamples; idx++) {
983  in_vec.push_back(*(in + idx));
984  }
985 
986  std::vector<double> out_vec = sg_smooth(in_vec, w, deg);
987 
988  for (size_t idx = 0; idx < nsamples; idx++) {
989  *(out + idx) = out_vec.at(idx);
990  }
991 }
992 
993 void beamlinereco::SGSmoothing::Derivative(size_t nsamples, double *in, double *out, const int w, const int deg, const double h) {
994  std::vector<double> in_vec;
995  for (size_t idx = 0; idx < nsamples; idx++) {
996  in_vec.push_back(*(in + idx));
997  }
998 
999  std::vector<double> out_vec = sg_derivative(in_vec, w, deg, h);
1000 
1001  for (size_t idx = 0; idx < nsamples; idx++) {
1002  *(out + idx) = out_vec.at(idx);
1003  }
1004 }
static float_vect lsqr_fprime(const float_vect &b, const int deg)
const T & GetPedestal() const
static constexpr Double_t deg
Definition: Munits.h:165
virtual bool BackwardFindingOfHitStart(size_t hitPeakTimeIndex, T hitPeakValue, T &hitStartTimeIndex, T &hitRiseTimeInIndex)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
static void Smooth(size_t nsamples, double *in, double *out, const int w, const int deg)
two dimensional floating point array
constexpr T pow(T x)
Definition: pow.h:75
virtual void FindRawHitLogic()
std::map< LEParams, T > _LEParamSet
static float_mat transpose(const float_mat &a)
returns the transposed matrix.
virtual T IntegrateWaveformInADC(size_t hitStartTimeIndex, size_t hitEndTimeIndex)
Float_t tmp
Definition: plot.C:36
std::vector< double > float_vect
comfortable array of doubles
void SetParams(const std::map< LEParams, double > &paramSet)
::xsd::cxx::tree::buffer< char > buffer
Definition: Database.h:179
static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false)
Perform backward substitution.
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
c2
Definition: demo5.py:33
std::vector< double > sg_derivative(const std::vector< double > &v, const int w, const int deg, const double h=1.0)
Double_t scale
Definition: plot.C:25
static float_mat invert(const float_mat &A)
Returns the inverse of a matrix using LU-decomposition.
static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true)
Perform forward substitution.
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]
std::map< T, hit_t< T > > _HitCollection
Float_t E
Definition: plot.C:20
const T & GetLEThreshold() const
Int_t col[ntarg]
Definition: Style.C:29
const double a
void Reset()
void permute(float_mat &A, int_vect &idx)
permute() orders the rows of A to match the integers in the index array.
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.
Float_t d
Definition: plot.C:236
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.
const double j
Definition: BetheBloch.cxx:29
void SetTimestamp(uint32_t timestamp)
double t2
static const double TINY_FLOAT
default convergence
std::vector< double > sg_smooth(const std::vector< double > &v, const int w, const int deg)
static const double A
Definition: Units.h:82
ifstream in
Definition: comparison.C:7
size_t nr_cols(void) const
get size
void SetWaveform(std::vector< uint16_t > &waveform, unsigned int channelNo, uint32_t timestamp)
float_mat operator*(const float_mat &a, const float_mat &b)
matrix multiplication.
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
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
void SetChannel(unsigned int channelNo)
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
const std::map< T, hit_t< T > > & GetHitCollection() const
void SetParam(LEParams param, T value)
std::vector< bool > _RawHitLogicNanosec
static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
Performs LU factorization in place.
static float_vect sg_coeff(const float_vect &b, const size_t deg)
calculate savitzky golay coefficients.
std::vector< T > _WaveformADCNanosec
double T
Definition: Xdiff_gwt.C:5
std::vector< T > _NonfilterWaveformADCNanosec
Double_t sum
Definition: plot.C:31
Float_t w
Definition: plot.C:20
float_mat()
disable the default constructor
virtual bool ForwardFindingOfHitFallTime(size_t hitPeakTimeIndex, T &hitFallTimeInIndex)