Merge pull request #52 from vlvovch/devel

Toward version 1.5.2

Author: vlvovch

CHANGELOG.md

@@ -4,6 +4,23 @@ All notable changes to this project will be documented in this file.
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
+## [Version 1.5.2]
+
+Date: 2025-07-23
+
+This update mainly fixes a couple of bugs in event generator mode (thanks to [A. Ohlson](https://github.com/aohlson)).
+Namely, it fixes a bug causing infinite loop in the event generator when the number of particles in the event is zero.
+This can happen when simulating low-multiplicity proton-proton collisions with local charge conservation.
+It also fixes an issue with initialization of Breit-Wigner distributions in the event generator, which could cause invalid_argument exceptions.
+
+## Bugfixes
+- Fixed a bug in the event generator causing infinite loop when the number of particles in the event is zero
+- Fixed an issue with initialization of Breit-Wigner distributions by properly checking if the width is zero
+
+## Enhancements  
+- Implemented $m \to 0$ limit for quantum statistics calculations using cluster expansion
+- Added additional flag checks for ThermalModelCanonical calculations
+
 ## [Version 1.5.1]
 
 Date: 2025-07-06
@@ -350,7 +367,10 @@ Date: 2018-08-02
 
 **The first public version of Thermal-FIST**
 
+[Version 1.5.2]: https://github.com/vlvovch/Thermal-FIST/compare/v1.5.1...v1.5.2
+
 [Version 1.5.1]: https://github.com/vlvovch/Thermal-FIST/compare/v1.5...v1.5.1
+
 [Version 1.5]: https://github.com/vlvovch/Thermal-FIST/compare/v1.4.2...v1.5
 
 [Version 1.4.2]: https://github.com/vlvovch/Thermal-FIST/compare/v1.4.1...v1.4.2

CMakeLists.txt

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

include/HRGBase/ThermalModelBase.h

@@ -1145,7 +1145,7 @@ namespace thermalfist {
     int  PdgToId(long long pdgid) const { return m_TPS->PdgToId(pdgid); }
 
     /// Reset all flags which correspond to a calculation status
-    void ResetCalculatedFlags();
+    virtual void ResetCalculatedFlags();
 
     /// Whether densities had been calculated with
     /// the CalculateDensities() method

include/HRGBase/ThermalModelCanonical.h

@@ -127,7 +127,7 @@ namespace thermalfist {
      * 
      * \param conserve true -- canonically, false -- grand-canonically
      */
-    virtual void ConserveBaryonCharge(bool conserve = true) { m_BCE = static_cast<int>(conserve); }
+    virtual void ConserveBaryonCharge(bool conserve = true) { m_BCE = static_cast<int>(conserve); m_PartialZCalculated = false;  }
 
     /**
      * \brief Specifies whether the electric charge is treated canonically.
@@ -136,7 +136,7 @@ namespace thermalfist {
      * 
      * \param conserve true -- canonically, false -- grand-canonically
      */
-    virtual void ConserveElectricCharge(bool conserve = true) { m_QCE = static_cast<int>(conserve); }
+    virtual void ConserveElectricCharge(bool conserve = true) { m_QCE = static_cast<int>(conserve); m_PartialZCalculated = false; }
 
     /**
      * \brief Specifies whether the strangeness charge is treated canonically.
@@ -145,7 +145,7 @@ namespace thermalfist {
      * 
      * \param conserve true -- canonically, false -- grand-canonically
      */
-    virtual void ConserveStrangeness(bool conserve = true) { m_SCE = static_cast<int>(conserve); }
+    virtual void ConserveStrangeness(bool conserve = true) { m_SCE = static_cast<int>(conserve); m_PartialZCalculated = false; }
 
     /**
      * \brief Specifies whether the charm charge is treated canonically.
@@ -154,7 +154,7 @@ namespace thermalfist {
      * 
      * \param conserve true -- canonically, false -- grand-canonically
      */
-    virtual void ConserveCharm(bool conserve = true) { m_CCE = static_cast<int>(conserve); }
+    virtual void ConserveCharm(bool conserve = true) { m_CCE = static_cast<int>(conserve); m_PartialZCalculated = false; }
 
     virtual bool IsConservedChargeCanonical(ConservedCharge::Name charge) const;
 
@@ -182,6 +182,27 @@ namespace thermalfist {
      */
     void SetIntegrationIterationsMultiplier(int multiplier) { (multiplier > 0 ? m_IntegrationIterationsMultiplier = multiplier : m_IntegrationIterationsMultiplier = 1); }
 
+    /* 
+     * \brief Reset all flags which correspond to a calculation status
+     */
+    virtual void ResetCalculatedFlags() {
+      ThermalModelBase::ResetCalculatedFlags();
+      ResetPartialZCalculated();
+    }
+    /*
+     * \brief Reset the flags indicating whether the partial partition functions are calculated
+    */
+    void ResetPartialZCalculated() { m_PartialZCalculated = false; m_PartialZCalculatedFlucts = false; }
+
+    /*
+     * \brief Whether the partial partition functions are calculated
+     */
+    bool IsPartialZCalculated() const { return m_PartialZCalculated; }
+    
+    /*
+     * \brief Whether the partial partition functions are calculated with fluctuations
+     */
+    bool IsPartialZCalculatedFlucts() const { return m_PartialZCalculatedFlucts; }
 
     // Override functions begin
 
@@ -310,6 +331,9 @@ namespace thermalfist {
      * \brief Flag indicating whether the analytical calculation of baryon fugacity is used.
      */
     bool m_Banalyt;
+
+    bool m_PartialZCalculated;
+    bool m_PartialZCalculatedFlucts;
   };
 
 } // namespace thermalfist

include/HRGBase/Utility.h

@@ -16,6 +16,7 @@
  */
 
 #include <string>
+#include <map>
 
 namespace thermalfist {
 
@@ -28,6 +29,43 @@ namespace thermalfist {
 
   };
 
+  // Read parameters from a file
+  /**
+   *  Records the parameters from a file into a map.
+   * 
+   *  The file should contain the parameters in the format:
+   *  
+   *  # comment
+   *  var1 val1  # comment
+   *  var2 val2
+   *  ...
+   * 
+   *  where var1, var2, ... are the parameter names and val1, val2, ... are the parameter values.
+   * 
+   *  The function returns the map with the parameters.
+   *  It overrides the parameters in @params with the parameters from the file.
+   * 
+   *  @param filename The name of the file with the parameters.
+   *  @param params The map with the existing/default parameters.
+   * 
+   */
+  std::map<std::string, std::string> ReadParametersFromFile(const std::string& filename, const std::map<std::string, std::string>& params = {});
+
+  /**
+   *  Records the parameters from command line arguments into a map.
+   * 
+   *  --var1 val1 --var2 val2 ...
+   * 
+   *  where var1, var2, ... are the parameter names and val1, val2, ... are the parameter values.
+   * 
+   *  The function writes the parameters to the map.
+   * 
+   *  @param argc The number of command line arguments.
+   *  @param argv The array of command line arguments.
+   *  @param params The map with the existing/default parameters.
+   * 
+   */
+  void ParametersFromArgs(int argc, char* argv[], std::map<std::string, std::string>& params);
 
   // For C99 compatibility
   long long stringToLongLong(const std::string &str);

src/examples/NeutronStars/NeutronStars-CSHRG.cpp

@@ -41,7 +41,7 @@ map<string, double> params = {
 };
 
 // Read the parameters from file
-void ReadParametersFromFile(const std::string &filename)
+void ReadParametersFromFileNS(const std::string &filename)
 {
 	if (filename == "")
 		return;
@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
   string params_file = "";
   if (argc > 1)
     params_file = argv[1];
-  ReadParametersFromFile(params_file);
+  ReadParametersFromFileNS(params_file);
 
   bool include_leptons = (params["include_leptons"] != 0.);
 

src/library/HRGBase/ThermalModelCanonical.cpp

@@ -37,6 +37,8 @@ namespace thermalfist {
     m_modelgce = NULL;
 
     m_Banalyt = false;
+    m_PartialZCalculated = false;
+    m_PartialZCalculatedFlucts = false;
   }
 
 
@@ -119,6 +121,7 @@ namespace thermalfist {
             ind++;
           }
 
+    m_PartialZCalculatedFlucts = computeFluctuations;
   }
 
   void ThermalModelCanonical::SetStatistics(bool stats) {
@@ -590,6 +593,8 @@ Obtained: %lf\n\
     for (size_t iN = 0; iN < m_PartialZ.size(); ++iN) {
       m_Corr[iN] = m_PartialZ[iN] / m_PartialZ[m_QNMap[QuantumNumbers(0, 0, 0, 0)]];
     }
+
+    m_PartialZCalculated = true;
   }
 
   double ThermalModelCanonical::ParticleScaledVariance(int part)
@@ -740,6 +745,11 @@ Obtained: %lf\n\
   }
 
   void ThermalModelCanonical::CalculateFluctuations() {
+    if (!m_PartialZCalculatedFlucts) {
+      CalculateQuantumNumbersRange(true);
+      CalculatePrimordialDensities();
+    }
+      
     CalculateTwoParticleCorrelations();
     CalculateSusceptibilityMatrix();
     CalculateTwoParticleFluctuationsDecays();

src/library/HRGBase/Utility.cpp

@@ -9,6 +9,8 @@
 
 #include <iostream>
 #include <sstream>
+#include <fstream>
+#include <limits>
 
 #include "ThermalFISTConfig.h"
 
@@ -156,6 +158,48 @@ namespace thermalfist {
     return ret * mn;
   }
 
+  std::map<std::string, std::string> ReadParametersFromFile(const std::string& filename, const std::map<std::string, std::string>& params) {
+    std::map<std::string, std::string> ret = params;
+    std::ifstream fin(filename);
+
+    if (!fin.is_open()) {
+      std::cout << "Cannot open parameters file!" << "\n";
+      return ret;
+    }
+
+    std::string var;
+    std::cout << "Reading input parameters from file " << filename << "\n";
+    while (fin >> var) {
+      if (var.size() == 0 || var[0] == '#') {
+        fin.ignore((std::numeric_limits<std::streamsize>::max)(), '\n');
+        continue;
+      }
+
+      std::cout << "Reading input parameter " << var << " = ";
+      std::string val;
+      fin >> val;
+      std::cout << val << std::endl;
+      ret[var] = val;
+    }
+    fin.close();
+    return ret;
+  }
+
+
+  void ParametersFromArgs(int argc, char *argv[], std::map<std::string, std::string> &params)
+  {
+    for(int i = 1; i < argc - 1; i++) {
+      std::string arg = argv[i];
+      if (arg.size() <= 2 || arg[0] != '-' || arg[1] != '-') {
+        continue;
+      }
+      std::string var = arg.substr(2);
+      std::string val = argv[i + 1];
+      params[var] = val;
+      i++;
+    }
+  }
+
   // Time keeping
   // Windows
   #ifdef _WIN32