Classes | Namespaces | Typedefs | Enumerations | Functions | Variables
LEHitFinderAlg.h File Reference
#include <map>
#include <iostream>
#include <vector>
#include <cstdio>
#include <cstddef>
#include <cmath>
#include <TH1D.h>
#include <TF1.h>
#include <TMath.h>

Go to the source code of this file.

Classes

struct  beamlinereco::hit_t< T >
 
class  beamlinereco::SGSmoothing
 
class  beamlinereco::LEHitFinder< T >
 
class  float_mat
 two dimensional floating point array More...
 

Namespaces

 beamlinereco
 

Typedefs

using float_vect = std::vector< double >
 comfortable array of doubles More...
 
using int_vect = std::vector< int >
 comfortable array of ints; More...
 

Enumerations

enum  beamlinereco::LEParams {
  beamlinereco::kADCNBits, beamlinereco::kADCDynamicRange, beamlinereco::kADCOffset, beamlinereco::kTimeSamplingInterval,
  beamlinereco::kNSamplingPoints, beamlinereco::kIsWaveformNegativePolarity, beamlinereco::kLEThresholdInNoiseBands, beamlinereco::kRawHitFinderThresholdInNoiseSigma,
  beamlinereco::kShortRawHitIgnoringDurationInTicks, beamlinereco::kConsecutiveHitSeperationDurationInTicks, beamlinereco::kGSFilter, beamlinereco::kGSFilterWindow,
  beamlinereco::kGSFilterDegree, beamlinereco::kIntergratedWindowFixed, beamlinereco::kIntergratedWindowLowerLimitIndex, beamlinereco::kIntergratedWindowUpperLimitIndex
}
 

Functions

std::vector< double > beamlinereco::sg_smooth (const std::vector< double > &v, const int w, const int deg)
 
std::vector< double > beamlinereco::sg_derivative (const std::vector< double > &v, const int w, const int deg, const double h=1.0)
 
void permute (float_mat &A, int_vect &idx)
 permute() orders the rows of A to match the integers in the index array. More...
 
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. More...
 
static void lu_backsubst (float_mat &A, float_mat &a, bool diag=false)
 Perform backward substitution. More...
 
static void lu_forwsubst (float_mat &A, float_mat &a, bool diag=true)
 Perform forward substitution. More...
 
static int lu_factorize (float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
 Performs LU factorization in place. More...
 
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. More...
 
static float_mat invert (const float_mat &A)
 Returns the inverse of a matrix using LU-decomposition. More...
 
static float_mat transpose (const float_mat &a)
 returns the transposed matrix. More...
 
float_mat operator* (const float_mat &a, const float_mat &b)
 matrix multiplication. More...
 
static float_vect sg_coeff (const float_vect &b, const size_t deg)
 calculate savitzky golay coefficients. More...
 
static float_vect lsqr_fprime (const float_vect &b, const int deg)
 

Variables

static const double TINY_FLOAT = 1.0e-300
 default convergence More...
 

Typedef Documentation

using float_vect = std::vector<double>

comfortable array of doubles

Definition at line 455 of file LEHitFinderAlg.h.

using int_vect = std::vector<int>

comfortable array of ints;

Definition at line 457 of file LEHitFinderAlg.h.

Function Documentation

static float_mat invert ( const float_mat A)
static

Returns the inverse of a matrix using LU-decomposition.

Definition at line 738 of file LEHitFinderAlg.h.

References E, MECModelEnuComparisons::i, lin_solve(), and getGoodRuns4SAM::n.

Referenced by lsqr_fprime(), and sg_coeff().

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 }
two dimensional floating point array
Float_t E
Definition: plot.C:20
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.
static float_mat lin_solve ( const float_mat A,
const float_mat a,
double  tol = TINY_FLOAT 
)
static

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.

Definition at line 715 of file LEHitFinderAlg.h.

References b, compare_h5_caf::idx, calib::j, lu_backsubst(), lu_factorize(), lu_forwsubst(), float_mat::nr_rows(), and permute().

