Classes | Namespaces | Typedefs | Enumerations | Functions | Variables
CFDHitFinderAlg.h File Reference
#include <map>
#include <iostream>
#include <vector>
#include <cstdio>
#include <cstddef>
#include <cmath>
#include <TH1D.h>
#include <TF1.h>
#include <TMath.h>
#include <TGraph.h>
#include <TFitResult.h>
#include <TFitResultPtr.h>
#include <TVirtualFitter.h>
#include "KalmanFilterAlg.h"

Go to the source code of this file.

Classes

struct  beamlinereco::hit_t< T >
 
class  beamlinereco::SGSmoothing
 
class  beamlinereco::CFDHitFinder< 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::CFDParams {
  beamlinereco::kADCNBits, beamlinereco::kADCDynamicRange, beamlinereco::kADCOffset, beamlinereco::kTimeSamplingInterval,
  beamlinereco::kNSamplingPoints, beamlinereco::kIsWaveformNegativePolarity, beamlinereco::kDiscriminationThreshold, beamlinereco::kRawHitFinderThresholdInNoiseSigma,
  beamlinereco::kRawHitFinderTicksFromEnd, beamlinereco::kShortRawHitIgnoringDurationInTicks, beamlinereco::kConsecutiveHitSeperationDurationInTicks, beamlinereco::kGSFilter,
  beamlinereco::kGSFilterWindow, beamlinereco::kGSFilterDegree, beamlinereco::kKalmanFilterProcessNoiseCovariance, beamlinereco::kKalmanFilterMeasurementNoiseCovariance,
  beamlinereco::kKalmanFilterGain, 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 630 of file CFDHitFinderAlg.h.

using int_vect = std::vector<int>

comfortable array of ints;

Definition at line 632 of file CFDHitFinderAlg.h.

Function Documentation

static float_mat invert ( const float_mat A)
static

Returns the inverse of a matrix using LU-decomposition.

Definition at line 913 of file CFDHitFinderAlg.h.

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

Referenced by lsqr_fprime(), and sg_coeff().

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 }
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.
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 
)
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 890 of file CFDHitFinderAlg.h.

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

Referenced by invert().

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

References genie::units::A, b, plot_validation_datamc::c, cols, d, make_syst_table_plots::h, MECModelEnuComparisons::i, makeTrainCVSamples::int, invert(), calib::j, cet::pow(), fillBadChanDBTables::rows, beamlinereco::sg_derivative(), transpose(), and allInOneTrainingPlots::window.

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 }
static constexpr Double_t deg
Definition: Munits.h:165
two dimensional floating point array
constexpr T pow(T x)
Definition: pow.h:75
std::vector< double > float_vect
comfortable array of doubles
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)
returns the transposed matrix.
static float_mat invert(const float_mat &A)
Returns the inverse of a matrix using LU-decomposition.
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 792 of file CFDHitFinderAlg.h.

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

Referenced by lin_solve().

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 }
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 842 of file CFDHitFinderAlg.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().

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 }
static const double TINY_FLOAT
default convergence
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.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void permute(float_mat &A, int_vect &idx)
permute() orders the rows of A to match the integers in the index array.
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
const double j
Definition: BetheBloch.cxx:29
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 818 of file CFDHitFinderAlg.h.

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

Referenced by lin_solve().

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 }
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 942 of file CFDHitFinderAlg.h.

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

Referenced by art::const_AssnsIter< L, R, D, Dir >::const_AssnsIter().

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 }
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 746 of file CFDHitFinderAlg.h.

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

Referenced by lu_factorize().

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 }
static const double TINY_FLOAT
default convergence
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
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 711 of file CFDHitFinderAlg.h.

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

Referenced by lin_solve(), and lu_factorize().

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 }
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 966 of file CFDHitFinderAlg.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.

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 }
static constexpr Double_t deg
Definition: Munits.h:165
two dimensional floating point array
constexpr T pow(T x)
Definition: pow.h:75
std::vector< double > float_vect
comfortable array of doubles
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)
returns the transposed matrix.
static float_mat invert(const float_mat &A)
Returns the inverse of a matrix using LU-decomposition.
static float_mat transpose ( const float_mat a)
static

returns the transposed matrix.

Definition at line 928 of file CFDHitFinderAlg.h.

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

Referenced by lsqr_fprime(), sg_coeff(), visualisationForPaperMasterPlot::vis_square(), visualisationForPaperMasterPlot::vis_square_kernel(), and visualisationForPaperMasterPlot::vis_square_kernel_zoom().

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 }
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 627 of file CFDHitFinderAlg.h.

Referenced by lu_factorize(), and partial_pivot().