WaveFitter.h

#ifndef WAVEFITTER_H
#define WAVEFITTER_H
 
#include <array>
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
 
namespace Belle2 {
  namespace SVD {
 
    class DefaultWave {
    public:
 
      DefaultWave(double tau = s_default_tau) : m_tau(tau) {}
 
      static double s_default_tau;  
      inline double getValue(double t) const
      {
        if (t < 0)
          return 0.0;
        else {
          double z = t / m_tau;
          return z * std::exp(1.0 - z);
        }
      }
 
      double getTau() const { return m_tau; }
 
    private:
      double m_tau; 
    };
 
 
    class WaveFitter {
    public:
 
      typedef std::function<double(double)> wave_function_type;
      typedef std::vector< std::array<double, 6> > strip_data_type;
 
      WaveFitter(wave_function_type wave, std::array<double, 6> times):
        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];
      }
 
      void addData(const std::array<double, 6>& inpData, double rmsNoise)
      {
        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;
      }
 
      void reset()
      {
        m_data.clear();
        m_noises.clear();
        m_fittedData.clear();
        m_fittedAmplitudes.clear();
        m_fittedAmplitudeErrors.clear();
        m_ndf = 0;
        m_hasFit = false;
      }
 
      double getFittedTime()
      {
        if (!m_hasFit) doFit();
        return m_fittedTime;
      }
 
      double getFittedTimeError()
      {
        if (!m_hasFit) doFit();
        return m_fittedTimeError;
      }
 
      const std::vector<double>& getFittedAmplitudes()
      {
        if (!m_hasFit) doFit();
        return m_fittedAmplitudes;
      }
 
      const std::vector<double>& getFittedAmplitudeErrors()
      {
        if (!m_hasFit) doFit();
        return m_fittedAmplitudeErrors;
      }
 
      double getFitLikelihood()
      {
        if (!m_hasFit) doFit();
        return m_fittedLik;
      }
 
      const strip_data_type& getData() const { return m_data; }
 
      const strip_data_type& getFitData()
      {
        if (!m_hasFit) doFit();
        return m_fittedData;
      }
 
      double getDt() const { return m_dt; }
      const std::array<double, 6>& getTimes() const { return m_times; }
 
      double Chi(double t);
      double negLogLikelihood(double t);
      double pSignal();
      double lrSignal(double a, double b);
 
      double wave(double t) const { return m_wave(t); }
      static void setPrecision(int ndigits);
 
    private:
 
      static int s_minimizer_precision;
 
      void doFit();
      void calculateAmplitudes();
      void calculateFitErrors();
      void calculateFittedData();
      double integral12(double lower, double upper, std::function<double(double)> f);
      double integral20(double lower, double upper, std::function<double(double)> f);
      bool m_hasFit;    
      double m_dt;  
      strip_data_type m_data;   
      std::vector<double> m_noises; 
      std::array<double, 6> m_times; 
      wave_function_type m_wave;  
      double m_fittedTime; 
      double m_fittedTimeError; 
      std::vector<double> m_fittedAmplitudes; 
      std::vector<double> m_fittedAmplitudeErrors;  
      double m_fittedLik; 
      int m_ndf; 
      strip_data_type m_fittedData; 
    };
 
  } // namespace SVD
} // namespace Belle2
#endif // WAVEFITTER_H