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().

{
    const int n = A.size();
    float_mat E(n, n, 0.0);
    float_mat B(A);
    int i;

    for (i = 0; i < n; ++i) {
        E[i][i] = 1.0;
    }

    return lin_solve(B, E);
}

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().

{
    float_mat B(A);
    float_mat b(a);
    int_vect idx(B.nr_rows());
    unsigned int j;

    for (j = 0; j < B.nr_rows(); ++j) {
        idx[j] = j;  // init row swap label array
    }
    lu_factorize(B,idx,tol); // get the lu-decomp.
    permute(b,idx);          // sort the inhomogeneity to match the lu-decomp
    lu_forwsubst(B,b);       // solve the forward problem
    lu_backsubst(B,b);       // solve the backward problem
    return b;
}

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.

{
    const int rows(b.size());
    const int cols(deg + 1);
    float_mat A(rows, cols);
    float_vect res(rows);

    // generate input matrix for least squares fit
    int i,j;
    for (i = 0; i < (int) rows; ++i) {
        for (j = 0; j < (int) cols; ++j) {
            A[i][j] = pow(double(i), double(j));
        }
    }

    float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));

    for (i = 0; i < (int) b.size(); ++i) {
        res[i] = c[1][0];
        for (j = 1; j < (int) deg; ++j) {
            res[i] += c[j + 1][0] * double(j+1)
                      * pow(double(i), double(j));
        }
    }
    return res;
}

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().

{
    int r,c;
    unsigned int k;

    for (r = (A.nr_rows() - 1); r >= 0; --r) {
        for (c = (A.nr_cols() - 1); c > r; --c) {
            for (k = 0; k < A.nr_cols(); ++k) {
                a[r][k] -= A[r][c] * a[c][k];
            }
        }
        if(!diag) {
            for (k = 0; k < A.nr_cols(); ++k) {
                a[r][k] /= A[r][r];
            }
        }
    }
}

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().

{
    if ( tol <= 0.0)
        tol = TINY_FLOAT;

    if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) {
        //sgs_error("lu_factorize(): cannot handle empty "
        //           "or nonsquare matrices.\n");

        return 0;
    }

    float_vect scale(A.nr_rows());  // implicit pivot scaling
    unsigned int i;
    int j;
    for (i = 0; i < A.nr_rows(); ++i) {
        double maxval = 0.0;
        for (j = 0; j < (int)A.nr_cols(); ++j) {
            if (fabs(A[i][j]) > maxval)
                maxval = fabs(A[i][j]);
        }
        if (maxval == 0.0) {
            //sgs_error("lu_factorize(): zero pivot found.\n");
            return 0;
        }
        scale[i] = 1.0 / maxval;
    }

    int swapNum = 1;
    unsigned int c,r;
    for (c = 0; c < A.nr_cols() ; ++c) {            // loop over columns
        swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal
        for(r = 0; r < A.nr_rows(); ++r) {      //  loop over rows
            int lim = (r < c) ? r : c;
            for (j = 0; j < lim; ++j) {
                A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c];
            }
            if (r > c)
                A[idx[r]][c] /= A[idx[c]][c];
        }
    }
    permute(A,idx);
    return swapNum;
}

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().

{
    unsigned int r, k, c;
    for (r = 0;r < A.nr_rows(); ++r) {
        for(c = 0; c < r; ++c) {
            for (k = 0; k < A.nr_cols(); ++k) {
                a[r][k] -= A[r][c] * a[c][k];
            }
        }
        if(!diag) {
            for (k = 0; k < A.nr_cols(); ++k) {
                a[r][k] /= A[r][r];
            }
        }
    }
}

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.

{
    float_mat res(a.nr_rows(), b.nr_cols());
    if (a.nr_cols() != b.nr_rows()) {
        //sgs_error("incompatible matrices in multiplication\n");
        return res;
    }

    unsigned int i,j,k;

    for (i = 0; i < a.nr_rows(); ++i) {
        for (j = 0; j < b.nr_cols(); ++j) {
            double sum(0.0);
            for (k = 0; k < a.nr_cols(); ++k) {
                sum += a[i][k] * b[k][j];
            }
            res[i][j] = sum;
        }
    }
    return res;
}

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().

{
    if (tol <= 0.0)
        tol = TINY_FLOAT;

    int swapNum = 1;

    // default pivot is the current position, [row,col]
    unsigned int pivot = row;
    double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]];

    // loop over possible pivots below current
    unsigned int j;
    for (j = row + 1; j < A.nr_rows(); ++j) {

        const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]];

        // if this elem is larger, then it becomes the pivot
        if (tmp > piv_elem) {
            pivot = j;
            piv_elem = tmp;
        }
    }

#if 0
    if(piv_elem < tol) {
      //sgs_error("partial_pivot(): Zero pivot encountered.\n")
#endif

    if(pivot > row) {           // bring the pivot to the diagonal
        j = idx[row];           // reorder swap array
        idx[row] = idx[pivot];
        idx[pivot] = j;
        swapNum = -swapNum;     // keeping track of odd or even swap
    }
    return swapNum;
}

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().

{
    int_vect i(idx.size());
    unsigned int j,k;

    for (j = 0; j < A.nr_rows(); ++j) {
        i[j] = j;
    }

    // loop over permuted indices
    for (j = 0; j < A.nr_rows(); ++j) {
        if (i[j] != idx[j]) {

            // search only the remaining indices
            for (k = j+1; k < A.nr_rows(); ++k) {
                if (i[k] ==idx[j]) {
                    std::swap(A[j],A[k]); // swap the rows and
                    i[k] = i[j];     // the elements of
                    i[j] = idx[j];   // the ordered index.
                    break; // next j
                }
            }
        }
    }
}

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.

{
    const size_t rows(b.size());
    const size_t cols(deg + 1);
    float_mat A(rows, cols);
    float_vect res(rows);

    // generate input matrix for least squares fit
    unsigned int i,j;
    for (i = 0; i < rows; ++i) {
        for (j = 0; j < cols; ++j) {
            A[i][j] = pow(double(i), double(j));
        }
    }

    float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));

    for (i = 0; i < b.size(); ++i) {
        res[i] = c[0][0];
        for (j = 1; j <= deg; ++j) {
            res[i] += c[j][0] * pow(double(i), double(j));
        }
    }
    return res;
}

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().

{
    float_mat res(a.nr_cols(), a.nr_rows());
    unsigned int i,j;

    for (i = 0; i < a.nr_rows(); ++i) {
        for (j = 0; j < a.nr_cols(); ++j) {
            res[j][i] = a[i][j];
        }
    }
    return res;
}

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().