SVDSimulationTools.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_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 waveform shapes
// ------------------------------------------------------------------------------
 
    typedef std::function<double(double)> WaveformShape;
 
    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);
    }
 
    inline double w_adjacentU(double t)
    {
 
      double p[8];
 
      // waveform parameters
      p[0] = -7.379;
      p[1] = -0.3594;
      p[2] =  0.07564;
      p[3] =  7.391;
      p[4] = -0.08004;
      p[5] = -2.139e-5;
      p[6] = -1.781e-7;
      p[7] =  0.01827;
 
      if (t < 0.0)
        return 0.0;
      else
        return (p[0] + p[1] * t) * exp(-p[2] * t) + (p[3] + p[4] * t + p[5] * pow(t, 2) + p[6] * pow(t, 3)) * exp(-p[7] * t);
    }
 
    inline double w_adjacentV(double t)
    {
 
      double p[8];
 
      // waveform parameters
      p[0] =  0.2043;
      p[1] = -0.0197;
      p[2] =  0.01448;
      p[3] = -0.2076;
      p[4] =  0.0387;
      p[5] =  0.01055;
      p[6] = -4.03e-6;
      p[7] =  0.07193;
 
      if (t < 0.0)
        return 0.0;
      else
        return (p[0] + p[1] * t) * exp(-p[2] * t) + (p[3] + p[4] * t + p[5] * pow(t, 2) + p[6] * pow(t, 3)) * exp(-p[7] * t);
    }
 
 
    class WaveGenerator {
    public:
      WaveGenerator(WaveformShape waveform = w_betaprime):
        m_samples({{0, 0, 0, 0, 0, 0}}), m_waveform(waveform)
      {}
      void setWaveform(WaveformShape waveform) { m_waveform = waveform; }
      const apvSamples& operator()(double t0, double tau)
      {
        std::transform(
          apvTimeBase.begin(), apvTimeBase.end(),
          m_samples.begin(),
          [this, t0, tau](double t)->double { return m_waveform((t - t0) / tau); }
        );
        return m_samples;
      }
    private:
      apvSamples m_samples; 
      WaveformShape m_waveform; 
    };
 
 
    // ==============================================================================
    // The 3-samples filter
    // ------------------------------------------------------------------------------
 
    template<typename T>
    inline void zeroSuppress(T& a, double thr)
    {
      std::replace_if(a.begin(), a.end(), std::bind(std::less<double>(), std::placeholders::_1, 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