Waveform fitter class. More...

#include <WaveFitter.h>

Collaboration diagram for WaveFitter:

[legend]

Public Types
typedef std::function< double(double)>	wave_function_type
	Function type for waveform.

typedef std::vector< std::array< double, 6 > >	strip_data_type
	Function type for cluster data.

Public Member Functions
	WaveFitter (wave_function_type wave, std::array< double, 6 > times)
	Constructor creates an empty fitter with default settings.

void	addData (const std::array< double, 6 > &inpData, double rmsNoise)
	Add strip data to the fitter.

void	reset ()
	Reset fitter data.

double	getFittedTime ()
	Retrieve fitted time shift value after a fit.

double	getFittedTimeError ()
	Retrieve error of time shift value after a fit.

const std::vector< double > &	getFittedAmplitudes ()
	Retrieve strip amplitudes after a fit.

const std::vector< double > &	getFittedAmplitudeErrors ()
	Retrieve errors of strip amplitudes after a fit.

double	getFitLikelihood ()
	Return negative log-likelihood for the fit.

const strip_data_type &	getData () const
	Get original fitted data.

const strip_data_type &	getFitData ()
	Get fitted values for all samples.

double	getDt () const
	Get sampling interval.

const std::array< double, 6 > &	getTimes () const
	Get sample times.

double	Chi (double t)
	Get standardized chi2 for a given time shift.

double	negLogLikelihood (double t)
	Get negative profile log-likelihood for a given time shift.

double	pSignal ()
	Calculate probability of the hit being a signal, with given priors for signal and background.

double	lrSignal (double a, double b)
	Calculate likelihood ratio Lsignal/Lbackground for acceptance (signal) window given as arguments.

double	wave (double t) const
	Get unit waveform value at a given time.

Static Public Member Functions
static void	setPrecision (int ndigits)
	Set precision of the minimzer.

Private Member Functions
void	doFit ()
	Perform fit on the data.

void	calculateAmplitudes ()
	Calculate amplitudes.

void	calculateFitErrors ()
	Calculate fit errors.

void	calculateFittedData ()
	Calculate fitted data.

double	integral12 (double lower, double upper, std::function< double(double)> f)
	High-order Gauss-Legendre quadrature for likelihood integrals.

double	integral20 (double lower, double upper, std::function< double(double)> f)
	High-order Gauss-Legendre quadrature for likelihood integrals.

Private Attributes
bool	m_hasFit
	Are fit results available?

double	m_dt
	Time interval between samples.

strip_data_type	m_data
	Vector of sextets of APV samples.

std::vector< double >	m_noises
	RMS noise for strips.

std::array< double, 6 >	m_times
	Vector of sample times.

wave_function_type	m_wave
	Waveform function.

double	m_fittedTime
	Time from interval (-dt,+dt).

double	m_fittedTimeError
	Error estimate of the fitted time.

std::vector< double >	m_fittedAmplitudes
	Fitted amplitudes of the sample.

std::vector< double >	m_fittedAmplitudeErrors
	Fitted amplitude errors.

double	m_fittedLik
	Chi2 from the current fit.

int	m_ndf
	Degrees of freedom of the fit.

strip_data_type	m_fittedData
	Fitted values from current fit.

Static Private Attributes
static int	s_minimizer_precision = 6
	minimizer precision

Detailed Description

Waveform fitter class.

Definition at line 59 of file WaveFitter.h.

Member Typedef Documentation

◆ strip_data_type

typedef std::vector< std::array<double, 6> > strip_data_type

Function type for cluster data.

Definition at line 65 of file WaveFitter.h.

◆ wave_function_type

typedef std::function<double(double)> wave_function_type

Function type for waveform.

Definition at line 63 of file WaveFitter.h.

Constructor & Destructor Documentation

◆ WaveFitter()

WaveFitter	(	wave_function_type	wave,
		std::array< double, 6 >	times
	)

inline

Constructor creates an empty fitter with default settings.

Parameters

wave	Wave function to use for fitting data
times	array of points in time at which data are taken

Definition at line 71 of file WaveFitter.h.

                                                                    :
        m_hasFit(false), m_times(times), m_wave(wave), m_fittedTime(0), m_fittedTimeError(), m_fittedLik(1.0e10), m_ndf(0)
      {
        m_dt = m_times[1] - m_times[0];
      }

Member Function Documentation

◆ addData()

void addData	(	const std::array< double, 6 > &	inpData,
		double	rmsNoise
	)

inline

Add strip data to the fitter.

