NNWaveFitTool.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
#ifndef _SVD_RECONSTRUCTION_NNWAVEFITTOOL_H
#define _SVD_RECONSTRUCTION_NNWAVEFITTOOL_H
 
#include <vector>
#include <memory>
#include <tuple>
#include <algorithm>
#include <numeric>
#include <functional>
#include <svd/simulation/SVDSimulationTools.h>
 
namespace Belle2 {
  namespace SVD {
 
    typedef std::vector<double> nnFitterBinData;
 
    typedef std::vector<double> nnFitterBins;
 
    class EmpiricalDistributionFunction {
    public:
      EmpiricalDistributionFunction(const nnFitterBinData& p, const nnFitterBins& bins):
        m_bins(bins)
      {
        // Create table of probabiilities for EmpiricalDistributionFunction calculation
        m_pTable.resize(m_bins.size());
        auto ibin = m_bins.begin();
        auto ipt = m_pTable.begin();
        double cp = 0.0;
        *ipt++ = std::make_pair(*ibin++, cp);
        for (auto prob : p) {
          cp += prob;
          *ipt++ = std::make_pair(*ibin++, cp);
        }
      }
 
      double operator()(double t)
      {
        // pTable is sorted by .first
        if (t > m_pTable.back().first) return 1.0;
        if (t < m_pTable[0].first) return 0.0;
        auto it = lower_bound(m_pTable.begin(), m_pTable.end(), std::make_pair(t, 0.0));
        if (it == m_pTable.begin()) return it->second;
        auto it2 = it;
        --it2;
        double result =  it2->second + (it->second - it2->second) * (t - it2->first) / (it->first - it2->first);
        return result;
      }
 
    private:
      const nnFitterBins& m_bins; 
      std::vector<std::pair<double, double> > m_pTable; 
    }; // class EmpiricalDistributionFunction
 
    // ========================================================================
    // NNWaveFitTool - object to handle computations with neural network fit.
    // ------------------------------------------------------------------------
    class NNWaveFitTool {
    public:
      NNWaveFitTool(const nnFitterBins& bins, const nnFitterBinData& binCenters,
                    const WaveGenerator& wave):
        m_bins(bins), m_binCenters(binCenters), m_waveGenerator(wave)
      {
        m_altBinData.resize(m_binCenters.size());
      }
 
      //------- Operations on probability arrays -------//
 
      const nnFitterBins& getBins() const { return m_bins; }
 
      const nnFitterBinData& getBinCenters() const { return m_binCenters; }
 
      void shiftInTime(nnFitterBinData& p, double timeShift);
 
      void multiply(nnFitterBinData& p, const nnFitterBinData& p1)
      {
        std::transform(p.begin(), p.end(),
                       p1.begin(), p.begin(), std::multiplies<double>());
        normalize(p);
      }
 
      std::shared_ptr<nnFitterBinData> pFromInterval(double left, double right);
 
      //------- Waveform parameters ------------------------------------//
 
      std::tuple<double, double> getTimeShift(const nnFitterBinData& p);
 
      std::tuple<double, double, double>
      getAmplitudeChi2(const apvSamples& samples, double timeShift, double tau);
 
      double pLessThan(nnFitterBinData p1, nnFitterBinData p2);
 
    protected:
 
      void normalize(nnFitterBinData& p)
      {
        double pnorm = std::accumulate(p.begin(), p.end(), 0.0);
        // If the norm is too small, return default distribution.
        if (pnorm < 1.0e-10) {
          double uniformP = 1.0 / std::distance(p.begin(), p.end());
          std::fill(p.begin(), p.end(), uniformP);
        } else {
          // Normalize if the norm makes sense.
          std::transform(p.begin(), p.end(), p.begin(),
                         std::bind2nd(std::divides<double>(), pnorm));
        }
      }
 
    private:
      const nnFitterBins& m_bins; 
      const nnFitterBinData& m_binCenters; 
      nnFitterBinData m_altBinData; 
      apvSamples m_altSamples; 
      WaveGenerator m_waveGenerator; 
    }; // class NNWaveFitTool
 
  } // SVD namespace
} // Belle2 namespace
 
#endif