Merge pull request #21 from vlvovch/devel

Version 1.2

CMakeLists.txt

@@ -7,8 +7,8 @@ project (ThermalFIST)
 
 # The version number.
 set (ThermalFIST_VERSION_MAJOR 1)
-set (ThermalFIST_VERSION_MINOR 1)
-set (ThermalFIST_VERSION_DEVEL 5)
+set (ThermalFIST_VERSION_MINOR 2)
+set (ThermalFIST_VERSION_DEVEL 0)
 
 # configure a header file to pass some of the CMake settings
 # to the source code

include/HRGBase/ThermalModelBase.h

@@ -454,6 +454,12 @@ namespace thermalfist {
     /// by the condition of charm neutrality
     void ConstrainMuC(bool constrain) { m_ConstrainMuC = constrain; }
 
+    /// Sets whether partial chemical equilibrium with additional chemical potentials is used
+    void UsePartialChemicalEquilibrium(bool usePCE) { m_PCE = usePCE; }
+
+    /// Whether partial chemical equilibrium with additional chemical potentials is used
+    bool UsePartialChemicalEquilibrium() { return m_PCE; }
+
     //@{
       /**
        * \brief The entropy per baryon ratio
@@ -1047,6 +1053,8 @@ namespace thermalfist {
     bool m_ConstrainMuS;
     bool m_ConstrainMuC;
 
+    bool m_PCE;
+
     bool m_useOpenMP;
 
     std::vector<double> m_densities;

include/HRGBase/xMath.h

@@ -61,6 +61,11 @@ namespace thermalfist {
     double BesselK1exp(double x);         // modified Bessel function K_1(x), divided by exponential factor
     double BesselKexp(int n, double x);    // integer order modified Bessel function K_n(x), divided by exponential factor
 
+    double BesselI0exp(double x);         // modified Bessel function I_0(x), divided by exponential factor
+    double BesselI1exp(double x);         // modified Bessel function I_1(x), divided by exponential factor
+    double BesselIexp(int n, double x);   // integer order modified Bessel function I_n(x), divided by exponential factor
+
+
     // Note that the functions Gamma and LogGamma are mutually dependent.
     double LogGamma(double);
     double Gamma(double);

include/HRGEventGenerator/CylindricalBlastWaveEventGenerator.h

@@ -29,31 +29,40 @@ namespace thermalfist {
      * \param TPS    A pointer to the particle list
      * \param config Event generator configuration
      * \param Tkin   The kinetic freeze-out temperature (in GeV) 
-     * \param beta   The transverse flow velocity
+     * \param betas  The transverse flow velocity at the surface
      * \param etamax The longitudinal space-time rapidity cut-off
      * \param npow   The power in the transverse flow profile function
      */
     CylindricalBlastWaveEventGenerator(ThermalParticleSystem *TPS, 
                                        const EventGeneratorConfiguration& config, 
                                        double T = 0.120, 
-                                       double beta = 0.5, 
+                                       double betas = 0.5, 
                                        double etamax = 0.5, 
                                        double npow = 1.);
 
     /// \deprecated
     /// \brief Old constructor. Included for backward compatibility.
-    CylindricalBlastWaveEventGenerator(ThermalModelBase *THM, double T = 0.120, double beta = 0.5, double etamax = 0.5, double npow = 1., bool onlyStable = false, EventGeneratorConfiguration::ModelType EV = EventGeneratorConfiguration::PointParticle, ThermalModelBase *THMEVVDW = NULL);
+    CylindricalBlastWaveEventGenerator(ThermalModelBase *THM, double T = 0.120, double betas = 0.5, double etamax = 0.5, double npow = 1., bool onlyStable = false, EventGeneratorConfiguration::ModelType EV = EventGeneratorConfiguration::PointParticle, ThermalModelBase *THMEVVDW = NULL);
 
     ~CylindricalBlastWaveEventGenerator() { }
 
     /// Sets the momentum distribution parameters
-    void SetParameters(double T, double beta, double etamax, double npow = 1.);
+    void SetParameters(double T, double betas, double etamax, double npow = 1.);
+
+    /**
+     * \brief Set the mean transverse flow velocity.
+     *
+     * Surface velocity is calculated as \\beta_s = \\langle \\beta_T \\rangle (2+n)/2 
+     * 
+     * \param betaT  The mean transverse flow velocity
+     */
+    void SetMeanBetaT(double betaT);
 
     /// Sets up the random generators of particle momenta
     /// and resonances masses
     void SetMomentumGenerators();
   private:
-    double m_T, m_Beta, m_EtaMax, m_n;
+    double m_T, m_BetaS, m_EtaMax, m_n;
   };
 
   /// For backward compatibility