Referenced by invert().

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 }
two dimensional floating point array
static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false)
Perform backward substitution.
static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true)
Perform forward substitution.
void permute(float_mat &A, int_vect &idx)
permute() orders the rows of A to match the integers in the index array.
const double j
Definition: BetheBloch.cxx:29
const hit & b
Definition: hits.cxx:21
std::vector< int > int_vect
comfortable array of ints;
static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
Performs LU factorization in place.
static float_vect lsqr_fprime ( const float_vect b,
const int  deg 
)
static

least squares fit a polynome of degree 'deg' to data in 'b'. then calculate the first derivative and return it.

Definition at line 899 of file LEHitFinderAlg.h.

References genie::units::A, b, plot_validation_datamc::c, cols, d, beamlinereco::SGSmoothing::Derivative(), make_syst_table_plots::h, MECModelEnuComparisons::i, compare_h5_caf::idx, in, makeTrainCVSamples::int, invert(), calib::j, confusionMatrixTree::out, cet::pow(), fillBadChanDBTables::rows, beamlinereco::sg_derivative(), beamlinereco::sg_smooth(), beamlinereco::SGSmoothing::Smooth(), transpose(), w, and allInOneTrainingPlots::window.

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 }
static constexpr Double_t deg
Definition: Munits.h:165
two dimensional floating point array
constexpr T pow(T x)
Definition: pow.h:75
static float_mat transpose(const float_mat &a)
returns the transposed matrix.
std::vector< double > float_vect
comfortable array of doubles
static float_mat invert(const float_mat &A)
Returns the inverse of a matrix using LU-decomposition.
const int cols[3]
const double j
Definition: BetheBloch.cxx:29
static const double A
Definition: Units.h:82
const hit & b
Definition: hits.cxx:21
static void lu_backsubst ( float_mat A,
float_mat a,
bool  diag = false 
)
static

Perform backward substitution.

Solves the system of equations A*b=a, ASSUMING that A is upper triangular. If diag==1, then the diagonal elements are additionally assumed to be 1. Note that the lower triangular elements are never checked, so this function is valid to use after a LU-decomposition in place. A is not modified, and the solution, b, is returned in a.

Definition at line 617 of file LEHitFinderAlg.h.

References plot_validation_datamc::c, float_mat::nr_cols(), float_mat::nr_rows(), and r().

Referenced by lin_solve().

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 }
size_t nr_rows(void) const
use default destructor
size_t nr_cols(void) const
get size
TRandom3 r(0)
static int lu_factorize ( float_mat A,
int_vect idx,
double  tol = TINY_FLOAT 
)
static

Performs LU factorization in place.

This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of swapped indeces is recorded in idx. The return value is +1 or -1 depending on whether the number of row swaps was even or odd respectively. idx must be preinitialized to a valid set of indices (e.g., {1,2, ... ,A.nr_rows()}).

Definition at line 667 of file LEHitFinderAlg.h.

References plot_validation_datamc::c, stan::math::fabs(), MECModelEnuComparisons::i, makeTrainCVSamples::int, calib::j, float_mat::nr_cols(), float_mat::nr_rows(), partial_pivot(), permute(), r(), scale, and TINY_FLOAT.

Referenced by lin_solve().

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 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< double > float_vect
comfortable array of doubles
Double_t scale
Definition: plot.C:25
size_t nr_rows(void) const
use default destructor
void permute(float_mat &A, int_vect &idx)
permute() orders the rows of A to match the integers in the index array.
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
static const double TINY_FLOAT
default convergence
size_t nr_cols(void) const
get size
TRandom3 r(0)
static void lu_forwsubst ( float_mat A,
float_mat a,
bool  diag = true 
)
static

Perform forward substitution.

Solves the system of equations A*b=a, ASSUMING that A is lower triangular. If diag==1, then the diagonal elements are additionally assumed to be 1. Note that the upper triangular elements are never checked, so this function is valid to use after a LU-decomposition in place. A is not modified, and the solution, b, is returned in a.

Definition at line 643 of file LEHitFinderAlg.h.

