SVDSimulationTools.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2010 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Peter Kvasnicka                                          *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#pragma once
#ifndef _SVD_SIMULATION_SIMULATIONTOOLS_H
#define _SVD_SIMULATION_SIMULATIONTOOLS_H
 
#include <cmath>
#include <array>
#include <algorithm>
#include <functional>
#include <sstream>
#include <framework/gearbox/Unit.h>
 
namespace Belle2 {
  namespace SVD {
 
// ==============================================================================
// APV25 sample data
// ------------------------------------------------------------------------------
// Constants:
    const std::size_t nAPVSamples = 6;
 
    const double dt_APV = 31.44 * Unit::ns;
 
    typedef double apvSampleBaseType;
    typedef std::array<apvSampleBaseType, nAPVSamples> apvSamples; 
    const apvSamples apvTimeBase =
    {{ -dt_APV, 0.0, dt_APV, 2 * dt_APV, 3 * dt_APV, 4 * dt_APV}};
 
// ==============================================================================
// APV25 waveforms
// ------------------------------------------------------------------------------
 
    typedef std::function<double(double)> waveFunction;
 
    inline double w_expo(double t)
    {
      if (t < 0.0) return 0.0;
      else return t * exp(1.0 - t);
    }
 
    inline double w_poly3(double t)
    {
      if (t < 0.0 || t > 1.0)
        return 0.0;
      else
        return 6.75 * t * (1.0 - t) * (1.0 - t);
    }
 
    inline double w_betaprime(double t)
    {
      if (t < 0.0)
        return 0.0;
      else
        return 149.012 * pow(t, 2) * pow(1.0 + t, -10);
    }
 
    class WaveGenerator {
    public:
      WaveGenerator(waveFunction wave = w_betaprime):
        m_samples( {{0, 0, 0, 0, 0, 0}}), m_wave(wave)
      {}
      void setWaveFunction(waveFunction wave) { m_wave  = wave; }
      const apvSamples& operator()(double t0, double tau)
      {
        std::transform(
          apvTimeBase.begin(), apvTimeBase.end(),
          m_samples.begin(),
          [this, t0, tau](double t)->double { return m_wave((t - t0) / tau); }
        );
        return m_samples;
      }
    private:
      apvSamples m_samples; 
      waveFunction m_wave; 
    };
 
 
    // ==============================================================================
    // The 3-samples filter
    // ------------------------------------------------------------------------------
 
    template<typename T>
    inline void zeroSuppress(T& a, double thr)
    {
      std::replace_if(a.begin(), a.end(), std::bind2nd(std::less<double>(), thr), 0.0);
    }
 
// ==============================================================================
// The 3-samples filter
// ------------------------------------------------------------------------------
 
    template<typename T>
    inline bool pass3Samples(const T& a, double thr)
    {
      return (search_n(a.begin(), a.end(), 3, thr, std::greater<double>()) != a.end());
    }
 
// ==============================================================================
// Tau (scale) conversion and encoding
// ------------------------------------------------------------------------------
 
    inline double tau_raw2real(double raw_tau) { return 7.32313 * raw_tau; }
 
    class TauEncoder {
    public:
      TauEncoder(double min_amplitude, double max_amplitude,
                 double min_tau, double max_tau):
        m_minAmplitude(min_amplitude), m_maxAmplitude(max_amplitude),
        m_minTau(min_tau), m_maxTau(max_tau)
      {
        m_ATRatio = (m_maxAmplitude - m_minAmplitude) / (m_maxTau - m_minTau);
      }
      TauEncoder(): TauEncoder(0, 100, 0, 100) {}
 
      void setBounds(double min_amplitude, double max_amplitude,
                     double min_tau, double max_tau)
      {
        m_minAmplitude = min_amplitude;
        m_maxAmplitude = max_amplitude;
        m_minTau = min_tau;
        m_maxTau = max_tau;
        // Re-initialize the range ratio setting
        m_ATRatio = (m_maxAmplitude - m_minAmplitude) / (m_maxTau - m_minTau);
      }
 
      double encodeTau(double tau) const
      { return m_minAmplitude + m_ATRatio * (tau - m_minTau); }
      double decodeTau(double scaledTau) const
      { return m_minTau + 1.0 / m_ATRatio * (scaledTau - m_minAmplitude); }
 
      void print(std::ostringstream& os) const
      {
        os << "TauEncoder object, "
           << "Amplitude range: " << m_minAmplitude << " to " << m_maxAmplitude
           << " Tau range: " << m_minTau << " to " << m_maxTau
           << " ATRatio: " << m_ATRatio << std::endl;
      }
 
    private:
      double m_minAmplitude; 
      double m_maxAmplitude; 
      double m_minTau; 
      double m_maxTau; 
      double m_ATRatio; 
    };
// ------------------------------------------------------------------------------
 
  } // namespace SVD
} // namespace Belle2
 
#endif