include/HRGEventGenerator/EventGeneratorBase.h

@@ -37,7 +37,7 @@ namespace thermalfist {
     /// Enumerates the statistical ensembles 
     enum Ensemble { 
       GCE, ///< Grand-canonical
-      CE,  ///< Caonical
+      CE,  ///< Canonical
       SCE, ///< Strangeness-canonical
       CCE  ///< Charm-canonical
     };
@@ -68,6 +68,15 @@ namespace thermalfist {
 
     /// The matrix of van der Waals attraction coefficients \f$ a_{ij} \f$
     std::vector< std::vector<double> > aij;
+
+    /// Whether partial chemical equilibrium is used
+    bool fUsePCE;
+
+    /// PCE chemical potentials
+    std::vector<double> fPCEChems;
+
+
+    EventGeneratorConfiguration();
   };
 
   /// \brief Base class for generating events with the Thermal Event Generator
@@ -127,6 +136,21 @@ namespace thermalfist {
     /// sampling used for canonical ensemble and/or eigenvolumes.
     static int fCEAccepted, fCETotal;
 
+    /**
+     * \brief Set system volume.
+     * 
+     * Can be used to include volume fluctuations
+     */
+    void SetVolume(double V);
+
+    /**
+     * \brief Rescale precalculated GCE means. 
+     * 
+     * Called when the system volume is changed
+     */
+    void RescaleCEMeans(double Vmod);
+
+
   protected:
     /**
      * \brief Sets the event generator configuration.

include/HRGEventGenerator/MomentumDistribution.h

@@ -151,14 +151,14 @@ namespace thermalfist {
     /**
      * \copydoc MomentumDistributionBase()
      * \param T      The kinetic temperature (in GeV)
-     * \param beta   The transverse flow velocity
+     * \param beta   The transverse flow velocity at the surface
      * \param etamax The longitudinal space-time rapidity cut-off
      * \param npow   The power in the transverse flow profile function
      * \param norm   Whether the momentum distribution should be normalized to unity 
      */
-    SSHDistribution(int pdgid = 0, double mass = 0., double T = 0.100, double beta = 0.5, double etamax = 0.5, double npow = 1., bool norm = false) :
+    SSHDistribution(int pdgid = 0, double mass = 0., double T = 0.100, double betas = 0.5, double etamax = 0.5, double npow = 1., bool norm = false) :
       MomentumDistributionBase(pdgid, mass),
-      m_T(T), m_Beta(beta), m_EtaMax(etamax), m_n(npow)
+      m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow)
     {
       m_NormY = m_NormPt = m_Norm = 1.;
       if (norm) Normalize();
@@ -172,16 +172,16 @@ namespace thermalfist {
      * \brief Set the parameters of the longitudinal blast-wave distribution
      * 
      * \param T      The kinetic temperature (in GeV)
-     * \param beta   The transverse flow velocity
+     * \param betas  The transverse flow velocity at the surface
      * \param etamax The longitudinal space-time rapidity cut-off
      * \param npow   The power in the transverse flow profile function
      * \param mass   Particle mass (in GeV)
      * \param pdgid  Particle PDG code
      * \param norm   Whether the momentum distribution should be normalized to unity 
      */
-    void SetParameters(double T, double beta, double etamax, double npow, double mass, int pdgid = 0, bool norm = true) {
+    void SetParameters(double T, double betas, double etamax, double npow, double mass, int pdgid = 0, bool norm = true) {
       m_T = T;
-      m_Beta = beta;
+      m_BetaS = betas;
       m_EtaMax = etamax;
       m_n = npow;
       m_Mass = mass;
@@ -192,6 +192,17 @@ namespace thermalfist {
       else Initialize();
     }
 
+    /**
+     * \brief Set the mean transverse flow velocity.
+     *
+     * Surface velocity is calculated as \\beta_s = \\langle \\beta_T \\rangle (2+n)/2 
+     * 
+     * \param betaT  The mean transverse flow velocity
+     */
+    void SetMeanBetaT(double betaT) {
+      m_BetaS = (2. + m_n) / 2. * betaT;
+    }
+
     void Normalize();
 
     /// Rapidity distribution at fixed pT
@@ -239,11 +250,11 @@ namespace thermalfist {
 
     double betar(double r) const {
       if (m_n == 1.)
-        return m_Beta * r;
+        return m_BetaS * r;
       else if (m_n == 2.)
-        return m_Beta * r * m_Beta * r;
+        return m_BetaS * r * r;
       else
-        return pow(m_Beta * r, m_n);
+        return m_BetaS * pow(r, m_n);
     }
 
     double rho(double r) const { return atanh(betar(r)); }
@@ -253,7 +264,7 @@ namespace thermalfist {
     double y2Av() const;
 
     double m_T;
-    double m_Beta, m_EtaMax;
+    double m_BetaS, m_EtaMax;
     double m_NormY, m_NormPt, m_Norm;
     double m_n;
     std::vector<double> m_xlag, m_wlag;

include/HRGEventGenerator/RandomGenerators.h

@@ -37,6 +37,59 @@ namespace thermalfist {
     /// \param rangen A Mersenne Twister random number generator to use
     int RandomPoisson(double mean, MTRand &rangen);
 
+    /// \brief Probability of a Skellam distributed random variable with Poisson means
+    ///        mu1 and mu2 to have the value of k.
+    double SkellamProbability(int k, double mu1, double mu2);
+
+
+    /// \brief Generator of a random number from the Bessel distribution (a, nu), nu is integer
+    ///        Uses methods from https://www.sciencedirect.com/science/article/pii/S016771520200055X
+    ///        Used in event generator with exact conservation of charges to
+    ///        generate two Poisson numbers with fixed difference, as described in https://arxiv.org/pdf/1609.01087.pdf
+    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); }
+    };
+
 
     /// \brief Base class for Monte Carlo sampling of particle momenta
     class ParticleMomentumGenerator
@@ -151,13 +204,13 @@ namespace thermalfist {
       /**
        * \brief Construct a new SSHGenerator object
        * \param T      The kinetic temperature (in GeV)
-       * \param beta   The transverse flow velocity
+       * \param betas  The transverse flow velocity at the surface
        * \param etamax The longitudinal space-time rapidity cut-off
        * \param npow   The power in the transverse flow profile function
        * \param mass   Particle mass (in GeV) 
        */
-      SSHMomentumGenerator(double T, double beta, double etamax, double npow, double mass) :m_T(T), m_Beta(beta), m_EtaMax(etamax), m_n(npow), m_Mass(mass) {
-        m_distr = SSHDistribution(0, m_Mass, m_T, m_Beta, m_EtaMax, m_n, false);
+      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();
@@ -170,18 +223,33 @@ namespace thermalfist {
        * \brief Sets the parameters of the distribution
        * 
        * \param T      The kinetic temperature (in GeV)
-       * \param beta   The transverse flow velocity
+       * \param betas  The transverse flow velocity at the surface
        * \param etamax The longitudinal space-time rapidity cut-off
        * \param npow   The power in the transverse flow profile function
        * \param mass   Particle mass (in GeV) 
        */
-      void SetParameters(double T, double beta, double etamax, double npow, double mass) {
+      void SetParameters(double T, double betas, double etamax, double npow, double mass) {
         m_T = T;
-        m_Beta = beta;
+        m_BetaS = betas;
         m_EtaMax = etamax;
         m_Mass = mass;
         m_n = npow;
-        m_distr = SSHDistribution(0, m_Mass, m_T, m_Beta, m_EtaMax, m_n);
+        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n);
+        m_dPt = 0.02;
+        m_dy = 0.05;
+        FixParameters2();
+      }
+
+      /**
+        * \brief Set the mean transverse flow velocity.
+        *
+        * Surface velocity is calculated as \\beta_s = \\langle \\beta_T \\rangle (2+n)/2 
+        * 
+        * \param betaT  The mean transverse flow velocity
+        */
+      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();
@@ -222,7 +290,7 @@ namespace thermalfist {
 
       std::pair<double, double> GetRandom2() const;
 
-      double m_T, m_Beta, m_EtaMax, m_n, m_Mass;
+      double m_T, m_BetaS, m_EtaMax, m_n, m_Mass;
       double m_MaxY, m_MaxPt;
       SSHDistribution m_distr;
       SplineFunction m_dndpt;

include/HRGEventGenerator/SimpleEvent.h

@@ -23,14 +23,22 @@ namespace thermalfist {
     /// Log of the event weight factor
     double logweight;
 
-    /// Vector of all particles in the event
+    /// Vector of all final particles in the event
     std::vector<SimpleParticle> Particles;
 
+     /// Vector of all particles which ever appeared in the event (including those that decay)
+    std::vector<SimpleParticle> AllParticles;
+
     /// Default constructor, empty event
-    SimpleEvent() { Particles.resize(0); weight = 1.; }
+    SimpleEvent() { Particles.resize(0); AllParticles.resize(0); weight = 1.; }
 
     /// Writes the event to an output file stream
     void writeToFile(std::ofstream & fout, int eventnumber = 1);
+
+    /// Rapidity boost for all particles
+    void RapidityBoost(double dY);
+
+    static SimpleEvent MergeEvents(const SimpleEvent &evt1, const SimpleEvent &evt2);
   };
 
 } // namespace thermalfist