References plot_validation_datamc::c, float_mat::nr_cols(), float_mat::nr_rows(), and r().

Referenced by lin_solve().

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 }
size_t nr_rows(void) const
use default destructor
size_t nr_cols(void) const
get size
TRandom3 r(0)
float_mat operator* ( const float_mat a,
const float_mat b 
)

matrix multiplication.

Definition at line 767 of file LEHitFinderAlg.h.

References MECModelEnuComparisons::i, calib::j, float_mat::nr_cols(), float_mat::nr_rows(), and sum.

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 }
two dimensional floating point array
size_t nr_rows(void) const
use default destructor
const double j
Definition: BetheBloch.cxx:29
size_t nr_cols(void) const
get size
Double_t sum
Definition: plot.C:31
static int partial_pivot ( float_mat A,
const size_t  row,
const size_t  col,
float_vect scale,
int_vect idx,
double  tol 
)
static

Implicit partial pivoting.

The function looks for pivot element only in rows below the current element, A[idx[row]][column], then swaps that row with the current one in the index map. The algorithm is for implicit pivoting (i.e., the pivot is chosen as if the max coefficient in each row is set to 1) based on the scaling information in the vector scale. The map of swapped indices is recorded in swp. The return value is +1 or -1 depending on whether the number of row swaps was even or odd respectively.

Definition at line 571 of file LEHitFinderAlg.h.

References stan::math::fabs(), calib::j, float_mat::nr_rows(), check_grl::row, TINY_FLOAT, and tmp.

Referenced by lu_factorize().

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 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t tmp
Definition: plot.C:36
Double_t scale
Definition: plot.C:25
size_t nr_rows(void) const
use default destructor
Int_t col[ntarg]
Definition: Style.C:29
const double j
Definition: BetheBloch.cxx:29
static const double TINY_FLOAT
default convergence
void permute ( float_mat A,
int_vect idx 
)

permute() orders the rows of A to match the integers in the index array.

Definition at line 536 of file LEHitFinderAlg.h.

References MECModelEnuComparisons::i, calib::j, float_mat::nr_rows(), and std::swap().

Referenced by lin_solve(), and lu_factorize().

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 }
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
size_t nr_rows(void) const
use default destructor
const double j
Definition: BetheBloch.cxx:29
std::vector< int > int_vect
comfortable array of ints;
static float_vect sg_coeff ( const float_vect b,
const size_t  deg 
)
static

calculate savitzky golay coefficients.

Definition at line 791 of file LEHitFinderAlg.h.

References genie::units::A, plot_validation_datamc::c, demo5::c1, demo5::c2, cols, Munits::deg, MECModelEnuComparisons::i, makeTrainCVSamples::int, invert(), calib::j, cet::pow(), fillBadChanDBTables::rows, beamlinereco::sg_smooth(), transpose(), registry_explorer::v, and allInOneTrainingPlots::window.

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 }
static constexpr Double_t deg
Definition: Munits.h:165
two dimensional floating point array
constexpr T pow(T x)
Definition: pow.h:75
static float_mat transpose(const float_mat &a)
returns the transposed matrix.
std::vector< double > float_vect
comfortable array of doubles
static float_mat invert(const float_mat &A)
Returns the inverse of a matrix using LU-decomposition.
const int cols[3]
const double j
Definition: BetheBloch.cxx:29
static const double A
Definition: Units.h:82
const hit & b
Definition: hits.cxx:21
static float_mat transpose ( const float_mat a)
static

returns the transposed matrix.

Definition at line 753 of file LEHitFinderAlg.h.

References MECModelEnuComparisons::i, calib::j, float_mat::nr_cols(), and float_mat::nr_rows().

Referenced by lsqr_fprime(), and sg_coeff().

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 }
two dimensional floating point array
size_t nr_rows(void) const
use default destructor
const double j
Definition: BetheBloch.cxx:29
size_t nr_cols(void) const
get size

Variable Documentation

const double TINY_FLOAT = 1.0e-300
static

default convergence

Definition at line 452 of file LEHitFinderAlg.h.

Referenced by lu_factorize(), and partial_pivot().