Parameters

inpData	array of 6 APV signals.
rmsNoise	RMS noise for the strip.

Definition at line 81 of file WaveFitter.h.

      {
        std::array<double, 6> v;
        std::copy(inpData.begin(), inpData.end(), v.begin());
        m_data.push_back(v);
        m_noises.push_back(rmsNoise);
        m_hasFit = false;
        m_ndf = 6 * m_data.size() - m_data.size() - 2;
      }

◆ calculateAmplitudes()

void calculateAmplitudes ( )

private

Calculate amplitudes.

Definition at line 115 of file WaveFitter.cc.

{
  // This we only need to calculate once
  array<double, 6> waves;
  double sumTemplateByTemplate = 0;
  for (int i = 0; i < 6; ++i) {
    waves[i] = m_wave(m_times[i] - m_fittedTime);
    sumTemplateByTemplate += waves[i] * waves[i];
  }
  // Calculate amplitudes
  for (unsigned int row = 0; row < m_data.size(); ++row) {
    double sumDataByTemplate = 0;
    for (int i = 0; i < 6; i++) {
      sumDataByTemplate += m_data[row][i] * waves[i];
    }
    m_fittedAmplitudes.push_back(sumDataByTemplate / sumTemplateByTemplate);
  }
}

◆ calculateFitErrors()

void calculateFitErrors ( )

private

Calculate fit errors.

This only calculates local (linearization) errors.

Definition at line 134 of file WaveFitter.cc.

{
  // Errors on amplitudes - simple from the linear model
  double sum_w2 = 0;
  for (int i = 0; i < 6; ++i) {
    double wave = m_wave(m_times[i] - m_fittedTime);
    sum_w2 += wave * wave;
  }
  double one_by_sqrtSumw2 = 1.0 / sqrt(sum_w2);
  for (unsigned int row = 0; row < m_data.size(); ++row)
    m_fittedAmplitudeErrors.push_back(one_by_sqrtSumw2 * m_noises[row]);
  // For time: ignore correlation with amplitudes
  // We will need derivatives, calculate them numerically
  double timestep = 0.1 * m_dt;
  double dwave_sq = 0;
  for (int i = 0; i < 6; ++i) {
    double dwave = (wave(m_times[i] - m_fittedTime + 0.5 * timestep) - wave(m_times[i] - m_fittedTime - 0.5 * timestep)) / timestep;
    dwave_sq += dwave * dwave;
  }
  double sn_sq = 0;
  for (unsigned int row = 0; row < m_data.size(); ++row) {
    double sn = m_fittedAmplitudes[row] / m_noises[row];
    sn_sq += sn * sn;
  }
  m_fittedTimeError = 1.0 / sqrt(sn_sq * dwave_sq);
}

◆ calculateFittedData()

void calculateFittedData ( )

private

Calculate fitted data.

Definition at line 161 of file WaveFitter.cc.

{
  m_fittedData.resize(m_data.size());
  for (int i = 0; i < 6; ++i) {
    double wave = m_wave(m_times[i] - m_fittedTime);
    for (unsigned int row = 0; row < m_data.size(); ++row)
      m_fittedData[row][i] = m_fittedAmplitudes[row] * wave;
  }
}

◆ Chi()

double Chi ( double t )

Get standardized chi2 for a given time shift.

Parameters

t	time shift

Returns: chi2/ndf at a given time shift

Definition at line 24 of file WaveFitter.cc.

{
  // This we only need to calculate once
  array<double, 6> waves;
  double sumTemplateByTemplate = 0;
  for (int i = 0; i < 6; ++i) {
    waves[i] = m_wave(m_times[i] - t0);
    sumTemplateByTemplate += waves[i] * waves[i];
  }
  // Calculate amplitudes
  vector<double> amplitudes(m_data.size());
  for (unsigned int row = 0; row < m_data.size(); ++row) {
    double sumDataByTemplate = 0;
    for (int i = 0; i < 6; i++) {
      sumDataByTemplate += m_data[row][i] * waves[i];
    }
    amplitudes[row] = sumDataByTemplate / sumTemplateByTemplate;
  }
  double result = 0;
  for (unsigned int row = 0; row < m_data.size(); ++row)
    for (int i = 0; i < 6; i++) {
      double diff =
        (m_data[row][i] - amplitudes[row] * waves[i]) / m_noises[row];
      result += diff * diff;
    }
  return result / m_ndf;
}

◆ doFit()

void doFit ( )

private

Perform fit on the data.

