RandomGenerators.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2015-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef RANDOMGENERATORS_H
#define RANDOMGENERATORS_H

#include <cmath>

#include "MersenneTwister.h"
#include "HRGEventGenerator/MomentumDistribution.h"
#include "HRGBase/ThermalParticle.h"

namespace thermalfist {

  namespace RandomGenerators {

    extern MTRand randgenMT;

    void SetSeed(const unsigned int seed);

    int RandomPoisson(double mean);

    int RandomPoisson(double mean, MTRand &rangen);

    double SkellamProbability(int k, double mu1, double mu2);


    struct BesselDistributionGenerator
    {
      static double pn(int n, double a, int nu);

      //static double R(double x, int nu) { return xMath::BesselI(nu + 1, x) / xMath::BesselI(nu, x); }
      static double R(double x, int nu);

      static double mu(double a, int nu) { return a * R(a, nu) / 2; }

      static double chi2(double a, int nu);

      static int    m(double a, int nu) { return static_cast<int>((sqrt(a*a+static_cast<double>(nu)*nu) - nu)/2.); }

      static double sig2(double a, int nu);

      static double Q2(double a, int nu);

      static int RandomBesselPoisson(double a, int nu, MTRand &rangen);

      static int RandomBesselPoisson(double a, int nu) { return RandomBesselPoisson(a, nu, randgenMT); }

      static double pmXmOverpm(int X, int tm, double a, int nu);

      static int RandomBesselDevroye1(double a, int nu, MTRand &rangen);

      static int RandomBesselDevroye1(double a, int nu) { return RandomBesselDevroye1(a, nu, randgenMT); }

      static int RandomBesselDevroye2(double a, int nu, MTRand &rangen);

      static int RandomBesselDevroye2(double a, int nu) { return RandomBesselDevroye2(a, nu, randgenMT); }

      static int RandomBesselDevroye3(double a, int nu, MTRand &rangen);

      static int RandomBesselDevroye3(double a, int nu) { return RandomBesselDevroye3(a, nu, randgenMT); }

      static int RandomBesselNormal(double a, int nu, MTRand &rangen);

      static int RandomBesselNormal(double a, int nu) { return RandomBesselNormal(a, nu, randgenMT); }

      static int RandomBesselCombined(double a, int nu, MTRand &rangen);

      static int RandomBesselCombined(double a, int nu) { return RandomBesselCombined(a, nu, randgenMT); }
    };


    class ParticleMomentumGenerator
    {
    public:
      ParticleMomentumGenerator() {}

      virtual ~ParticleMomentumGenerator() { }

      virtual std::vector<double> GetMomentum(double mass = -1.) const = 0;
    };



    class ThermalMomentumGenerator
    {
    public:
      ThermalMomentumGenerator(double mass = 0.938, int statistics = 0, double T = 0.100, double mu = 0.) :
        m_Mass(mass), m_T(T), m_Mu(mu), m_Statistics(statistics)
      {
        //FixParameters();
        m_Max = ComputeMaximum(m_Mass);
      }

      double GetP(double mass = -1.) const;

    private:
      double g(double x, double mass = -1.) const;

      double ComputeMaximum(double mass) const;

      void FixParameters();

      double m_Mass, m_T, m_Mu;
      int m_Statistics;

      double m_Max;
    };



    class SiemensRasmussenMomentumGenerator
      : public ParticleMomentumGenerator
    {
    public:

      SiemensRasmussenMomentumGenerator() { }

      SiemensRasmussenMomentumGenerator(double T, double beta, double mass) :m_T(T), m_Beta(beta), m_Mass(mass) {
        m_Gamma = 1. / sqrt(1. - m_Beta * m_Beta);
        FixParameters();
      }

      ~SiemensRasmussenMomentumGenerator() { }

      void SetParameters(double T, double beta, double mass) {
        m_T = T;
        m_Beta = beta;
        m_Mass = mass;
        m_Gamma = 1. / sqrt(1. - m_Beta * m_Beta);
        FixParameters();
      }

      double GetBeta() const { return m_Beta; }
      double GetMass() const { return m_Mass; }

      // Override functions begin

      virtual std::vector<double> GetMomentum(double mass = -1.) const;

      // Override functions end


    private:
      double w(double p) const {
        return sqrt(p*p + m_Mass * m_Mass);
      }

      double alpha(double p) const {
        return m_Gamma * m_Beta * p / m_T;
      }

      double g(double x, double mass = -1.) const;
      double g2(double x, double tp, double mass = -1.) const;

      void FixParameters();

      double GetRandom(double mass = -1.) const;

      double m_T;
      double m_Beta;
      double m_Mass;
      double m_Gamma;
      double m_Norm;
      double m_Max;
    };


    class SiemensRasmussenMomentumGeneratorGeneralized
      : public SiemensRasmussenMomentumGenerator
    {
    public:
      SiemensRasmussenMomentumGeneratorGeneralized() { }

      SiemensRasmussenMomentumGeneratorGeneralized(double T, double beta, double mass, int statistics = 0, double mu = 0) :
        SiemensRasmussenMomentumGenerator(T, beta, mass),
        m_Generator(mass, statistics, T, mu)
      {

      }

      // Override functions begin

      virtual std::vector<double> GetMomentum(double mass = -1.) const;

      // Override functions end

    private:
      ThermalMomentumGenerator m_Generator;
    };



