HypersurfaceSampler.h

Go to the documentation of this file.

#ifndef HYPERSURFACESAMPLER_H
#define HYPERSURFACESAMPLER_H
 
/*
 * Thermal-FIST package
 *
 * Copyright (c) 2021 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
 
#include "HRGBase/SplineFunction.h"
#include "HRGEventGenerator/MomentumDistribution.h"
#include "HRGEventGenerator/RandomGenerators.h"
#include "HRGEventGenerator/SphericalBlastWaveEventGenerator.h"
#include "HRGEventGenerator/CylindricalBlastWaveEventGenerator.h"
 
namespace thermalfist {
 
  struct ParticlizationHypersurfaceElement {
    double tau, x, y, eta;
    double dsigma[4];
    double u[4];
    double T, muB, muQ, muS;
    double edens, rhoB;
    double pi[10];
    double Pi;
    double P;
  };
 
  typedef std::vector<ParticlizationHypersurfaceElement> ParticlizationHypersurface;
 
  namespace RandomGenerators {

    class VolumeElementSampler {
      std::vector<double> m_CumulativeProbabilities;
    public:
      VolumeElementSampler(const ParticlizationHypersurface* Hypersurface = NULL);
      VolumeElementSampler(const std::vector<double>& Weights);
      void FillProbabilities(const ParticlizationHypersurface* Hypersurface);
      void FillProbabilities(const std::vector<double>& Weights);
      void SetProbabilities(const std::vector<double>& CumulativeProbabilities) { m_CumulativeProbabilities = CumulativeProbabilities; }
      int SampleVolumeElement(MTRand& rangen = RandomGenerators::randgenMT) const;
    };
 

    class HypersurfaceMomentumGenerator
      : public ParticleMomentumGenerator
    {
    public:

      static std::vector<double> SamplePhaseSpaceCoordinateFromElement(
        const ParticlizationHypersurfaceElement* elem,
        const ThermalParticle* particle,
        const double& mass = -1.,
        const double& etasmear = 0.,
        const bool shear_correction = false,
        const bool bulk_correction = false,
        const double speed_of_sound_squared = 0.15
      );

      struct HypersurfaceMomentumGeneratorConfiguration {
        double etaSmear;
        bool shearCorrection;
        bool bulkCorrection;
        double speedOfSoundSquared;
        HypersurfaceMomentumGeneratorConfiguration(double etasmear = 0.0, bool shear_correction = false, bool bulk_correction = false, double speed_of_sound_squared = 0.15) :
                etaSmear(etasmear), shearCorrection(shear_correction), bulkCorrection(bulk_correction), speedOfSoundSquared(speed_of_sound_squared) {}
      };

      HypersurfaceMomentumGenerator(
        const ParticlizationHypersurface* hypersurface = NULL,
        const ThermalParticle* particle = NULL,
        const VolumeElementSampler* positionsampler = NULL,
        const HypersurfaceMomentumGeneratorConfiguration& config = HypersurfaceMomentumGeneratorConfiguration());

      virtual ~HypersurfaceMomentumGenerator() {}
 
      double EtaSmear() const { return m_EtaSmear; }
      double Mass() const { return m_Particle->Mass(); }
      bool ShearCorrection() const { return m_ShearCorrection; }
      bool BulkCorrection() const { return m_BulkCorrection; }
      double SpeedOfSoundSquared() const { return m_SpeedOfSoundSquared; }
 
      // Override functions begin
 
      virtual std::vector<double> GetMomentum(double mass = -1.) const;
 
      // Override functions end
 
 
    protected:
 
    private:
      const ParticlizationHypersurface* m_ParticlizationHypersurface;
      const ThermalParticle* m_Particle;
      const VolumeElementSampler* m_VolumeElementSampler;
      double m_EtaSmear;
      bool m_ShearCorrection;
      bool m_BulkCorrection;
      double m_SpeedOfSoundSquared;
    };
 

    class BoostInvariantHypersurfaceMomentumGenerator
      : public ParticleMomentumGenerator
    {
    public:
 

      BoostInvariantHypersurfaceMomentumGenerator(
        const ParticlizationHypersurface* hypersurface = NULL,
        VolumeElementSampler* positionsampler = NULL,
        double Tkin = 0.100, double etamax = 3.0, double mass = 0.938, int statistics = 0, double mu = 0);

      virtual ~BoostInvariantHypersurfaceMomentumGenerator() {}
 
      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:
 
    private:
      const ParticlizationHypersurface* m_ParticlizationHypersurface;
      VolumeElementSampler* m_VolumeElementSampler;
      ThermalMomentumGenerator m_Generator;
      double m_Tkin;
      double m_EtaMax;
      double m_Mass;
    };
 
  }
 
  

  class HypersurfaceEventGenerator : public EventGeneratorBase
  {
  public:

    HypersurfaceEventGenerator(
      const ParticlizationHypersurface* hypersurface = NULL,
      ThermalModelBase* model = NULL,
      double etasmear = 0.0, bool shear_correction = false, bool bulk_correction = false, double speed_of_sound_squared = 0.15) : EventGeneratorBase()
    {
      SetHypersurface(hypersurface);
      SetEtaSmear(etasmear);
      SetRescaleTmu();
      SetShearCorrection(shear_correction);
      SetBulkCorrection(bulk_correction);
      SetSpeedOfSoundSquared(speed_of_sound_squared);
      m_THM = model;
      //SetParameters(hypersurface, model, etasmear);
    }

    HypersurfaceEventGenerator(
            ThermalParticleSystem* TPS,
            const EventGeneratorConfiguration& config = EventGeneratorConfiguration(),
            const ParticlizationHypersurface* hypersurface = NULL,
            const RandomGenerators::HypersurfaceMomentumGenerator::HypersurfaceMomentumGeneratorConfiguration& configMomentumGenerator = RandomGenerators::HypersurfaceMomentumGenerator::HypersurfaceMomentumGeneratorConfiguration());
 

    HypersurfaceEventGenerator(
      ThermalParticleSystem* TPS,
      const EventGeneratorConfiguration& config,
      const ParticlizationHypersurface* hypersurface,
      double etasmear, bool shear_correction = false, bool bulk_correction = false, double speed_of_sound_squared = 0.15);
 
    virtual ~HypersurfaceEventGenerator() {}

 
    //virtual std::pair< std::vector<int>, double > SampleYields() const { exit(1); return std::pair< std::vector<int>, double >(); }
    //virtual SimpleEvent GetEvent(bool PerformDecays = true) const { exit(1); return SimpleEvent(); }
    virtual std::vector<double> GCEMeanYields() const;
    virtual std::vector<double>& GCEMeanYields();

 
    
 
    void SetModel(ThermalModelBase* model) { m_THM = model; m_ParametersSet = false; }
 
    void SetHypersurface(const ParticlizationHypersurface* hypersurface) { m_ParticlizationHypersurface = hypersurface; m_ParametersSet = false; }
 
    void SetMomentumGeneratorConfig(const RandomGenerators::HypersurfaceMomentumGenerator::HypersurfaceMomentumGeneratorConfiguration& config) { m_MomentumGeneratorConfig = config; m_ParametersSet = false; }
    RandomGenerators::HypersurfaceMomentumGenerator::HypersurfaceMomentumGeneratorConfiguration GetMomentumGeneratorConfig() const { return m_MomentumGeneratorConfig; }
 
    void SetEtaSmear(double etaSmear) { m_MomentumGeneratorConfig.etaSmear = etaSmear; m_ParametersSet = false; }
    double GetEtaSmear() const { return m_MomentumGeneratorConfig.etaSmear; }
 
    void SetShearCorrection(bool shear_correction) { m_MomentumGeneratorConfig.shearCorrection = shear_correction; }
    bool GetShearCorrection() { return m_MomentumGeneratorConfig.shearCorrection; }
 
    void SetBulkCorrection(bool bulk_correction) { m_MomentumGeneratorConfig.bulkCorrection = bulk_correction; }
    bool GetBulkCorrection() { return m_MomentumGeneratorConfig.bulkCorrection; }
 
    void SetSpeedOfSoundSquared(double speed_of_sound_squared) { m_MomentumGeneratorConfig.speedOfSoundSquared = speed_of_sound_squared; }
    double GetSpeedOfSoundSquared() { return m_MomentumGeneratorConfig.speedOfSoundSquared; }
 
    void SetRescaleTmu(bool rescale = false, double edens = 0.26);

    //void SetParameters(const ParticlizationHypersurface* hypersurface, ThermalModelBase* model, double etasmear = 0.0);
    //virtual void SetParameters();

    virtual SimpleEvent GetEvent(bool PerformDecays = true) const;
 

    //virtual void CheckSetParameters() { if (!m_ParametersSet) SetParameters(); }
 
  protected:
    void SetMomentumGenerators();

    void ProcessVolumeElements();

    static std::vector<std::vector<double>> CalculateTMuMap(ThermalModelBase* model, double edens, double rhomin = 0.0, double rhomax = 0.27, double drho = 0.001);

    static void RescaleHypersurfaceParametersEdens(ParticlizationHypersurface *hypersurface, ThermalModelBase* model, double edens, double rhomin = 0.0, double rhomax = 0.27, double drho = 0.001, double rhocrit = 0.16);

    //void SetParameters(const ParticlizationHypersurface* hypersurface, ThermalModelBase* model, double etasmear = 0.0);
    virtual void SetParameters();

    const std::vector<double>& FullSpaceYields() const { return m_FullSpaceYields; }
    
    //bool m_ParametersSet;
 
  private:
    const ParticlizationHypersurface* m_ParticlizationHypersurface;
    std::vector<RandomGenerators::VolumeElementSampler> m_VolumeElementSamplers;
    std::vector<double> m_FullSpaceYields;
    double m_Tav;
    std::vector<double> m_Musav;
    bool m_RescaleTmu;
    RandomGenerators::HypersurfaceMomentumGenerator::HypersurfaceMomentumGeneratorConfiguration m_MomentumGeneratorConfig;
    double m_edens;
    std::vector<SplineFunction> m_SplinesTMu;
 
    // Find T and muB = 0 to match energy density
    class BroydenEquationsTen : public BroydenEquations
    {
    public:
      BroydenEquationsTen(ThermalModelBase* model, double edens = 0.5) : BroydenEquations(),
        m_THM(model), m_edens(edens) {
        m_N = 1;
      }
 
      std::vector<double> Equations(const std::vector<double>& x) {
        const double& T = x[0];
 
        m_THM->SetTemperature(T);
        m_THM->SetBaryonChemicalPotential(0.);
        m_THM->SetElectricChemicalPotential(0.);
        m_THM->SetStrangenessChemicalPotential(0.);
        m_THM->CalculatePrimordialDensities();
 
        double en = m_THM->EnergyDensity();
 
        std::vector<double> ret(1, 0);
        ret[0] = (en / m_edens - 1.);
 
        return ret;
      }
    private:
      ThermalModelBase* m_THM;
      double m_edens;
    };
 
    // Find T,muB to match energy and baryon densities
    class BroydenEquationsTmuB : public BroydenEquations
    {
    public:
      BroydenEquationsTmuB(ThermalModelBase* model, double edens = 0.5, double rhoB = 0.16, int zeroMuBmode = 0) : BroydenEquations(),
        m_THM(model), m_edens(edens), m_rhoB(rhoB), m_zeroMuBmode(zeroMuBmode)
      {
        //m_N = 2;
        m_N = 4;
      }
      std::vector<double> Equations(const std::vector<double>& x) {
        const double& T = x[0];
        const double& muB = x[1];
        const double& muS = x[2];
        const double& muQ = x[3];
 
        //m_THM->SetQoverB(0.4);
        //m_THM->ConstrainMuQ(true);
        //m_THM->ConstrainMuS(true);
 
        m_THM->SetTemperature(T);
        if (!m_zeroMuBmode) {
          //m_THM->SetBaryonChemicalPotential(muB);
          //m_THM->ConstrainChemicalPotentials(false);
          m_THM->SetBaryonChemicalPotential(muB);
          m_THM->SetStrangenessChemicalPotential(muS);
          m_THM->SetElectricChemicalPotential(muQ);
          m_THM->CalculatePrimordialDensities();
        }
        else {
          m_THM->SetBaryonChemicalPotential(0.);
          m_THM->SetElectricChemicalPotential(0.);
          m_THM->SetStrangenessChemicalPotential(0.);
          m_THM->CalculatePrimordialDensities();
        }
 
        double en = m_THM->EnergyDensity();
        double rhoB = m_THM->BaryonDensity();
 
        std::vector<double> ret(4, 0);
        ret[0] = (en / m_edens - 1.);
 
        if (!m_zeroMuBmode) {
          ret[1] = rhoB / m_rhoB - 1.;
          ret[2] = m_THM->StrangenessDensity() / m_THM->AbsoluteStrangenessDensity();
          ret[3] = m_THM->ElectricChargeDensity() / m_THM->BaryonDensity() / 0.4 - 1.;
        }
        else {
          ret[1] = 0.;
          ret[2] = 0.;
          ret[3] = 0.;
        }
        return ret;
      }
    private:
      ThermalModelBase* m_THM;
      double m_edens, m_rhoB;
      int m_zeroMuBmode;
    };
 
    // returns (T,muB,muS,muQ) from given energy and baryon densities, assuming Q/B = 0.4, S = 0
    static std::vector<double> MatchEnergyBaryonDensities(ThermalModelBase* model, double edens, double rhoB) {
      if (edens == 0.0)
        return { 0.,0.,0.,0. };
 
      int zeromuBmode = 0;
      if (rhoB == 0.0)
        zeromuBmode = 1;
 
      if (!zeromuBmode) {
        BroydenEquationsTmuB eqs(model, edens, rhoB, zeromuBmode);
        Broyden broydn(&eqs);
        if (model->Parameters().muB == 0.0)
          model->SetBaryonChemicalPotential(0.010);
        //broydn.Solve({ model->Parameters().T, model->Parameters().muB });
        broydn.Solve({ model->Parameters().T, model->Parameters().muB, model->Parameters().muS, model->Parameters().muQ });
      }
      else {
        BroydenEquationsTen eqs(model, edens);
        Broyden broydn(&eqs);
        broydn.Solve({ model->Parameters().T });
      }
 
      return {
        model->Parameters().T,
        model->Parameters().muB,
        model->Parameters().muS,
        model->Parameters().muQ,
        model->Pressure()
      };
    }
  };

  class HypersurfaceEventGeneratorEVHRG : public HypersurfaceEventGenerator
  {
  public:
    HypersurfaceEventGeneratorEVHRG(
      ThermalParticleSystem* TPS,
      const EventGeneratorConfiguration& config = EventGeneratorConfiguration(),
      const ParticlizationHypersurface* hypersurface = NULL,
      double etasmear = 0.0, bool shear_correction = false) : HypersurfaceEventGenerator(TPS, config, hypersurface, etasmear, shear_correction) {
      m_b = 0.0;
      m_rad = 0.0;
      m_MeanB = m_MeanAB = m_VEff = 0.0;
    }

    virtual SimpleEvent GetEvent(bool DoDecays = true) const;
 
    void SetExcludedVolume(double b) { m_b = b; m_ParametersSet = false; }
    void SetBaryonRadius(double r) { m_rad = r; m_ParametersSet = false; }
  protected:
    static double EVHRGWeight(int sampledN, double meanN, double V, double b);
    virtual void SetParameters();
 
    virtual SimpleEvent SampleParticles(const std::vector<int>& yields) const;
 
  private:
 
    std::pair<int, int> ComputeNBNBbar(const std::vector<int>& yields) const;
 
    double m_b, m_rad;
    double m_MeanB, m_MeanAB, m_VEff;
  };

  class BoostInvariantHypersurfaceEventGenerator : public EventGeneratorBase
  {
  public:

    BoostInvariantHypersurfaceEventGenerator(ThermalParticleSystem* TPS = NULL,
      const EventGeneratorConfiguration& config = EventGeneratorConfiguration(),
      double etamax = 0.5,
      const ParticlizationHypersurface* hypersurface = NULL);
 
    virtual ~BoostInvariantHypersurfaceEventGenerator();

    void SetParameters(double etamax, const ParticlizationHypersurface* hypersurface);
 
    double GetEtaMax() const { return m_EtaMax; }
 
  protected:
    void SetMomentumGenerators();
 
  private:
    double m_EtaMax;
    const ParticlizationHypersurface* m_ParticlizationHypersurface;
    RandomGenerators::VolumeElementSampler* m_VolumeElementSampler;
  };
 
} // namespace thermalfist
 
#endif