Definition at line 73 of file WaveFitter.cc.

{
  // Degrees of freedom:
  // # samples - # strips (amplitudes) - 1 (timeshift) - 1 (sum = 0)
  m_ndf = 6 * m_data.size() - m_data.size() - 2;
  if (m_ndf < 1) { // We have no data, so there'll be no fit and an error
    B2ERROR("SVD::WaveFitterter : Fit required with no data provided.");
    m_fittedTime = 0;
    m_fittedTimeError = 0;
    m_fittedAmplitudes.clear();
    m_fittedAmplitudeErrors.clear();
    m_fittedLik = 1.0e10;
    m_fittedData.clear();
    // We don't need to repeat this
    m_hasFit = true;
    return;
  }
  // Normal operation: minimize wrt. t0 in (-dt,0) and (0,dt),
  // take the solution giving a better fit
  auto result_low = boost::math::tools::brent_find_minima(
                      [this](double t)->double { return negLogLikelihood(t); }, -m_dt, 0.0, s_minimizer_precision
                    );
  double t0_low = result_low.first;
  double lik_low = result_low.second;
  auto result_hi = boost::math::tools::brent_find_minima(
                     [this](double t)->double { return negLogLikelihood(t); }, 0.0, m_dt, s_minimizer_precision
                   );
  double t0_hi = result_hi.first;
  double lik_hi = result_hi.second;
  if (lik_low < lik_hi) {
    m_fittedTime = t0_low;
    m_fittedLik = lik_low;
  } else {
    m_fittedTime = t0_hi;
    m_fittedLik = lik_hi;
  }
  calculateAmplitudes();
  calculateFitErrors();
  calculateFittedData();
  m_hasFit = true;
}

◆ getData()

const strip_data_type & getData ( ) const

inline

Get original fitted data.

Returns: vector of arrays of strip data used in the fit.

Definition at line 151 of file WaveFitter.h.

151{ return m_data; }

◆ getDt()

double getDt ( ) const

inline

Get sampling interval.

Definition at line 163 of file WaveFitter.h.

163{ return m_dt; }

◆ getFitData()

const strip_data_type & getFitData ( )

inline

Get fitted values for all samples.

Returns: vector of arrays of fit estimates of fitted data.

Definition at line 156 of file WaveFitter.h.

      {
        if (!m_hasFit) doFit();
        return m_fittedData;
      }

◆ getFitLikelihood()

double getFitLikelihood ( )

inline

Return negative log-likelihood for the fit.

Returns: negative log-likelihood for the fit

Definition at line 142 of file WaveFitter.h.

      {
        if (!m_hasFit) doFit();
        return m_fittedLik;
      }

◆ getFittedAmplitudeErrors()

const std::vector< double > & getFittedAmplitudeErrors ( )

inline

Retrieve errors of strip amplitudes after a fit.

Returns: vector<double> of errors of strip amplitude estimates

Definition at line 133 of file WaveFitter.h.

      {
        if (!m_hasFit) doFit();
        return m_fittedAmplitudeErrors;
      }

◆ getFittedAmplitudes()

const std::vector< double > & getFittedAmplitudes ( )

inline

Retrieve strip amplitudes after a fit.

Returns: vector<double> of strip amplitude estimates

Definition at line 124 of file WaveFitter.h.

      {
        if (!m_hasFit) doFit();
        return m_fittedAmplitudes;
      }

◆ getFittedTime()

double getFittedTime ( )

inline

Retrieve fitted time shift value after a fit.

Returns: time shift estimate

Definition at line 106 of file WaveFitter.h.

      {
        if (!m_hasFit) doFit();
        return m_fittedTime;
      }

◆ getFittedTimeError()

double getFittedTimeError ( )

inline

Retrieve error of time shift value after a fit.

Returns: time shift estimate error

Definition at line 115 of file WaveFitter.h.

      {
        if (!m_hasFit) doFit();
        return m_fittedTimeError;
      }

◆ getTimes()

const std::array< double, 6 > & getTimes ( ) const

inline

Get sample times.

Definition at line 165 of file WaveFitter.h.

165{ return m_times; }

◆ integral12()

double integral12	(	double	lower,
		double	upper,
		std::function< double(double)>	f
	)

private

High-order Gauss-Legendre quadrature for likelihood integrals.

Definition at line 171 of file WaveFitter.cc.

