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>
#include <framework/utilities/MathHelpers.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 {
        // Replace pow(1.+t, -10) by a division and a sequential calculation of the tenth power for improved execution time
        return 149.012 * square(t) / (pow4(1.0 + t) * pow4(1.0 + t) * square(1.0 + t));
      }
    }

    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] * square(t) + p[6] * cube(t)) * 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] * square(t) + p[6] * cube(t)) * 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