GaussianDip.hh

Go to the documentation of this file.

#ifndef NPSTAT_GAUSSIANDIP_HH_
#define NPSTAT_GAUSSIANDIP_HH_
 
/*!
// \file GaussianDip.hh
//
// \brief "Gaussian dip" function, intended for use in generating
//        random numbers with acceptance-rejection
//
// Author: I. Volobouev
//
// June 2022
*/
 
#include <cmath>
#include <stdexcept>
 
#include "npstat/nm/SimpleFunctors.hh"
 
namespace npstat {
    class GaussianDip : public Functor1<double, double>
    {
    public:
        inline GaussianDip(const double mu, const double sigma, const double a)
            : mu_(mu), sigma_(sigma), amplitude_(a)
        {
            if (sigma_ <= 0.0) throw std::invalid_argument(
                "In npstat::GaussianDip constructor: sigma must be positive");
            if (amplitude_ < -1.0) throw std::invalid_argument(
                "In npstat::GaussianDip constructor: amplitude out of range");
        }
 
        inline virtual ~GaussianDip() {}
 
        inline double operator()(const double& x) const
        {
            const double delta = (x - mu_)/sigma_;
            const double f = 1.0 + amplitude_*std::exp(-delta*delta/2.0);
            if (amplitude_ > 0.0)
                return 1.0/f;
            else
                return f;
        }
 
    private:
        double mu_;
        double sigma_;
        double amplitude_;
    };
}
 
#endif // NPSTAT_GAUSSIANDIP_HH_