{
  const unsigned int half_order = 6;
  const double knots[half_order] = {
    0.1252334085114689154724414,
    0.3678314989981801937526915,
    0.5873179542866174472967024,
    0.7699026741943046870368938,
    0.9041172563704748566784659,
    0.9815606342467192506905491
  };
  const double weights[half_order] = {
    0.2491470458134027850005624,
    0.2334925365383548087608499,
    0.2031674267230659217490645,
    0.1600783285433462263346525,
    0.1069393259953184309602547,
    0.0471753363865118271946160
  };
  double span = 0.5 * (upper - lower);
  double center = 0.5 * (upper + lower);
  double result = 0;
  for (unsigned int iknot = 0; iknot < half_order; ++iknot) {
    result += weights[iknot] * f(center + span * knots[iknot]);
    result += weights[iknot] * f(center - span * knots[iknot]);
  }
  result *= span;
  return result;
}

◆ integral20()

double integral20	(	double	lower,
		double	upper,
		std::function< double(double)>	f
	)

private

High-order Gauss-Legendre quadrature for likelihood integrals.

Definition at line 201 of file WaveFitter.cc.

{
  const unsigned int half_order = 10;
  const double knots[half_order] = {
    0.0765265211,
    0.2277858511,
    0.3737060887,
    0.510867002,
    0.6360536807,
    0.7463319065,
    0.8391169718,
    0.9122344283,
    0.9639719273,
    0.9931285992
  };
  const double weights[half_order] = {
    0.1527533871,
    0.1491729865,
    0.1420961093,
    0.1316886384,
    0.118194532,
    0.1019301198,
    0.0832767416,
    0.0626720483,
    0.0406014298,
    0.0176140071
  };
  double span = 0.5 * (upper - lower);
  double center = 0.5 * (upper + lower);
  double result = 0;
  for (unsigned int iknot = 0; iknot < half_order; ++iknot) {
    result += weights[iknot] * f(center + span * knots[iknot]);
    result += weights[iknot] * f(center - span * knots[iknot]);
  }
  result *= span;
  return result;
}

◆ lrSignal()

double lrSignal	(	double	a,
		double	b
	)

Calculate likelihood ratio Lsignal/Lbackground for acceptance (signal) window given as arguments.

Parameters

a	lower bound of the acceptance window
b	upper bound of the acceptance window

Returns: L(in acc. w.)/L(out of acc. w.)

Definition at line 257 of file WaveFitter.cc.

{
  const double outer = m_dt;
  if (!m_hasFit) doFit();
  auto integrand = [this](double t)->double { return exp(-negLogLikelihood(t));};
  // Calculate total integral and integral between a and b:
  double int_total = integral20(-outer, 0, integrand);
  int_total += integral20(0, outer, integrand);
  double int_inner = 0;
  if (low * high < 0) {
    int_inner += integral20(low, 0, integrand);
    int_inner += integral20(0, high, integrand);
  } else {
    int_inner += integral20(low, high, integrand);
  }
  double int_outer = int_total - int_inner;
  double LR = int_inner / int_outer;
  //cout << L_background << " <><> " << L_signal << " <> " << prob << endl;
  return LR;
}

◆ negLogLikelihood()

double negLogLikelihood ( double t )

Get negative profile log-likelihood for a given time shift.

Parameters

t	the time shift at which to calculate the likelihood

Returns: log-likelihood integrated over amplitudes at the given t

Definition at line 52 of file WaveFitter.cc.

{
  // wave values
  array<double, 6> w;
  for (int i = 0; i < 6; ++i) w[i] = m_wave(m_times[i] - t);
  double sum_w2 = 0;
  for (int i = 0; i < 6; ++i) sum_w2 += w[i] * w[i];
  double lognorm = 0.5 * m_data.size() * log(sum_w2);
  double sum_ksi = 0.0;
  for (unsigned int row = 0; row < m_data.size(); ++row) {
    double ksi = 0;
    for (int i = 1; i < 6; ++i)
      for (int j = 0; j < i; ++j) {
        double cross = w[i] * m_data[row][j] - w[j] * m_data[row][i];
        ksi += cross * cross;
      }
    sum_ksi += ksi;
  }
  return lognorm + 0.5 * sum_ksi / sum_w2;
}

◆ pSignal()

double pSignal ( )

Calculate probability of the hit being a signal, with given priors for signal and background.

Definition at line 239 of file WaveFitter.cc.