    class BoostInvariantMomentumGenerator
      : public ParticleMomentumGenerator
    {
    public:
      BoostInvariantMomentumGenerator(BoostInvariantFreezeoutParametrization* FreezeoutModel = NULL,
        double Tkin = 0.100, double etamax = 3.0, double mass = 0.938, int statistics = 0, double mu = 0);

      virtual ~BoostInvariantMomentumGenerator();

      double EtaMax() const { return m_EtaMax; }
      double Mass() const { return m_Mass; }

      // Override functions begin

      virtual std::vector<double> GetMomentum(double mass = -1.) const;

      // Override functions end

    protected:
     virtual double GetRandomZeta(MTRand& rangen = RandomGenerators::randgenMT) const;

    private:
      BoostInvariantFreezeoutParametrization* m_FreezeoutModel;
      ThermalMomentumGenerator m_Generator;
      double m_Tkin;
      double m_EtaMax;
      double m_Mass;
    };


    class BreitWignerGenerator
    {
    public:

      BreitWignerGenerator() { }

      BreitWignerGenerator(double M, double gamma, double mthr) : m_M(M), m_Gamma(gamma), m_Mthr(mthr) { FixParameters(); }
      
      ~BreitWignerGenerator() { }

      void SetParameters(double M, double gamma, double mthr);

      // Generates random resonance mass m from relativistic Breit-Wigner distribution
      double GetRandom() const;

    private:
      double f(double x) const;

      void FixParameters();

      double m_M;
      double m_Gamma;
      double m_Mthr;
      double m_Norm;
      double m_Max;
    };


    class ThermalBreitWignerGenerator
    {
    public:
      ThermalBreitWignerGenerator() { }

      ThermalBreitWignerGenerator(ThermalParticle *part, double T, double Mu) : m_part(part), m_T(T), m_Mu(Mu) { FixParameters(); }
      
      virtual ~ThermalBreitWignerGenerator() { }

      void SetParameters(ThermalParticle *part, double T, double Mu);

      double GetRandom() const;

    protected:
      virtual void FixParameters();

      virtual double f(double M) const;

      ThermalParticle *m_part;
      double m_T, m_Mu;
      double m_Xmin, m_Xmax;
      double m_Max;
    };

    // Class for generating mass of resonance in accordance with its fixed Breit-Wigner distribution multiplied by the Boltzmann factor
    class ThermalEnergyBreitWignerGenerator
      : public ThermalBreitWignerGenerator
    {
    public:

      ThermalEnergyBreitWignerGenerator() : ThermalBreitWignerGenerator() { }

      ThermalEnergyBreitWignerGenerator(ThermalParticle *part, double T, double Mu) :
        ThermalBreitWignerGenerator(part, T, Mu) {
        FixParameters();
      }

      ~ThermalEnergyBreitWignerGenerator() { }

    protected:
      virtual void FixParameters();
      
      virtual double f(double M) const;
    };
  }

} // namespace thermalfist


namespace thermalfist {

  namespace RandomGenerators {

    // Class for generating momentum of particle with mass m in accordance with Schnedermann-Sollfrank-Heinz formula with temperature T, transverse flow beta, and longitudinal flow etamax
    class SSHMomentumGenerator
      : public ParticleMomentumGenerator
    {
    public:
      SSHMomentumGenerator() { }

      SSHMomentumGenerator(double T, double betas, double etamax, double npow, double mass) :m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow), m_Mass(mass) {
        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n, false);
        m_dPt = 0.02;
        m_dy = 0.05;
        FixParameters2();

      }

      ~SSHMomentumGenerator() { }

      void SetParameters(double T, double betas, double etamax, double npow, double mass) {
        m_T = T;
        m_BetaS = betas;
        m_EtaMax = etamax;
        m_Mass = mass;
        m_n = npow;
        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n);
        m_dPt = 0.02;
        m_dy = 0.05;
        FixParameters2();
      }

      void SetMeanBetaT(double betaT) {
        m_BetaS = (2. + m_n) / 2. * betaT;
        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n);
        m_dPt = 0.02;
        m_dy = 0.05;
        FixParameters2();
      }

      double GetTkin() const { return m_T; }
      double GetBetaSurface() const { return m_BetaS; }
      double GetMass() const { return m_Mass; }
      double GetNPow() const { return m_n; }
      double GetEtaMax() const { return m_EtaMax; }

      // Override functions begin

      std::vector<double> GetMomentum(double mass = -1.) const;

      // Override functions end

    private:
      double w(double p) const {
        return sqrt(p * p + m_Mass * m_Mass);
      }

      double g(double x) const {
        return m_distr.dndpt(-log(x)) / x;
      }

      double g2(double x) const {
        return m_dndpt.f(x);
      }

      void FixParameters();

      void FixParameters2();

      // Finds the maximum of g(x) using the ternary search
      void FindMaximumPt();

      void FindMaximumY(double pt);

      void FindMaximumY2(double pt);

      // Generates random pt and y
      std::pair<double, double> GetRandom(double mass = -1.);

      std::pair<double, double> GetRandom2(double mass = -1.) const;

      double m_T, m_BetaS, m_EtaMax, m_n, m_Mass;
      double m_MaxY, m_MaxPt;
      SSHDistribution m_distr;
      SplineFunction m_dndpt;
      std::vector<SplineFunction> m_dndy;
      std::vector<double> m_MaxYs;
      double m_dPt;
      double m_dy;
    };

  }

} // namespace thermalfist

#endif