{
  const double outer = m_dt;
  const double inner = 5; // ns
  if (!m_hasFit) doFit();
  auto integrand = [this](double t)->double { return exp(-negLogLikelihood(t));};
 
  double i1 = integral20(-outer, -inner, integrand);
  double i2 = integral20(-inner, 0, integrand);
  double i3 = integral20(0, inner, integrand);
  double i4 = integral20(inner, outer, integrand);
  double L_background = 1.0 / (2 * outer) * (i1 + i2 + i3 + i4);
  double L_signal = 1.0 / (2 * inner) * (i2 + i3);
  double prob = (L_signal / (L_background + L_signal));
  //cout << L_background << " <><> " << L_signal << " <> " << prob << endl;
  return prob;
}

◆ reset()

void reset ( )

inline

Reset fitter data.

Definition at line 92 of file WaveFitter.h.

      {
        m_data.clear();
        m_noises.clear();
        m_fittedData.clear();
        m_fittedAmplitudes.clear();
        m_fittedAmplitudeErrors.clear();
        m_ndf = 0;
        m_hasFit = false;
      }

◆ setPrecision()

void setPrecision ( int ndigits )

static

Set precision of the minimzer.

Definition at line 22 of file WaveFitter.cc.

22{ s_minimizer_precision = ndigits; }

◆ wave()

double wave ( double t ) const

inline

Get unit waveform value at a given time.

Definition at line 188 of file WaveFitter.h.

188{ return m_wave(t); }

Member Data Documentation

◆ m_data

strip_data_type m_data

private

Vector of sextets of APV samples.

Definition at line 214 of file WaveFitter.h.

◆ m_dt

double m_dt

private

Time interval between samples.

Definition at line 213 of file WaveFitter.h.

◆ m_fittedAmplitudeErrors

std::vector<double> m_fittedAmplitudeErrors

private

Fitted amplitude errors.

Definition at line 221 of file WaveFitter.h.

◆ m_fittedAmplitudes

std::vector<double> m_fittedAmplitudes

private

Fitted amplitudes of the sample.

Definition at line 220 of file WaveFitter.h.

◆ m_fittedData

strip_data_type m_fittedData

private

Fitted values from current fit.

Definition at line 224 of file WaveFitter.h.

◆ m_fittedLik

double m_fittedLik

private

Chi2 from the current fit.

Definition at line 222 of file WaveFitter.h.

◆ m_fittedTime

double m_fittedTime

private

Time from interval (-dt,+dt).

Definition at line 218 of file WaveFitter.h.

◆ m_fittedTimeError

double m_fittedTimeError

private

Error estimate of the fitted time.

Definition at line 219 of file WaveFitter.h.

◆ m_hasFit

bool m_hasFit

private

Are fit results available?

Definition at line 212 of file WaveFitter.h.

◆ m_ndf

int m_ndf

private

Degrees of freedom of the fit.

Definition at line 223 of file WaveFitter.h.

◆ m_noises

std::vector<double> m_noises

private

RMS noise for strips.

Definition at line 215 of file WaveFitter.h.

◆ m_times

std::array<double, 6> m_times

private

Vector of sample times.

Definition at line 216 of file WaveFitter.h.

◆ m_wave

wave_function_type m_wave

private

Waveform function.

Definition at line 217 of file WaveFitter.h.

◆ s_minimizer_precision

int s_minimizer_precision = 6

staticprivate

minimizer precision

Definition at line 195 of file WaveFitter.h.

The documentation for this class was generated from the following files:

svd/reconstruction/include/WaveFitter.h
svd/reconstruction/src/WaveFitter.cc

Public Types

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Member Typedef Documentation

◆ strip_data_type

◆ wave_function_type

Constructor & Destructor Documentation

◆ WaveFitter()

Member Function Documentation

◆ addData()

◆ calculateAmplitudes()

◆ calculateFitErrors()

◆ calculateFittedData()

◆ Chi()

◆ doFit()

◆ getData()

◆ getDt()

◆ getFitData()

◆ getFitLikelihood()

◆ getFittedAmplitudeErrors()

◆ getFittedAmplitudes()

◆ getFittedTime()

◆ getFittedTimeError()

◆ getTimes()

◆ integral12()

◆ integral20()

◆ lrSignal()

◆ negLogLikelihood()

◆ pSignal()

◆ reset()

◆ setPrecision()

◆ wave()

Member Data Documentation

◆ m_data

◆ m_dt

◆ m_fittedAmplitudeErrors

◆ m_fittedAmplitudes

◆ m_fittedData

◆ m_fittedLik

◆ m_fittedTime

◆ m_fittedTimeError

◆ m_hasFit

◆ m_ndf

◆ m_noises

◆ m_times

◆ m_wave

◆ s_minimizer_precision