diff --git a/.travis.yml b/.travis.yml index 51a1424..7c0524f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,11 +1,10 @@ language: cpp compiler: g++ sudo: required -dist: trusty +dist: bionic before_install: - - sudo add-apt-repository -y ppa:ubuntu-toolchain-r/test - - sudo apt-get update -qq + - sudo apt-get update install: - sudo apt-get install g++-6 cmake qt5-default diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..b9e69bf --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,194 @@ +# Changelog + +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.3.1] + +Date: 2021-05-26 + +Minor enhancements and bugfixes + +## Enhancements +- GUI: Output moments in addition to susceptiblities in the correlations dialog, output nudyn and D measure in equation of state tab +- GUI: Show lattice data for ratios of quantities for the equation of state tab in the GUI, where available +- GUI: Display mean pT in the event generator tab +- Event generator: Sampling of the space-time coordinates in addition to the momenta (*not yet supported for the Siemens-Rasmussen model*) +- Event generator: Write the output in the format tailor for UrQMD afterburner from https://github.com/jbernhard/urqmd-afterburner +- Event generator: Output particle energy (p0) by default when writing to file + + +## Bugfixes +- Fix the flag setting in ThermalModelPCE for the Saha equation mode and the freeze-out of long-lived resonances +- GUI: Fix the output of strongly intensive quantities for net-particle correlations +- Event generator: Correct resonance mass-energy rescaling when masses of decay products exceed the mass of a resonance +- Add contributions of the K0S decays to the inclusive (weak decay included) yields of pions. Note that decay contributions from K0L are not included + +## [Version 1.3] + +Date: 2020-11-07 + +Version 1.3 contains a number of new features, improvements, and bugfixes + +## Partial chemical equilibrium (PCE) + +This release implements HRG in partial chemical equilibrium (PCE) [Bebie et al, [Nucl. Phys. B **378**, 95 (1992)](https://inspirehep.net/literature/316591)], as appropriate for modeling the hadronic phase of heavy-ion collisions. PCE is implemented as a new **ThermalModelPCE** class, and has been used in two recent papers, [Phys. Lett. B **800**, 135131 (2020)](https://inspirehep.net/literature/1726464) and [Phys. Rev. C **102**, 024909 (2020)](https://inspirehep.net/literature/1751960) + +See [PCE-Saha-LHC.cpp](./src/examples/PCE/PCE-Saha-LHC.cpp) for a usage example and [this link](https://fias.uni-frankfurt.de/~vovchenko/project/thermal-fist/doc/classthermalfist_1_1_thermal_model_p_c_e.html) for documentation. + +The PCE-HRG model is available in the thermal fitting routine, with Tkin being an extra fitting parameter + +The PCE-HRG is also accessible for use in the graphical user interface program through a button `PCE/Saha/Other...` + +## Particle lists + +- The default particle list is now based on the 2020 edition of the PDG listing. As before, the list contains only hadrons and resonances with an established status. Changes compared to the old 2014 PDG list are mild. + +- Version 1.3 now fully supports charm. The charmed hadrons are available in an extended input list file [list-withcharm.dat](./input/list/PDG2020/list-withcharm.dat) + +- The excited nuclei with mass number up to A = 5 have been added with this release, they are available in an extended input list file [list-withexcitednuclei.dat](./input/list/PDG2020/list-withexcitednuclei.dat). The reference for this list of excited nuclei is [Phys. Lett. B **809**, 135746 (2020)](https://inspirehep.net/literature/1790683) + +- The hadron list converted from [SMASH](https://smash-transport.github.io/) transport model (version 1.8) is available [here](./input/list/SMASH-1.8) + +## Other changes + +- Mixed-canonical ensembles in Monte Carlo event generator +- Reworked cylindrical blast-wave event generator with much faster initialization +- Cracow freeze-out generator +- Fermi-Dirac and Bose-Einstein momentum samplers for the event generator with quantum statistics +- Correct decay kinematics involving leptons and photons, the decay leptons/photons, if any, are included in event generator output +- Option for modular loading of particle lists +- Functions for computing net-particle/net-charge correlators (see documentation) +- Leptons and photons will appear in decays of some hadrons with correct kinematics, even if they are not included in the particle list + +## Changes to the graphical user interface +- Calculation of matrices of various 2nd-order correlators (accessible on the `Thermal model` tab via `Correlations...` button +- For quantum van der Waals and X-terms excluded volume models: the EV/vdW interactions by default are now only for baryon-baryon and antibaryon-antibaryon pairs, in order to make it consistent with the literature +- Enable PCE mode via `PCE/Saha/Other...`, then possible to fit Tkin in the thermal fits mode +- Option for quantum statistical momentum distribution in the event generator mode (note that multiplicity distributions are still Poisson!) +- Mixed-canonical ensembles in the event generator mode +- Cracow model in the event generator mode + +## Bugfixes +- Fix for Bose-Einstein condensation issue in event generator with finite resonance widths +- Fix for Fermi-Dirac integrals in case of massless particles +- Fix for decay properties of some charmed hadrons +- Fix for momentum sampling in event generators in case of zero collective/radial flow + +## [Version 1.2] + +Date: 2019-07-22 + +**The version of the code published in [Computer Physics Communications](https://doi.org/10.1016/j.cpc.2019.06.024)** + +### Changes +- Improved performace of the event generator. Exact baryon number conservation now taken into account using [Devroye's method](https://doi.org/10.1016/S0167-7152(02)00055-X) of sampling the Bessel distribution +- Possibility to set volume on an event-by-event basis in the event generator +- The beta parameter in cylindrical blast-wave model now corresponds to the surface velocity parameter \beta_s +- Mixed-canonical ensembles are now available through GUI (use the canonical ensemble and conservation laws dialogs) + +### Bugfixes +- Mixed-canonical ensembles now work at non-zero densities of conserved charges +- Correct analytic calculation of fluctuations in mixed-canonical ensembles +- Small violation of electric charge conservation in event generator within full canonical ensemble +- Proxy net charge susceptibility now corresponds again to net-charge instead of net-pion number. This fixes the output of the cpc4 example, a bug introduced in the pre-release of version 1.2 + +## [Version 1.1] + +Date: 2019-01-25 + +- A pdg code is now held in a 64-bit integer (was 32-bit) +- A pdg code with a trailing zero in thermal fit/yields applications can now be interpreted as a sum of yields of particles and antiparticles, which correspond to the pdg code without this trailing zero +- The fit procedure is improved to catch possible issues with a Bose-Einstein condensation +- cpc2 and cpc3 examples slightly reworked and optimized, calculation of fit number of dof in cpc3 corrected +- Cosmetic changes to the code, mainly in response to compiler warnings + +## [Version 1.0] + +Date: 2019-01-16 + +Matches the code released with the [arXiv documentation paper](https://arxiv.org/abs/1901.05249v1) + +- Charm particles +- Possibility to set seed in the random generator + +## [Version 0.9] + +Date: 2019-01-09 + +**Library** + +- Possibility to constrain baryochemical potential by entropy-per-baryon ratio +- Support for charm-canonical ensemble in event generator +- Code structure optimization +- Doxygen documentation + +**GUI frontend** + +- More convenient HRG model and thermal fit configuration in GUI +- Thermal fit plots +- New equation of state tab for calculating temperature dependencies of various observables, +- Possibility to increase/decrease font size + +**Other** + +- More files with various experimental data +- New examples + +## [Version 0.8] + +Date: 2018-12-26 + +#### Library +- New, more human-readable format of decays input file (old one still supported) +- Broyden method related routines refactored. Non-linear equations are now solved when needed using the generic routines, resulting in a cleaner code +- Fix for calculations of the ideal gas functions in the case of zero particle mass +- Decay feeddown contributions can now be separated into strong/electromagnetic/weak decays +- Fixed bug regarding light quark content calculations for hypernuclei (only relevant if both chemical non-eq. of light quarks and light hypernuclei are considered) +- All third-party libraries/code moved into the `thirdparty` folder +- Resonance widths: eBW scheme with constant branching ratios added (very similar results to the standard eBW scheme but much faster) + +#### GUI +- New tab for editing the particle list on the fly +- Tool-tips +- Possibility to do calculations in the Thermal Model tab using the parameters extracted from a thermal fit in the other tab +- *About Thermal-FIST* dialog window +- Menu, with an option to increase/decrease font size +- Splash screen, app icon (Windows only) + +## [Version 0.7] + +Date: 2018-12-10 + +- eBW scheme now also implemented in event generator +- Quantum statistics by default now computed using quadratures (was cluster expansion) +- Zero-width scheme is now the default one (was energy-independent Breit-Wigner) +- Spectral functions and (energy-dependent) branching ratios visualization in gui +- Some improvements in numerics at low temperatures + +## [Version 0.6] + +Date: 2018-08-02 + +**The first public version of Thermal-FIST** + +[Version 1.3.1]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v1.3.1 + +[Version 1.3]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v1.3 + +[Version 1.2]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v1.2 + +[Version 1.1]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v1.1 + +[Version 1.0]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v1.0 + +[Version 0.9]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v0.9 + +[Version 0.8]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v0.8 + +[Version 0.7]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v0.7 + +[Version 0.6]: https://github.com/vlvovch/Thermal-FIST/releases/tag/v0.6 + diff --git a/CMakeLists.txt b/CMakeLists.txt index 18dfc51..f6d23ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,7 @@ project (ThermalFIST) # The version number. set (ThermalFIST_VERSION_MAJOR 1) set (ThermalFIST_VERSION_MINOR 3) -set (ThermalFIST_VERSION_DEVEL 0) +set (ThermalFIST_VERSION_DEVEL 1) # configure a header file to pass some of the CMake settings # to the source code diff --git a/include/HRGBase/SplineFunction.h b/include/HRGBase/SplineFunction.h index 9e2224d..3b8006c 100644 --- a/include/HRGBase/SplineFunction.h +++ b/include/HRGBase/SplineFunction.h @@ -1,12 +1,4 @@ /* -/* -/* -/* -/* -/* -/* -/* -/* /* * GNU General Public License (GPLv3 or later) */ #ifndef SPLINEFUNCTION_H diff --git a/include/HRGBase/ThermalModelBase.h b/include/HRGBase/ThermalModelBase.h index b35cccb..ba20641 100644 --- a/include/HRGBase/ThermalModelBase.h +++ b/include/HRGBase/ThermalModelBase.h @@ -956,6 +956,7 @@ namespace thermalfist { /// in accordance with the 0-based index of each specie. /// PdgToId() maps PDG ID to the 0-based index. const std::vector& Densities() const { return m_densities; } + std::vector& Densities() { return m_densities; } /// A vector with total particle number densities, /// which include the feeddown contribution in accordance diff --git a/include/HRGBase/ThermalParticleSystem.h b/include/HRGBase/ThermalParticleSystem.h index 9a74d9f..855fb5c 100644 --- a/include/HRGBase/ThermalParticleSystem.h +++ b/include/HRGBase/ThermalParticleSystem.h @@ -560,6 +560,15 @@ namespace thermalfist { void Initialize(const std::string& InputFile = "", const std::string& DecayFile = "", bool GenAntiP = true, double mcut = -1.); + void FinalizeListLoad(); + + void FinalizeDecaysLoad(); + + /// Check is the particle list is in the iSS format + bool CheckListIsiSS(const std::string &filename); + + /// Load particle list from the iSS sampler format used at https://github.com/chunshen1987/iSS + void LoadListiSS(const std::string& filename, const std::set& flags = std::set(), double mcut = 1.e9); private: std::vector m_Particles; diff --git a/include/HRGEV/ThermalModelEVCrossterms.h b/include/HRGEV/ThermalModelEVCrossterms.h index 0ac9e83..00c5b00 100644 --- a/include/HRGEV/ThermalModelEVCrossterms.h +++ b/include/HRGEV/ThermalModelEVCrossterms.h @@ -98,6 +98,11 @@ namespace thermalfist { // Override functions end + + const std::vector< std::vector >& EVComponentIndices() const { return m_EVComponentIndices; } + virtual double DeltaMu(int i) const { return MuShift(i); } + const std::vector< std::vector >& VirialMatrix() const { return m_Virial; } + protected: /** * \brief Solves the system of transcdental equations @@ -197,6 +202,11 @@ namespace thermalfist { double m_Pressure; /**< The (solved) total pressure */ double m_TotalEntropyDensity; /**< The (solved) entropy pressure */ + + std::vector m_MapToEVComponent; + std::vector m_MapFromEVComponent; + std::vector< std::vector > m_EVComponentIndices; + private: class BroydenEquationsCRS : public BroydenEquations { diff --git a/include/HRGEventGenerator/CylindricalBlastWaveEventGenerator.h b/include/HRGEventGenerator/CylindricalBlastWaveEventGenerator.h index fa196e5..cf74130 100644 --- a/include/HRGEventGenerator/CylindricalBlastWaveEventGenerator.h +++ b/include/HRGEventGenerator/CylindricalBlastWaveEventGenerator.h @@ -38,7 +38,8 @@ namespace thermalfist { double T = 0.120, double betas = 0.5, double etamax = 0.5, - double npow = 1.); + double npow = 1., + double Rperp = 6.5); /// \deprecated /// \brief Old constructor. Included for backward compatibility. @@ -66,8 +67,11 @@ namespace thermalfist { double GetBetaSurface() const { return m_BetaS; } double GetNPow() const { return m_n; } double GetEtaMax() const { return m_EtaMax; } + double GetRperp() const { return m_Rperp; } private: - double m_T, m_BetaS, m_EtaMax, m_n; + double GetVeffIntegral() const; + + double m_T, m_BetaS, m_EtaMax, m_n, m_Rperp; }; /// For backward compatibility diff --git a/include/HRGEventGenerator/EventGeneratorBase.h b/include/HRGEventGenerator/EventGeneratorBase.h index f7c318d..589f074 100644 --- a/include/HRGEventGenerator/EventGeneratorBase.h +++ b/include/HRGEventGenerator/EventGeneratorBase.h @@ -32,6 +32,9 @@ namespace thermalfist { return os.str(); } + /// Lorentz boost + std::vector LorentzBoost(const std::vector& fourvector, double vx, double vy, double vz); + /// \brief Structure containing the thermal event generator configuration. struct EventGeneratorConfiguration { /// Enumerates the statistical ensembles @@ -133,7 +136,7 @@ namespace thermalfist { * The first element is a vector of the sampled yields. * The second element is the weight. */ - std::pair< std::vector, double > SampleYields() const; + virtual std::pair< std::vector, double > SampleYields() const; /** * \brief Samples the momenta of the particles and returns the sampled list of particles as an event. @@ -168,7 +171,12 @@ namespace thermalfist { */ static SimpleEvent PerformDecays(const SimpleEvent& evtin, ThermalParticleSystem* TPS); - + /** + * \brief The grand-canonical mean yields. + * + * \return std::vector The computed grand-canonical mean yields. + */ + virtual std::vector GCEMeanYields() const; /// Helper variable to monitor the Acceptance rate of the rejection /// sampling used for canonical ensemble and/or eigenvolumes. @@ -195,6 +203,7 @@ namespace thermalfist { double ComputeWeight(const std::vector& totals) const; + double ComputeWeightNew(const std::vector& totals) const; protected: /** diff --git a/include/HRGEventGenerator/RandomGenerators.h b/include/HRGEventGenerator/RandomGenerators.h index 2c3c39a..9d2120b 100644 --- a/include/HRGEventGenerator/RandomGenerators.h +++ b/include/HRGEventGenerator/RandomGenerators.h @@ -104,7 +104,9 @@ namespace thermalfist { /// Samples the 3-momentum of a particle /// \param mass The mass of a particle. If negative value provided, defaults to the pole/vacuum mass /// \return std::vector A vector containing the sampled - /// \f$p_x\f$, \f$p_y\f$, \f$p_z\f$ components of the three-momentum + /// \f$p_x\f$, \f$p_y\f$, \f$p_z\f$ components of the three-momentum, + /// and, additionally, the space-time Cartesian coordinates \f$r_0\f$, \f$r_x\f$, \f$r_y\f$, \f$r_z\f$ + /// in the collision center-of-mass frame virtual std::vector GetMomentum(double mass = -1.) const = 0; }; @@ -149,7 +151,7 @@ namespace thermalfist { double ComputeMaximum(double mass) const; - void FixParameters(); + //void FixParameters(); double m_Mass, m_T, m_Mu; int m_Statistics; diff --git a/include/HRGEventGenerator/SimpleEvent.h b/include/HRGEventGenerator/SimpleEvent.h index 350a0e4..298c083 100644 --- a/include/HRGEventGenerator/SimpleEvent.h +++ b/include/HRGEventGenerator/SimpleEvent.h @@ -60,11 +60,17 @@ namespace thermalfist { /// Output photons and leptons, if any bool printPhotonsLeptons; - /// Print the number of succesive decays before the particle was produced + /// Print the number of successive decays before the particle was produced bool printDecayEpoch; + /// Print the space-time coordinates of the particles + bool printCoordinates; + + /// Print the event weight for importance sampling + bool printWeight; + EventOutputConfig() : - printEnergy(false), printMotherPdg(false), printPhotonsLeptons(false), printDecayEpoch(false) { } + printEnergy(true), printMotherPdg(false), printPhotonsLeptons(false), printDecayEpoch(false), printCoordinates(false), printWeight(true) { } }; /// Writes the event to an output file stream @@ -72,6 +78,9 @@ namespace thermalfist { /// Writes the event to an output file stream void writeToFile(std::ofstream& fout, int eventnumber = 1) { writeToFile(fout, EventOutputConfig(), eventnumber); } + + /// Writes the event in a format suitable for UrQMD afterburner, as described here https://github.com/jbernhard/urqmd-afterburner + void writeToFileForUrqmd(std::ofstream& fout); }; } // namespace thermalfist diff --git a/include/HRGEventGenerator/SimpleParticle.h b/include/HRGEventGenerator/SimpleParticle.h index 0d802f5..0d55302 100644 --- a/include/HRGEventGenerator/SimpleParticle.h +++ b/include/HRGEventGenerator/SimpleParticle.h @@ -13,22 +13,28 @@ namespace thermalfist { /// Structure holding information about a single particle in the event generator. struct SimpleParticle { - double px, py, pz; ///< 3-momentum components (in GeV) - double m; ///< Mass (in GeV) - double p0; ///< Energy (in GeV) + double px, py, pz; ///< 3-momentum components (in GeV) + double m; ///< Mass (in GeV) + double p0; ///< Energy (in GeV) long long PDGID; ///< PDG code long long MotherPDGID; ///< PDG code of a mother particle, if applicable int epoch; ///< 0 - primary particle, 1 - after decay of primary particles, 2 - after a casacde of two decays and so on... bool processed; ///< Used in event generator + double r0, rx, ry, rz; ///< Space-time coordinates + /// Constructs a particle from provided three-momentum, mass, and PDG code - SimpleParticle(double inPx = 0., double inPy = 0., double inPz = 0., double inM = 0., long long inPDGID = 0., long long inMotherPDGID = 0) : + SimpleParticle(double inPx = 0., double inPy = 0., double inPz = 0., double inM = 0., long long inPDGID = 0., long long inMotherPDGID = 0, + double inR0 = 0., double inRx = 0., double inRy = 0., double inRz = 0.) : px(inPx), py(inPy), pz(inPz), m(inM), p0(sqrt(m*m + px * px + py * py + pz * pz)), - PDGID(inPDGID), MotherPDGID(inMotherPDGID), epoch(0), processed(false) { } + PDGID(inPDGID), MotherPDGID(inMotherPDGID), + epoch(0), processed(false), + r0(inR0), rx(inRx), ry(inRy), rz(inRz) + { } /// Absolute value of the 3-momentum (in GeV) double GetP() const { @@ -56,10 +62,27 @@ namespace thermalfist { return 0.5*log((GetP() + pz) / (GetP() - pz)); } + /// The longitudinal proper time + double GetTau() const { + return sqrt(r0 * r0 - rz * rz); + } + + /// The longitudinal space-time rapidity + double GetEtaS() const { + return 0.5 * log((r0 + rz) / (r0 - rz)); + } + /// Rapidity boost void RapidityBoost(double dY) { pz = GetMt() * sinh(GetY() + dY); p0 = sqrt(m*m + px * px + py * py + pz * pz); + + double tau = GetTau(); + double etas = GetEtaS(); + etas += dY; + + r0 = tau * cosh(etas); + rz = tau * sinh(etas); } }; diff --git a/src/examples/PCE/PCE-Saha-LHC.cpp b/src/examples/PCE/PCE-Saha-LHC.cpp index beab7f0..7c911b0 100644 --- a/src/examples/PCE/PCE-Saha-LHC.cpp +++ b/src/examples/PCE/PCE-Saha-LHC.cpp @@ -77,7 +77,7 @@ int main(int argc, char *argv[]) // The list of yield ratios to output vector< pair > ratios; // First the nuclei - ratios.push_back(make_pair(1000010020, 2212)); // d/p + ratios.push_back(make_pair(1000010020, 2212)); // d/p ratios.push_back(make_pair(1000020030, 2212)); // He3/p ratios.push_back(make_pair(1000010030, 2212)); // H3/p ratios.push_back(make_pair(1000020040, 2212)); // He4/p diff --git a/src/gui/QtThermalFIST/aboutdialog.cpp b/src/gui/QtThermalFIST/aboutdialog.cpp index f36e60b..c755e29 100644 --- a/src/gui/QtThermalFIST/aboutdialog.cpp +++ b/src/gui/QtThermalFIST/aboutdialog.cpp @@ -48,7 +48,7 @@ AboutDialog::AboutDialog(QWidget *parent) : labelVersion->setFont(fontDefault); layout->addWidget(labelVersion, 0, Qt::AlignCenter); layout->addSpacing(20); - QLabel *labelCC = new QLabel(tr("Copyright (c) 2014-2020 Volodymyr Vovchenko")); + QLabel *labelCC = new QLabel(tr("Copyright (c) 2014-2021 Volodymyr Vovchenko")); labelCC->setFont(fontDefault); layout->addWidget(labelCC, 0, Qt::AlignCenter); QLabel *labelLic = new QLabel(tr("GNU General Public License (GPLv3 or later)")); diff --git a/src/gui/QtThermalFIST/correlationsdialog.cpp b/src/gui/QtThermalFIST/correlationsdialog.cpp index ab347a2..dc5ee50 100644 --- a/src/gui/QtThermalFIST/correlationsdialog.cpp +++ b/src/gui/QtThermalFIST/correlationsdialog.cpp @@ -51,6 +51,7 @@ CorrelationsDialog::CorrelationsDialog(QWidget* parent, ThermalModelBase* mod) : QLabel* labelQuantity = new QLabel(tr("Quantity:")); comboQuantity = new QComboBox(); comboQuantity->addItem(tr("Susceptibility")); + comboQuantity->addItem(tr("Moment")); comboQuantity->addItem(tr("Delta[N1,N2]")); comboQuantity->addItem(tr("Sigma[N1,N2]")); comboQuantity->setCurrentIndex(0); @@ -146,6 +147,8 @@ void CorrelationsDialog::recalculate() if (comboQuantity->currentIndex() == 0) tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr))); + else if (comboQuantity->currentIndex() == 1) + tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3)))); else { if (pdg1 == pdg2) { tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(0.))); @@ -179,12 +182,12 @@ void CorrelationsDialog::recalculate() - model->TwoParticleSusceptibilityPrimordialByPdg(pdg1, -pdg1) - model->TwoParticleSusceptibilityPrimordialByPdg(-pdg1, pdg1) + model->TwoParticleSusceptibilityPrimordialByPdg(-pdg1, -pdg1)) / N1; - wn1 *= model->Volume() * pow(model->Parameters().T, 3); + wn1 *= model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3); wn2 = (model->TwoParticleSusceptibilityPrimordialByPdg(pdg2, pdg2) - model->TwoParticleSusceptibilityPrimordialByPdg(pdg2, -pdg2) - model->TwoParticleSusceptibilityPrimordialByPdg(-pdg2, pdg2) + model->TwoParticleSusceptibilityPrimordialByPdg(-pdg2, -pdg2)) / N2; - wn2 *= model->Volume() * pow(model->Parameters().T, 3); + wn2 *= model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3); } else { N1 = model->GetDensity(pdg1, Feeddown::StabilityFlag) * model->Volume() @@ -207,7 +210,7 @@ void CorrelationsDialog::recalculate() corr *= model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3); // Delta - if (comboQuantity->currentIndex() == 1) { + if (comboQuantity->currentIndex() == 2) { double DeltaN1N2 = (N1 * wn2 - N2 * wn1) / (N2 - N1); tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(DeltaN1N2))); continue; @@ -252,6 +255,9 @@ void CorrelationsDialog::recalculate() if (comboQuantity->currentIndex() == 0) { tableCorr->setItem(i1, i2, new QTableWidgetItem(QString::number(model->Susc((ConservedCharge::Name)i, (ConservedCharge::Name)j)))); } + else if (comboQuantity->currentIndex() == 1) { + tableCorr->setItem(i1, i2, new QTableWidgetItem(QString::number(model->Susc((ConservedCharge::Name)i, (ConservedCharge::Name)j) * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3)))); + } else { double N1 = model->ConservedChargeDensity((ConservedCharge::Name)i) * model->Volume(); double N2 = model->ConservedChargeDensity((ConservedCharge::Name)j) * model->Volume(); @@ -259,7 +265,7 @@ void CorrelationsDialog::recalculate() double wn2 = model->Susc((ConservedCharge::Name)j, (ConservedCharge::Name)j) * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3) / N2; double corr = model->Susc((ConservedCharge::Name)i, (ConservedCharge::Name)j) * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3); // Delta - if (comboQuantity->currentIndex() == 1) { + if (comboQuantity->currentIndex() == 2) { double DeltaN1N2 = (N1 * wn2 - N2 * wn1) / (N2 - N1); tableCorr->setItem(i1, i2, new QTableWidgetItem(QString::number(DeltaN1N2))); } @@ -333,6 +339,8 @@ void CorrelationsDialog::recalculate() if (comboQuantity->currentIndex() == 0) tableCorr->setItem(i, i2, new QTableWidgetItem(QString::number(corr))); + else if (comboQuantity->currentIndex() == 1) + tableCorr->setItem(i, i2, new QTableWidgetItem(QString::number(corr * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3)))); else { double N1 = 0., N2 = 0.; double wn1 = 0., wn2 = 0.; @@ -360,7 +368,7 @@ void CorrelationsDialog::recalculate() - model->TwoParticleSusceptibilityPrimordialByPdg(pdg1, -pdg1) - model->TwoParticleSusceptibilityPrimordialByPdg(-pdg1, pdg1) + model->TwoParticleSusceptibilityPrimordialByPdg(-pdg1, -pdg1)) / N1; - wn1 *= model->Volume() * pow(model->Parameters().T, 3); + wn1 *= model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3); wn2 = model->Susc((ConservedCharge::Name)cc, (ConservedCharge::Name)cc) * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3) / N2; } else { @@ -379,7 +387,7 @@ void CorrelationsDialog::recalculate() corr *= model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3); // Delta - if (comboQuantity->currentIndex() == 1) { + if (comboQuantity->currentIndex() == 2) { double DeltaN1N2 = (N1 * wn2 - N2 * wn1) / (N2 - N1); tableCorr->setItem(i, i2, new QTableWidgetItem(QString::number(DeltaN1N2))); } diff --git a/src/gui/QtThermalFIST/equationofstatetab.cpp b/src/gui/QtThermalFIST/equationofstatetab.cpp index f55d173..a376705 100644 --- a/src/gui/QtThermalFIST/equationofstatetab.cpp +++ b/src/gui/QtThermalFIST/equationofstatetab.cpp @@ -641,9 +641,23 @@ void EquationOfStateTab::plotLatticeData() } // WB data - if (mapWB.count(paramname) > 0) { + if (mapWB.count(paramname) > 0 + || (CBratio->isChecked() + && mapWB.count(paramnames[index]) && mapWB.count(paramnames[comboQuantity2->currentIndex()]) + && dataWBx[mapWB[paramnames[index]]].size() == dataWBx[mapWB[paramnames[comboQuantity2->currentIndex()]]].size()) + ) { - int indWB = mapWB[paramname]; + int indWB = 0; + if (mapWB.count(paramname)) + indWB = mapWB[paramname]; + int indWB2 = 0; + if (CBratio->isChecked()) { + int index2 = comboQuantity2->currentIndex(); + indWB = mapWB[paramnames[index]]; + indWB2 = mapWB[paramnames[index2]]; + if (dataWBx[indWB].size() != dataWBx[indWB2].size()) + return; + } plotDependence->addGraph(); plotDependence->graph(graphStart)->setName("LQCD (Wuppertal-Budapest)"); plotDependence->graph(graphStart)->setPen(QPen(QColor(255, 255, 255, 0))); @@ -659,8 +673,18 @@ void EquationOfStateTab::plotLatticeData() QVector yConfUpper(dataWBx[indWB].size()), yConfLower(dataWBx[indWB].size()); for (int i = 0; i < x0.size(); ++i) { x0[i] = dataWBx[indWB][i]; - yConfUpper[i] = dataWBy[indWB][i] + dataWByerrp[indWB][i]; - yConfLower[i] = dataWBy[indWB][i] - dataWByerrm[indWB][i]; + if (mapWB.count(paramname)) { + yConfUpper[i] = dataWBy[indWB][i] + dataWByerrp[indWB][i]; + yConfLower[i] = dataWBy[indWB][i] - dataWByerrm[indWB][i]; + } + else if (CBratio->isChecked()) { + double mean = dataWBy[indWB][i] / dataWBy[indWB2][i]; + double error = mean * sqrt((dataWByerrp[indWB][i]/dataWBy[indWB][i])*(dataWByerrp[indWB][i]/dataWBy[indWB][i]) + + (dataWByerrp[indWB2][i]/dataWBy[indWB2][i])*(dataWByerrp[indWB2][i]/dataWBy[indWB2][i])); + + yConfUpper[i] = mean + error; + yConfLower[i] = mean - error; + } if (x0[i] >= spinTMin->value() && x0[i] <= spinTMax->value()) { tmin = std::min(tmin, yConfLower[i]); tmax = std::max(tmax, yConfUpper[i]); @@ -676,8 +700,23 @@ void EquationOfStateTab::plotLatticeData() } // HotQCD data - if (mapHotQCD.count(paramname) > 0) { - int indHotQCD = mapHotQCD[paramname]; + if (mapHotQCD.count(paramname) > 0 || (CBratio->isChecked() + && mapHotQCD.count(paramnames[index]) && mapHotQCD.count(paramnames[comboQuantity2->currentIndex()]) + && dataHotQCDx[mapHotQCD[paramnames[index]]].size() == dataHotQCDx[mapWB[paramnames[comboQuantity2->currentIndex()]]].size()) + ) + { + int indHotQCD = 0; + if (mapHotQCD.count(paramname)) + indHotQCD = mapHotQCD[paramname]; + int indHotQCD2 = 0; + if (CBratio->isChecked()) { + int index2 = comboQuantity2->currentIndex(); + indHotQCD = mapHotQCD[paramnames[index]]; + indHotQCD2 = mapHotQCD[paramnames[index2]]; + if (dataHotQCDx[indHotQCD].size() != dataHotQCDx[indHotQCD2].size()) + return; + } + plotDependence->addGraph(); plotDependence->graph(graphStart)->setName("LQCD (HotQCD)"); plotDependence->graph(graphStart)->setPen(QPen(QColor(255, 255, 255, 0))); @@ -693,8 +732,20 @@ void EquationOfStateTab::plotLatticeData() QVector yConfUpper(dataHotQCDx[indHotQCD].size()), yConfLower(dataHotQCDx[indHotQCD].size()); for (int i = 0; i < x0.size(); ++i) { x0[i] = dataHotQCDx[indHotQCD][i]; - yConfUpper[i] = dataHotQCDy[indHotQCD][i] + dataHotQCDyerrp[indHotQCD][i]; - yConfLower[i] = dataHotQCDy[indHotQCD][i] - dataHotQCDyerrm[indHotQCD][i]; + + if (mapHotQCD.count(paramname)) { + yConfUpper[i] = dataHotQCDy[indHotQCD][i] + dataHotQCDyerrp[indHotQCD][i]; + yConfLower[i] = dataHotQCDy[indHotQCD][i] - dataHotQCDyerrm[indHotQCD][i]; + } + else if (CBratio->isChecked()) { + double mean = dataHotQCDy[indHotQCD][i] / dataHotQCDy[indHotQCD2][i]; + double error = mean * sqrt((dataHotQCDyerrp[indHotQCD][i]/dataHotQCDy[indHotQCD][i])*(dataHotQCDyerrp[indHotQCD][i]/dataHotQCDy[indHotQCD][i]) + + (dataHotQCDyerrp[indHotQCD2][i]/dataHotQCDy[indHotQCD2][i])*(dataHotQCDyerrp[indHotQCD2][i]/dataHotQCDy[indHotQCD2][i])); + + yConfUpper[i] = mean + error; + yConfLower[i] = mean - error; + } + if (x0[i] >= spinTMin->value() && x0[i] <= spinTMax->value()) { tmin = std::min(tmin, yConfLower[i]); tmax = std::max(tmax, yConfUpper[i]); diff --git a/src/gui/QtThermalFIST/modeltab.cpp b/src/gui/QtThermalFIST/modeltab.cpp index 2b29894..fc135ae 100644 --- a/src/gui/QtThermalFIST/modeltab.cpp +++ b/src/gui/QtThermalFIST/modeltab.cpp @@ -521,6 +521,7 @@ void ModelTab::performCalculation(const ThermalModelConfig & config) dbgstrm << "Net strangeness\t= " << model->CalculateStrangenessDensity() * model->Volume() << endl; if (model->TPS()->hasCharmed()) dbgstrm << "Net charm\t= " << model->CalculateCharmDensity() * model->Volume() << endl; + dbgstrm << "Absolute baryon number\t= " << model->AbsoluteBaryonDensity() * model->Volume() << endl; dbgstrm << "E/N\t\t= " << model->CalculateEnergyDensity() / model->CalculateHadronDensity() << endl; if (fabs(nb) > 1.e-10) dbgstrm << "E/Nb\t\t= " << model->CalculateEnergyDensity() / nb << endl; diff --git a/src/gui/QtThermalFIST/particledialog.cpp b/src/gui/QtThermalFIST/particledialog.cpp index 6cc8d70..51c6640 100644 --- a/src/gui/QtThermalFIST/particledialog.cpp +++ b/src/gui/QtThermalFIST/particledialog.cpp @@ -106,7 +106,8 @@ ParticleDialog::ParticleDialog(QWidget* parent, ThermalModelBase* mod, int Parti if (model->IsCalculated()) { std::vector< std::pair > sources = model->TPS()->DecayContributionsByFeeddown()[Feeddown::StabilityFlag][pid]; - for (int i = 0; i < sources.size(); ++i) sources[i].first *= model->Densities()[sources[i].second]; + for (int i = 0; i < sources.size(); ++i) + sources[i].first *= model->Densities()[sources[i].second]; qSort(sources.begin(), sources.end()); diff --git a/src/gui/QtThermalFIST/particlespectra.cpp b/src/gui/QtThermalFIST/particlespectra.cpp index 61b483b..2cf4128 100644 --- a/src/gui/QtThermalFIST/particlespectra.cpp +++ b/src/gui/QtThermalFIST/particlespectra.cpp @@ -24,6 +24,11 @@ void ParticleSpectrum::AddParticle(const SimpleParticle &part) { dndmt.insert(part.GetMt()); dndpt.insert(part.GetPt()); d2ndptdy.insert(part.GetY(), part.GetPt()); + + double pT = part.GetPt(); + pTsum_event += pT; + pT2sum_event += pT * pT; + //pTcnt += 1.; } void ParticleSpectrum::FinishEvent(double weight) { @@ -38,6 +43,13 @@ void ParticleSpectrum::FinishEvent(double weight) { n8 += weight*tmn*tmn*tmn*tmn*tmn*tmn*tmn*tmn; wsum += weight; w2sum += weight*weight; + + pTsum += pTsum_event * weight; + pT2sum += pT2sum_event * weight; + pTcnt += tmn * weight; + pTcnt2 += tmn * weight * weight; + pTsum_event = pT2sum_event = 0.; + tmpn = 0; dndp.updateEvent(weight); dndy.updateEvent(weight); @@ -142,6 +154,19 @@ double ParticleSpectrum::GetScaledVarianceError() const { return GetVarianceError() / GetMean(); } +double ParticleSpectrum::GetMeanPt() const +{ + return pTsum / pTcnt; +} + +double ParticleSpectrum::GetMeanPtError() const +{ + double pTmean = pTsum / pTcnt; + double pT2mean = pT2sum / pTcnt; + double nEpT = pTcnt * pTcnt / pTcnt2; + return sqrt((pT2mean - pTmean * pTmean) / (nEpT - 1.)); +} + ParticlesSpectra::ParticlesSpectra(ThermalModelBase *model, double T, double beta, int distrtype, double etamax) { fEtaMax = etamax; fDistributionType = distrtype; diff --git a/src/gui/QtThermalFIST/particlespectra.h b/src/gui/QtThermalFIST/particlespectra.h index 8d12c8a..6ed95c2 100644 --- a/src/gui/QtThermalFIST/particlespectra.h +++ b/src/gui/QtThermalFIST/particlespectra.h @@ -231,6 +231,9 @@ class ParticleSpectrum std::vector xs; std::vector ys; bool acc; + double pTsum_event, pT2sum_event; + double pTsum, pT2sum; + double pTcnt, pTcnt2; public: density dndp; density dndy; @@ -246,7 +249,10 @@ class ParticleSpectrum d2ndptdy = density2d(-3. - etamax, 3. + etamax, 40, 0., 2., 40); acc = false; isAveragesCalculated = false; + tmpn = 0; means.resize(9); + pTsum = pT2sum = pTcnt = pTcnt2 = 0.; + pTsum_event = pT2sum_event = 0.; } ~ParticleSpectrum() { } @@ -262,6 +268,8 @@ class ParticleSpectrum dndmt = density(mass, 2.+mass, 500); dndpt = density(0., 3. + mass, 500); d2ndptdy = density2d(-3., 3., 40, 0., 2., 40); + pTsum = pT2sum = pTcnt = pTcnt2 = 0.; + pTsum_event = pT2sum_event = 0.; } void SetDistribution(thermalfist::MomentumDistributionBase *distr) { fDistribution = distr; @@ -291,6 +299,10 @@ class ParticleSpectrum double GetVarianceError() const; double GetScaledVarianceError() const; + double GetMeanPt() const; + double GetMeanPtError() const; + + double GetModeldNdp(double p) const { if (fDistribution!=NULL && !fDistribution->isNormalized()) fDistribution->Normalize(); if (fDistribution==NULL) return 0.; diff --git a/src/gui/QtThermalFIST/resultdialog.cpp b/src/gui/QtThermalFIST/resultdialog.cpp index 3c70013..b9b2bd7 100644 --- a/src/gui/QtThermalFIST/resultdialog.cpp +++ b/src/gui/QtThermalFIST/resultdialog.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include "HRGBase/ThermalModelBase.h" #include "HRGBase/xMath.h" @@ -374,6 +375,33 @@ QString ResultDialog::GetResults() { ret += "\r\n"; + { + double Np = model->ChargedMultiplicity(1); + double Nm = model->ChargedMultiplicity(-1); + double wplus = model->ChargedScaledVariance(1); + double wminus = model->ChargedScaledVariance(-1); + double covNpm = + (Np * wplus + Nm * wminus - (model->Susc(ConservedCharge::ElectricCharge, ConservedCharge::ElectricCharge) + * model->Volume() * pow(model->Parameters().T * xMath::GeVtoifm(), 3))) / 2.; + + double nudyn = wplus / Np + wminus / Nm - 2. * covNpm / Np / Nm - (Np + Nm) / Np / Nm; + + sprintf(cc, "%-25s = ", "Primordial *nu_dyn"); + ret += QString(cc); + ret += QString::number((Np + Nm) * nudyn) + "\r\n"; + + double Dprim = 4. * model->Susc(ConservedCharge::ElectricCharge, ConservedCharge::ElectricCharge) + * model->Volume() * pow(model->Parameters().T * xMath::GeVtoifm(), 3) / + model->ChargedMultiplicity(0); + + sprintf(cc, "%-25s = ", "Primordial D"); + ret += QString(cc); + ret += QString::number(Dprim) + "\r\n"; + + ret += "\r\n"; + } + + sprintf(cc, "%-25s = ", "Final Nch"); ret += QString(cc); ret += QString::number(model->ChargedMultiplicityFinal(0)) + "\r\n"; @@ -399,6 +427,34 @@ QString ResultDialog::GetResults() { ret += QString(cc); ret += QString::number(model->ChargedScaledVarianceFinal(-1)) + "\r\n"; + ret += "\r\n"; + + { + double Np = model->ChargedMultiplicityFinal(1); + double Nm = model->ChargedMultiplicityFinal(-1); + double wplus = model->ChargedScaledVarianceFinal(1); + double wminus = model->ChargedScaledVarianceFinal(-1); + double covNpm = + (Np * wplus + Nm * wminus - (model->Susc(ConservedCharge::ElectricCharge, ConservedCharge::ElectricCharge) + * model->Volume() * pow(model->Parameters().T * xMath::GeVtoifm(), 3))) / 2.; + + double nudyn = wplus / Np + wminus / Nm - 2. * covNpm / Np / Nm - (Np + Nm) / Np / Nm; + + sprintf(cc, "%-25s = ", "Final *nu_dyn"); + ret += QString(cc); + ret += QString::number((Np + Nm) * nudyn) + "\r\n"; + + double Dtot = 4. * model->Susc(ConservedCharge::ElectricCharge, ConservedCharge::ElectricCharge) + * model->Volume() * pow(model->Parameters().T * xMath::GeVtoifm(), 3) / + model->ChargedMultiplicityFinal(0); + + sprintf(cc, "%-25s = ", "Final D"); + ret += QString(cc); + ret += QString::number(Dtot) + "\r\n"; + + ret += "\r\n"; + } + if (flucts != NULL && flucts->flag) { ret += "\r\n"; diff --git a/src/gui/QtThermalFIST/spectramodel.cpp b/src/gui/QtThermalFIST/spectramodel.cpp index 3c07e2f..b3637ef 100644 --- a/src/gui/QtThermalFIST/spectramodel.cpp +++ b/src/gui/QtThermalFIST/spectramodel.cpp @@ -123,7 +123,9 @@ QVariant SpectraModel::data(const QModelIndex &index, int role) const if (col == 7) return QString::number(spectra->fNegativeCharges[tind].GetKurtosis(), 'f', 3) + " ± " + QString::number(spectra->fNegativeCharges[tind].GetKurtosisError(), 'f', 3); } } - if (col == 8 && row < spectra->fParticles.size() && spectra->fParticles[RowToParticle[row]].GetAcceptance()) return "*"; + if (col == 8 && row < spectra->fParticles.size()) + return QString::number(spectra->fParticles[RowToParticle[row]].GetMeanPt(), 'f', 3) + " ± " + QString::number(spectra->fParticles[RowToParticle[row]].GetMeanPtError(), 'f', 3); + //if (col == 8 && row < spectra->fParticles.size() && spectra->fParticles[RowToParticle[row]].GetAcceptance()) return "*"; break; case Qt::FontRole: // if (model->TPS()->Particles()[RowToParticle[row]].IsStable()) //change font only for cell(0,0) @@ -172,7 +174,7 @@ QVariant SpectraModel::headerData(int section, Qt::Orientation orientation, int case 0: return tr("Name"); case 2: - return tr("Mass"); + return tr("m [GeV]"); case 3: return tr("Multiplicity"); case 4: @@ -184,7 +186,7 @@ QVariant SpectraModel::headerData(int section, Qt::Orientation orientation, int case 7: return tr("Kurtosis"); case 8: - return tr("Acceptance"); + return tr(" [GeV]"); } } else if (orientation == Qt::Vertical) return section + 1; diff --git a/src/library/HRGBase/ThermalModelBase.cpp b/src/library/HRGBase/ThermalModelBase.cpp index 1747e28..5a62073 100644 --- a/src/library/HRGBase/ThermalModelBase.cpp +++ b/src/library/HRGBase/ThermalModelBase.cpp @@ -743,7 +743,21 @@ namespace thermalfist { if (feeddown != Feeddown::Primordial && !m_FeeddownCalculated) CalculateFeeddown(); - return GetDensity(PDGID, dens); + double ret = GetDensity(PDGID, dens); + + // Weak decay contributions from K0S if this particle is not in the list + if (feeddown == Feeddown::Weak && m_TPS->PdgToId(310) == -1) { + // pi0 + if (PDGID == 111) { + ret += 2. * 0.308 * GetDensity(310, dens); + } + // pi+,- + if (PDGID == 211 || PDGID == -211) { + ret += 0.692 * GetDensity(310, dens); + } + } + + return ret; } @@ -943,7 +957,6 @@ namespace thermalfist { m_TotalCorrel[i][j] += m_densities[r] / m_Parameters.T * dnij; } } - } } } diff --git a/src/library/HRGBase/ThermalParticle.cpp b/src/library/HRGBase/ThermalParticle.cpp index 7b0baea..b453c70 100644 --- a/src/library/HRGBase/ThermalParticle.cpp +++ b/src/library/HRGBase/ThermalParticle.cpp @@ -26,7 +26,8 @@ namespace thermalfist { ThermalParticle::ThermalParticle(bool Stable, std::string Name, long long PDGID, double Deg, int Stat, double Mass, int Strange, int Baryon, int Charge, double AbsS, double Width, double Threshold, int Charm, double AbsC, int Quark) : m_Stable(Stable), m_AntiParticle(false), m_Name(Name), m_PDGID(PDGID), m_Degeneracy(Deg), m_Statistics(Stat), m_StatisticsOrig(Stat), m_Mass(Mass), - m_Strangeness(Strange), m_Baryon(Baryon), m_ElectricCharge(Charge), m_Charm(Charm), m_ArbitraryCharge(Baryon), m_AbsS(AbsS), m_AbsC(AbsC), m_Width(Width), m_Threshold(Threshold), m_Quark(Quark), m_Weight(1.) + m_Strangeness(Strange), m_Baryon(Baryon), m_ElectricCharge(Charge), m_Charm(Charm), m_ArbitraryCharge(Baryon), m_AbsS(AbsS), m_AbsC(AbsC), m_Width(Width), m_Threshold(Threshold), m_Quark(Quark), m_Weight(1.), + m_ResonanceWidthShape(ThermalParticle::RelativisticBreitWigner), m_ResonanceWidthIntegrationType(ThermalParticle::ZeroWidth) { if (!Disclaimer::DisclaimerPrinted) Disclaimer::DisclaimerPrinted = Disclaimer::PrintDisclaimer(); diff --git a/src/library/HRGBase/ThermalParticleSystem.cpp b/src/library/HRGBase/ThermalParticleSystem.cpp index ccdec01..0455a44 100644 --- a/src/library/HRGBase/ThermalParticleSystem.cpp +++ b/src/library/HRGBase/ThermalParticleSystem.cpp @@ -486,6 +486,11 @@ namespace thermalfist { m_NumBaryons = m_NumCharged = m_NumStrange = m_NumCharmed = 0; + if (ListFiles.size() == 1 && CheckListIsiSS(ListFiles[0])) { + LoadListiSS(ListFiles[0], flags, mcut); + return; + } + for(int i = 0; i < ListFiles.size(); ++i) AddParticlesToListFromFile(ListFiles[i], flags, mcut); @@ -529,11 +534,7 @@ namespace thermalfist { LoadDecays(tDecayFiles, flags); - SetResonanceWidthShape(m_ResonanceWidthShape); - SetResonanceWidthIntegrationType(m_ResonanceWidthIntegrationType); - SetCalculationType(m_QStatsCalculationType); - - CheckDecayChannelsAreSpecified(); + FinalizeListLoad(); } void ThermalParticleSystem::LoadList(const std::string& InputFile, const std::string& DecayFile, bool GenAntiP, double mcut) { @@ -846,12 +847,7 @@ namespace thermalfist { } } - for (size_t i = 0; i < m_Particles.size(); ++i) - m_Particles[i].SetDecaysOriginal(m_Particles[i].Decays()); - - FillDecayProperties(); - FillDecayThresholds(); - ProcessDecays(); + FinalizeDecaysLoad(); } void ThermalParticleSystem::LoadDecays(const std::string& DecaysFile, bool GenerateAntiParticles) @@ -950,6 +946,8 @@ namespace thermalfist { m_NumberOfParticles = 0; m_Particles.resize(0); m_PDGtoID.clear(); + m_ResonanceWidthShape = ThermalParticle::RelativisticBreitWigner; + m_ResonanceWidthIntegrationType = ThermalParticle::ZeroWidth; m_SortMode = ThermalParticleSystem::SortByMass; @@ -970,6 +968,171 @@ namespace thermalfist { Initialize(vector(1, InputFile), vector(1, DecayFile), flags, mcut); } + void ThermalParticleSystem::FinalizeListLoad() + { + SetResonanceWidthShape(m_ResonanceWidthShape); + SetResonanceWidthIntegrationType(m_ResonanceWidthIntegrationType); + SetCalculationType(m_QStatsCalculationType); + + CheckDecayChannelsAreSpecified(); + } + + void ThermalParticleSystem::FinalizeDecaysLoad() + { + for (size_t i = 0; i < m_Particles.size(); ++i) + m_Particles[i].SetDecaysOriginal(m_Particles[i].Decays()); + + FillDecayProperties(); + FillDecayThresholds(); + ProcessDecays(); + } + + bool ThermalParticleSystem::CheckListIsiSS(const std::string& filename) + { + ifstream fin(filename.c_str()); + if (fin.is_open()) { + string tmp; + fin >> tmp; + return (tmp == "iSS"); + } + return false; + } + + void ThermalParticleSystem::LoadListiSS(const std::string& filename, const std::set& flags, double mcut) + { + m_NumberOfParticles = 0; + m_Particles.resize(0); + m_PDGtoID.clear(); + + m_NumBaryons = m_NumCharged = m_NumStrange = m_NumCharmed = 0; + + ifstream fin(filename.c_str()); + + if (fin.is_open()) { + char cc[2000]; + // Skip the header line + fin.getline(cc, 2000); + // Read all the particles + while (!fin.eof()) { + fin.getline(cc, 2000); + string tmp = string(cc); + vector elems = CuteHRGHelper::split(tmp, '#'); + if (elems.size() < 1 || elems[0].size() == 0) + continue; + + istringstream iss(elems[0]); + + int stable, stat, str, bary, chg, charm, itmp, tdecaynumber; + long long pdgid, ltmp; + double mass, degeneracy, width, threshold, abss, absc; + string name; + + if (iss >> pdgid + >> name + >> mass + >> width + >> degeneracy + >> bary + >> str + >> charm + >> itmp + >> itmp + >> chg + >> tdecaynumber + ) { + + stable = 0; + abss = std::abs(str); + if (pdgid == 333) + abss = 2.; + absc = std::abs(charm); + threshold = 0.; + stat = (std::abs(bary) % 2 == 0 ? -1 : 1); + + ThermalParticle part_candidate = ThermalParticle( + (bool)stable, name, pdgid, degeneracy, stat, mass, str, bary, chg, abss, width, threshold, charm, absc); + + if (pdgid < 0) + part_candidate.SetAntiParticle(true); + + if (!AcceptParticle(part_candidate, flags, mcut) || m_PDGtoID.count(pdgid) != 0) + continue; + + if (bary != 0) m_NumBaryons++; + if (chg != 0) m_NumCharged++; + if (str != 0) m_NumStrange++; + if (charm != 0) m_NumCharmed++; + + m_Particles.push_back(part_candidate); + m_PDGtoID[pdgid] = m_Particles.size() - 1; + m_NumberOfParticles++; + + // now read the decays + ThermalParticle::ParticleDecaysVector tdecays(0); + for (size_t idec = 0; idec < tdecaynumber; ++idec) { + fin.getline(cc, 2000); + string tmp = string(cc); + elems = CuteHRGHelper::split(tmp, '#'); + istringstream issdec(elems[0]); + + ParticleDecayChannel tdecay; + issdec >> ltmp; + int ndecayparticles; + double BR; + issdec >> ndecayparticles >> tdecay.mBratio; + for (size_t idaughter = 0; idaughter < 5; ++idaughter) { + issdec >> ltmp; + if (idaughter < ndecayparticles) + tdecay.mDaughters.push_back(ltmp); + } + tdecays.push_back(tdecay); + } + + if (static_cast(tdecays.size()) == tdecaynumber + && static_cast(tdecays.size()) != 0 + && tdecays[0].mDaughters.size() != 1) { + m_Particles.back().SetDecays(tdecays); + } + else { + stable = 1; + m_Particles.back().SetStable(true); + } + + // Add antibaryon + if (bary > 0 && (m_PDGtoID.count(-pdgid) == 0)) { + + if (bary != 0) m_NumBaryons++; + if (chg != 0) m_NumCharged++; + if (str != 0) m_NumStrange++; + if (charm != 0) m_NumCharmed++; + + if (bary == 0 && name[name.size() - 1] == '+') + name[name.size() - 1] = '-'; + else if (bary == 0 && name[name.size() - 1] == '-') + name[name.size() - 1] = '+'; + else + name = "anti-" + name; + m_Particles.push_back(ThermalParticle((bool)stable, name, -pdgid, degeneracy, stat, mass, -str, -bary, -chg, abss, width, threshold, -charm, absc)); + m_Particles.back().SetAntiParticle(true); + m_PDGtoID[pdgid] = m_Particles.size() - 1; + } + + } + } + + FinalizeList(); + + for (size_t i = 0; i < m_Particles.size(); ++i) { + if (m_Particles[i].BaryonCharge() < 0) + m_Particles[i].SetDecays(GetDecaysFromAntiParticle(m_Particles[m_PDGtoID[-m_Particles[i].PdgId()]].Decays())); + } + + FinalizeDecaysLoad(); + + FinalizeListLoad(); + } + } + void ThermalParticleSystem::WriteDecaysToFile(const std::string& OutputFile, bool WriteAntiParticles) { std::ofstream fout(OutputFile.c_str()); diff --git a/src/library/HRGBase/Utility.cpp b/src/library/HRGBase/Utility.cpp index 7d640bd..dd9cec3 100644 --- a/src/library/HRGBase/Utility.cpp +++ b/src/library/HRGBase/Utility.cpp @@ -88,7 +88,7 @@ namespace thermalfist { email += "."; email += "uni-frankfurt.de"; - tmpstr = "Copyright (c) 2020 Volodymyr Vovchenko <" + email + ">"; + tmpstr = "Copyright (c) 2021 Volodymyr Vovchenko <" + email + ">"; tmpstr = OutputString(tmpstr); diff --git a/src/library/HRGEV/ThermalModelEVCrossterms.cpp b/src/library/HRGEV/ThermalModelEVCrossterms.cpp index 732bea9..d8280f4 100644 --- a/src/library/HRGEV/ThermalModelEVCrossterms.cpp +++ b/src/library/HRGEV/ThermalModelEVCrossterms.cpp @@ -240,6 +240,37 @@ namespace thermalfist { void ThermalModelEVCrossterms::CalculatePrimordialDensities() { m_FluctuationsCalculated = false; + map< vector, int> m_MapEVcomponent; + + { + int NN = m_densities.size(); + m_MapToEVComponent.resize(NN); + m_MapFromEVComponent.clear(); + m_MapEVcomponent.clear(); + m_EVComponentIndices.clear(); + + int tind = 0; + for (int i = 0; i < NN; ++i) { + vector EVParam(0); + for (int j = 0; j < NN; ++j) { + EVParam.push_back(m_Virial[i][j]); + EVParam.push_back(m_Virial[j][i]); + } + + if (m_MapEVcomponent.count(EVParam) == 0) { + m_MapEVcomponent[EVParam] = tind; + m_MapToEVComponent[i] = tind; + m_MapFromEVComponent.push_back(i); + m_EVComponentIndices.push_back(vector(1, i)); + tind++; + } + else { + m_MapToEVComponent[i] = m_MapEVcomponent[EVParam]; + m_EVComponentIndices[m_MapEVcomponent[EVParam]].push_back(i); + } + } + } + // Pressure SolvePressure(); vector tN(m_densities.size()); diff --git a/src/library/HRGEventGenerator/CracowFreezeoutEventGenerator.cpp b/src/library/HRGEventGenerator/CracowFreezeoutEventGenerator.cpp index 5960f9c..80404b7 100644 --- a/src/library/HRGEventGenerator/CracowFreezeoutEventGenerator.cpp +++ b/src/library/HRGEventGenerator/CracowFreezeoutEventGenerator.cpp @@ -42,11 +42,15 @@ namespace thermalfist { { ClearMomentumGenerators(); m_BWGens.resize(0); + + // Find \tau_H from Veff / (\delta \eta) = \pi \tau_H R^2 where R = m_RoverTauH * \tau_H + double tauH = pow(m_THM->Volume() / (2. * m_EtaMax) / xMath::Pi() / m_RoverTauH / m_RoverTauH, 1. / 3.); + if (m_THM != NULL) { for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) { const ThermalParticle& part = m_THM->TPS()->Particles()[i]; //m_MomentumGens.push_back(new RandomGenerators::CracowFreezeoutMomentumGenerator(m_T, m_RoverTauH, m_EtaMax, part.Mass(), part.Statistics(), m_THM->FullIdealChemicalPotential(i))); - m_MomentumGens.push_back(new RandomGenerators::BoostInvariantMomentumGenerator(new CracowFreezeoutParametrization(m_RoverTauH), GetTkin(), GetEtaMax(), part.Mass(), part.Statistics(), m_THM->FullIdealChemicalPotential(i))); + m_MomentumGens.push_back(new RandomGenerators::BoostInvariantMomentumGenerator(new CracowFreezeoutParametrization(m_RoverTauH, tauH), GetTkin(), GetEtaMax(), part.Mass(), part.Statistics(), m_THM->FullIdealChemicalPotential(i))); double T = m_THM->Parameters().T; double Mu = m_THM->FullIdealChemicalPotential(i); diff --git a/src/library/HRGEventGenerator/CylindricalBlastWaveEventGenerator.cpp b/src/library/HRGEventGenerator/CylindricalBlastWaveEventGenerator.cpp index 3990b5f..f3472f0 100644 --- a/src/library/HRGEventGenerator/CylindricalBlastWaveEventGenerator.cpp +++ b/src/library/HRGEventGenerator/CylindricalBlastWaveEventGenerator.cpp @@ -10,6 +10,7 @@ #include #include "HRGBase/xMath.h" +#include "HRGBase/NumericalIntegration.h" #include "HRGBase/ThermalModelBase.h" #include "HRGEV/ExcludedVolumeModel.h" #include "HRGEventGenerator/ParticleDecaysMC.h" @@ -20,15 +21,15 @@ namespace thermalfist { // m_THM = NULL; //} - CylindricalBlastWaveEventGenerator::CylindricalBlastWaveEventGenerator(ThermalParticleSystem * TPS, const EventGeneratorConfiguration & config, double T, double betas, double etamax, double npow) : - m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow) + CylindricalBlastWaveEventGenerator::CylindricalBlastWaveEventGenerator(ThermalParticleSystem * TPS, const EventGeneratorConfiguration & config, double T, double betas, double etamax, double npow, double Rperp) : + m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow), m_Rperp(Rperp) { SetConfiguration(TPS, config); SetMomentumGenerators(); } - CylindricalBlastWaveEventGenerator::CylindricalBlastWaveEventGenerator(ThermalModelBase *THM, double T, double betas, double etamax, double npow, bool /*onlyStable*/, EventGeneratorConfiguration::ModelType EV, ThermalModelBase *THMEVVDW) :m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow) { + CylindricalBlastWaveEventGenerator::CylindricalBlastWaveEventGenerator(ThermalModelBase *THM, double T, double betas, double etamax, double npow, bool /*onlyStable*/, EventGeneratorConfiguration::ModelType EV, ThermalModelBase *THMEVVDW) :m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow), m_Rperp(6.5) { EventGeneratorConfiguration::ModelType modeltype = EV; EventGeneratorConfiguration::Ensemble ensemble = EventGeneratorConfiguration::GCE; if (THM->Ensemble() == ThermalModelBase::CE) @@ -89,10 +90,14 @@ namespace thermalfist { { ClearMomentumGenerators(); m_BWGens.resize(0); + + // Find \tau_H from Veff / (\delta \eta) = \pi \tau R^2 * I where I is an integral over transverse velocity profile computed numerically + double tau = m_THM->Volume() / (2. * GetEtaMax()) / (2. * xMath::Pi()) / GetRperp() / GetRperp() / GetVeffIntegral(); + if (m_THM != NULL) { for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) { const ThermalParticle& part = m_THM->TPS()->Particles()[i]; - m_MomentumGens.push_back(new RandomGenerators::BoostInvariantMomentumGenerator(new CylindricalBlastWaveParametrization(GetBetaSurface(), GetNPow()), GetTkin(), GetEtaMax(), part.Mass(), part.Statistics(), m_THM->FullIdealChemicalPotential(i))); + m_MomentumGens.push_back(new RandomGenerators::BoostInvariantMomentumGenerator(new CylindricalBlastWaveParametrization(GetBetaSurface(), GetNPow(), tau, GetRperp()), GetTkin(), GetEtaMax(), part.Mass(), part.Statistics(), m_THM->FullIdealChemicalPotential(i))); double T = m_THM->Parameters().T; double Mu = m_THM->FullIdealChemicalPotential(i); @@ -104,4 +109,21 @@ namespace thermalfist { } } + double CylindricalBlastWaveEventGenerator::GetVeffIntegral() const + { + double ret = 0.0; + + std::vector xleg, wleg; + NumericalIntegration::GetCoefsIntegrateLegendre32(0., 1., &xleg, &wleg); + + for (int iint = 0; iint < xleg.size(); ++iint) { + double zeta = xleg[iint]; + double w = wleg[iint]; + double betar = pow(zeta, GetNPow()) * GetBetaSurface(); + ret += w * zeta / sqrt(1. - betar * betar); + } + + return ret; + } + } // namespace thermalfist \ No newline at end of file diff --git a/src/library/HRGEventGenerator/EventGeneratorBase.cpp b/src/library/HRGEventGenerator/EventGeneratorBase.cpp index b32ab85..8434206 100644 --- a/src/library/HRGEventGenerator/EventGeneratorBase.cpp +++ b/src/library/HRGEventGenerator/EventGeneratorBase.cpp @@ -26,6 +26,29 @@ namespace thermalfist { int EventGeneratorBase::fCEAccepted, EventGeneratorBase::fCETotal; double EventGeneratorBase::m_LastWeight = 1., EventGeneratorBase::m_LastLogWeight = 0., EventGeneratorBase::m_LastNormWeight = 1.; + std::vector LorentzBoost(const std::vector& fourvector, double vx, double vy, double vz) + { + std::vector ret(4, 0); + double v2 = vx * vx + vy * vy + vz * vz; + if (v2 == 0.0) + return fourvector; + double gamma = 1. / sqrt(1. - v2); + + const double& r0 = fourvector[0]; + const double& rx = fourvector[1]; + const double& ry = fourvector[2]; + const double& rz = fourvector[3]; + + ret[0] = gamma * r0 - gamma * (vx * rx + vy * ry + vz * rz); + ret[1] = -gamma * vx * r0 + (1. + (gamma - 1.) * vx * vx / v2) * rx + + (gamma - 1.) * vx * vy / v2 * ry + (gamma - 1.) * vx * vz / v2 * rz; + ret[2] = -gamma * vy * r0 + (1. + (gamma - 1.) * vy * vy / v2) * ry + + (gamma - 1.) * vy * vx / v2 * rx + (gamma - 1.) * vy * vz / v2 * rz; + ret[3] = -gamma * vz * r0 + (1. + (gamma - 1.) * vz * vz / v2) * rz + + (gamma - 1.) * vz * vx / v2 * rx + (gamma - 1.) * vz * vy / v2 * ry; + + return ret; + } EventGeneratorBase::~EventGeneratorBase() { @@ -374,7 +397,8 @@ namespace thermalfist { else totals = GenerateTotalsGCE(); - double weight = ComputeWeight(totals); + double weight = ComputeWeightNew(totals); + //std::cout << weight << " " << ComputeWeightNew(totals) << "\n"; if (weight < 0.) continue; @@ -1187,7 +1211,8 @@ namespace thermalfist { std::vector momentum = m_MomentumGens[i]->GetMomentum(tmass); //std::vector momentum = m_MomentumGens[i]->GetMomentum(0.99999 * m_THM->TPS()->Particles()[i].Mass()); - primParticles[i].push_back(SimpleParticle(momentum[0], momentum[1], momentum[2], tmass, species.PdgId())); + primParticles[i].push_back(SimpleParticle(momentum[0], momentum[1], momentum[2], tmass, species.PdgId(), 0, + momentum[3], momentum[4], momentum[5], momentum[6])); } } @@ -1463,6 +1488,15 @@ namespace thermalfist { return ret; } + std::vector EventGeneratorBase::GCEMeanYields() const + { + std::vector ret = m_THM->Densities(); + for(size_t i = 0; i < ret.size(); ++i) { + ret[i] *= m_THM->Volume(); + } + return ret; + } + void EventGeneratorBase::SetVolume(double V) { if (m_Config.fEnsemble != EventGeneratorConfiguration::GCE) @@ -1561,16 +1595,17 @@ namespace thermalfist { double normweight = 1.; double weightev = 1.; bool fl = 1; - for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) { + int Nspecies = m_THM->TPS()->Particles().size(); + for (size_t i = 0; i < Nspecies; ++i) { double VVN = m_THM->Volume(); - for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) + for (size_t j = 0; j < Nspecies; ++j) VVN -= model->VirialCoefficient(j, i) * totals[j]; if (VVN < 0.) { fl = false; break; } double VVNev = m_THM->Volume(); - for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) + for (size_t j = 0; j < Nspecies; ++j) VVNev -= model->VirialCoefficient(j, i) * densities[j] * V; weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]); @@ -1647,6 +1682,192 @@ namespace thermalfist { return ret; } + double EventGeneratorBase::ComputeWeightNew(const std::vector& totals) const + { + // Compute the normlaized weight factor due to EV/vdW interactions + // If V - bN < 0, returns -1. + std::vector* densitiesid = NULL; + std::vector tmpdens; + const std::vector& densities = m_THM->Densities(); + if (m_THM->InteractionModel() != ThermalModelBase::Ideal) { + tmpdens = m_DensitiesIdeal; + densitiesid = &tmpdens; + } + + double ret = 1.; + + if (m_THM->InteractionModel() == ThermalModelBase::DiagonalEV) { + ThermalModelEVDiagonal* model = static_cast(m_THM); + double V = m_THM->Volume(); + double VVN = m_THM->Volume(); + + for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) + VVN -= model->ExcludedVolume(i) * totals[i]; + + if (VVN < 0.) + return -1.; + + double weight = 1.; + double logweight = 0.; + + double normweight = 1.; + double weightev = 1.; + double VVNev = m_THM->Volume(); + for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) + VVNev -= model->ExcludedVolume(i) * densities[i] * m_THM->Volume(); + + for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) { + weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]); + if (densitiesid->operator[](i) > 0. && densities[i] > 0.) + logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]); + + weightev *= pow(VVNev / V * densitiesid->operator[](i) / densities[i], densities[i] * V); + + if (densitiesid->operator[](i) > 0. && densities[i] > 0.) + //normweight *= pow(VVN / V, totals[i]) / pow(VVNev / V, densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V)); + normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V)); + } + + m_LastWeight = weight; + m_LastLogWeight = logweight; + m_LastNormWeight = normweight; + + ret = normweight; + } + + + if (m_THM->InteractionModel() == ThermalModelBase::CrosstermsEV) { + ThermalModelEVCrossterms* model = static_cast(m_THM); + double V = m_THM->Volume(); + + double weight = 1.; + double logweight = 0.; + double normweight = 1.; + double weightev = 1.; + bool fl = true; + int Nspecies = m_THM->TPS()->Particles().size(); + + int NEVcomp = model->EVComponentIndices().size(); + std::vector Nscomp(NEVcomp, 0); + std::vector Nevscomp(NEVcomp, 0.); + std::vector bns(NEVcomp, 0.), bnevs(NEVcomp, 0.), dmuTs(NEVcomp, 0.); + const std::vector< std::vector >& virial = model->VirialMatrix(); + + for (size_t icomp = 0; icomp < NEVcomp; ++icomp) { + const std::vector& indis = model->EVComponentIndices()[icomp]; + int Nlocal = indis.size(); + for (size_t ilocal = 0; ilocal < Nlocal; ++ilocal) { + int ip = indis[ilocal]; + Nscomp[icomp] += totals[ip]; + Nevscomp[icomp] += densities[ip] * V; + } + + if (indis.size()) { + int i1 = indis[0]; + + for (size_t j = 0; j < Nspecies; ++j) { + //bns[icomp] += model->VirialCoefficient(j, i1) * totals[j] / V; + //bnevs[icomp] += model->VirialCoefficient(j, i1) * densities[j]; + bns[icomp] += virial[j][i1] * totals[j];// / V; + bnevs[icomp] += virial[j][i1] * densities[j]; + } + bns[icomp] /= V; + + if (bns[icomp] > 1.) + fl = false; + + dmuTs[icomp] = model->DeltaMu(i1) / model->Parameters().T; + } + + normweight *= pow((1. - bns[icomp]) / (1. - bnevs[icomp]), Nscomp[icomp]) * exp(-dmuTs[icomp] * (Nscomp[icomp] - Nevscomp[icomp])); + + } + + //for (size_t i = 0; i < Nspecies; ++i) { + // double VVN = m_THM->Volume(); + + // for (size_t j = 0; j < Nspecies; ++j) + // VVN -= model->VirialCoefficient(j, i) * totals[j]; + + // if (VVN < 0.) { fl = false; break; } + + // double VVNev = m_THM->Volume(); + // for (size_t j = 0; j < Nspecies; ++j) + // VVNev -= model->VirialCoefficient(j, i) * densities[j] * V; + + // weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]); + // if (densitiesid->operator[](i) > 0. && densities[i] > 0.) + // logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]); + + // weightev *= pow(VVNev / m_THM->Volume() * densitiesid->operator[](i) / densities[i], densities[i] * V); + // if (densitiesid->operator[](i) > 0. && densities[i] > 0.) + // normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V)); + //} + + if (!fl) + return -1.; + + m_LastWeight = normweight; + m_LastLogWeight = log(normweight); + m_LastNormWeight = normweight; + + ret = normweight; + } + + if (m_THM->InteractionModel() == ThermalModelBase::QvdW) { + ThermalModelVDW* model = static_cast(m_THM); + double V = m_THM->Volume(); + + double weight = 1.; + double logweight = 0.; + double normweight = 1.; + double weightvdw = 1.; + bool fl = 1; + for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) { + double VVN = m_THM->Volume(); + + for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) + VVN -= model->VirialCoefficient(j, i) * totals[j]; + + if (VVN < 0.) { fl = false; break; } + + double VVNev = m_THM->Volume(); + for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) + VVNev -= model->VirialCoefficient(j, i) * densities[j] * V; + + weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]); + if (densitiesid->operator[](i) > 0. && densities[i] > 0.) + logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]); + + for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) { + double aij = model->AttractionCoefficient(i, j); + weight *= exp(aij * totals[j] / m_THM->Parameters().T / m_THM->Volume() * totals[i]); + logweight += totals[i] * aij * totals[j] / m_THM->Parameters().T / m_THM->Volume(); + } + + weightvdw *= pow(VVNev / m_THM->Volume() * densitiesid->operator[](i) / densities[i], densities[i] * V); + if (densitiesid->operator[](i) > 0. && densities[i] > 0.) + normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V)); + + for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) { + double aij = model->AttractionCoefficient(i, j); + weightvdw *= exp(aij * densities[j] / m_THM->Parameters().T * densities[i] * V); + normweight *= exp(aij * totals[j] / m_THM->Parameters().T / m_THM->Volume() * totals[i] - aij * densities[j] / m_THM->Parameters().T * densities[i] * V); + } + + } + if (!fl) + return -1.; + + m_LastWeight = weight; + m_LastLogWeight = logweight; + m_LastNormWeight = normweight; + + ret = normweight; + } + + return ret; + } EventGeneratorConfiguration::EventGeneratorConfiguration() { fEnsemble = GCE; diff --git a/src/library/HRGEventGenerator/ParticleDecaysMC.cpp b/src/library/HRGEventGenerator/ParticleDecaysMC.cpp index 60a3210..3a29418 100644 --- a/src/library/HRGEventGenerator/ParticleDecaysMC.cpp +++ b/src/library/HRGEventGenerator/ParticleDecaysMC.cpp @@ -77,6 +77,7 @@ namespace thermalfist { } std::vector TwoBodyDecay(const SimpleParticle & Mother, double m1, long long pdg1, double m2, long long pdg2) { + std::vector ret(0); ret.push_back(Mother); ret.push_back(Mother); @@ -103,6 +104,8 @@ namespace thermalfist { ret[1].pz = -ret[0].pz; ret[1].p0 = sqrt(m2*m2 + ret[1].px*ret[1].px + ret[1].py*ret[1].py + ret[1].pz*ret[1].pz); + double ten2 = ret[1].p0; + ret[0] = LorentzBoost(ret[0], -vx, -vy, -vz); ret[1] = LorentzBoost(ret[1], -vx, -vy, -vz); @@ -117,6 +120,14 @@ namespace thermalfist { printf("**WARNING** Issue in a two-body decay!\n"); } +#ifdef DEBUGDECAYS + if (abs(Mother.p0 - (ret[0].p0 + ret[1].p0)) > 1.e-9) { + printf("Two-body decay energy conservation issue: %lf %lf\n", + Mother.p0 - (ret[0].p0 + ret[1].p0), sqrt(vx*vx+vy*vy+vz*vz)); + printf("%lf %lf\n", Mother.m, ten1 + ten2); + } +#endif + return ret; } @@ -138,8 +149,15 @@ namespace thermalfist { SimpleParticle Mother2 = Mother; // Mass validation double tmasssum = 0.; - for (size_t i = 0; i < masses.size(); ++i) tmasssum += masses[i]; - if (Mother2.m < tmasssum) Mother2.m = tmasssum + 1e-7; + for (size_t i = 0; i < masses.size(); ++i) + tmasssum += masses[i]; + + // If sum of decay product masses larger than the mother particle mass (happens with zero-width resonances), + // adjust the mass and the energy of the mother particle + if (Mother2.m < tmasssum) { + Mother2.m = tmasssum + 1.e-7; + Mother2.p0 = sqrt(Mother2.px * Mother2.px + Mother2.py * Mother2.py + Mother2.pz * Mother2.pz + Mother2.m * Mother2.m); + } if (masses.size() == 2) return TwoBodyDecay(Mother2, masses[0], pdgs[0], masses[1], pdgs[1]); double tmin = 0.; diff --git a/src/library/HRGEventGenerator/RandomGenerators.cpp b/src/library/HRGEventGenerator/RandomGenerators.cpp index e984f0e..e62e26e 100644 --- a/src/library/HRGEventGenerator/RandomGenerators.cpp +++ b/src/library/HRGEventGenerator/RandomGenerators.cpp @@ -10,6 +10,7 @@ #include "HRGBase/xMath.h" #include "HRGEventGenerator/SimpleParticle.h" #include "HRGEventGenerator/ParticleDecaysMC.h" +#include "HRGEventGenerator/EventGeneratorBase.h" namespace thermalfist { @@ -189,6 +190,11 @@ namespace thermalfist { ret.push_back(tp*cos(tphi)*sthe); //px ret.push_back(tp*sin(tphi)*sthe); //py ret.push_back(tp*cthe); //pz + // TODO: proper Cartesian coordinates + ret.push_back(0.); //r0 + ret.push_back(0.); //rx + ret.push_back(0.); //ry + ret.push_back(0.); //rz return ret; } @@ -202,27 +208,28 @@ namespace thermalfist { return tp * tp / (texp + m_Statistics) / x; } - void ThermalMomentumGenerator::FixParameters() - { - double eps = 1e-8; - double l = 0., r = 1.; - double m1 = l + (r - l) / 3.; - double m2 = r - (r - l) / 3.; - int MAXITERS = 200; - int iter = 0; - while (fabs(m2 - m1) > eps && iter < MAXITERS) { - if (g(m1) < g(m2)) { - l = m1; - } - else { - r = m2; - } - m1 = l + (r - l) / 3.; - m2 = r - (r - l) / 3.; - iter++; - } - m_Max = g((m1 + m2) / 2.); - } + //void ThermalMomentumGenerator::FixParameters() + //{ + // //double eps = 1e-8; + // //double l = 0., r = 1.; + // //double m1 = l + (r - l) / 3.; + // //double m2 = r - (r - l) / 3.; + // //int MAXITERS = 200; + // //int iter = 0; + // //while (fabs(m2 - m1) > eps && iter < MAXITERS) { + // // if (g(m1) < g(m2)) { + // // l = m1; + // // } + // // else { + // // r = m2; + // // } + // // m1 = l + (r - l) / 3.; + // // m2 = r - (r - l) / 3.; + // // iter++; + // //} + // //m_Max = g((m1 + m2) / 2.); + // m_Max = ComputeMaximum(m_Mass); + //} double ThermalMomentumGenerator::ComputeMaximum(double mass) const { @@ -293,24 +300,44 @@ namespace thermalfist { if (mass < 0.) mass = Mass(); - std::vector ret(3, 0.); - while (1) { - double zetacand = GetRandomZeta(RandomGenerators::randgenMT); - double eta = -EtaMax() + 2. * EtaMax() * RandomGenerators::randgenMT.rand(); - double ph = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand(); - double betar = m_FreezeoutModel->tanhetaperp(zetacand); - double cosheta = cosh(eta); - double sinheta = sinh(eta); + double zetacand = GetRandomZeta(RandomGenerators::randgenMT); + double eta = -EtaMax() + 2. * EtaMax() * RandomGenerators::randgenMT.rand(); + double ph = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand(); + + double betar = m_FreezeoutModel->tanhetaperp(zetacand); + double cosheta = cosh(eta); + double sinheta = sinh(eta); - double cosphi = cos(ph); - double sinphi = sin(ph); + double cosphi = cos(ph); + double sinphi = sin(ph); - double vx = betar * cosphi / cosheta; - double vy = betar * sinphi / cosheta; - double vz = tanh(eta); + double vx = betar * cosphi / cosheta; + double vy = betar * sinphi / cosheta; + double vz = tanh(eta); - SimpleParticle part(0., 0., 0., mass, 0); + double dRdZeta = m_FreezeoutModel->dRdZeta(zetacand); + double dtaudZeta = m_FreezeoutModel->dtaudZeta(zetacand); + + double coshetaperp = m_FreezeoutModel->coshetaperp(zetacand); + double sinhetaperp = m_FreezeoutModel->sinhetaperp(zetacand); + + std::vector dsigma_lab; + dsigma_lab.push_back(dRdZeta * cosheta); + dsigma_lab.push_back(dtaudZeta * cosphi); + dsigma_lab.push_back(dtaudZeta * sinphi); + dsigma_lab.push_back(dRdZeta * sinheta); + + // dsigma^\mu in the local rest frame + std::vector dsigma_loc = LorentzBoost(dsigma_lab, vx, vy, vz); + + // Maximum weight for the rejection sampling of the momentum + double maxWeight = 1. + std::abs(dsigma_loc[1] / dsigma_loc[0]) + std::abs(dsigma_loc[2] / dsigma_loc[0]) + std::abs(dsigma_loc[3] / dsigma_loc[0]); + maxWeight += 1.e-5; // To stabilize models where maxWeight is equal to unity, like Cracow model + + SimpleParticle part(0., 0., 0., mass, 0); + + while (true) { double tp = m_Generator.GetP(mass); double tphi = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand(); @@ -323,30 +350,48 @@ namespace thermalfist { double p0LRF = part.p0; - if (betar != 0.0 || eta != 0.0) - part = ParticleDecaysMC::LorentzBoost(part, -vx, -vy, -vz); + double dsigmamu_pmu_loc = dsigma_loc[0] * part.p0 + - dsigma_loc[1] * part.px - dsigma_loc[2] * part.py - dsigma_loc[3] * part.pz; - //double prob = (cosheta * part.p0 - sinheta * part.pz) / - // (2. * (1. / (1 - betar * betar)) * (cosheta * part.p0 - sinheta * part.pz - betar * (part.px * cos(ph) + part.py * sin(ph)))); - double dRdZeta = m_FreezeoutModel->dRdZeta(zetacand); - double dtaudZeta = m_FreezeoutModel->dtaudZeta(zetacand); + double dsigmamu_umu_loc = dsigma_loc[0]; - double coshetaperp = m_FreezeoutModel->coshetaperp(zetacand); - double sinhetaperp = m_FreezeoutModel->sinhetaperp(zetacand); + double dumu_pmu_loc = p0LRF; - double prob = (dRdZeta * (cosheta * part.p0 - sinheta * part.pz) - - dtaudZeta * (cosphi * part.px + sinphi * part.py)) / - (coshetaperp * dRdZeta - sinhetaperp * dtaudZeta) / p0LRF - / 2.; + double Weight = dsigmamu_pmu_loc / dsigmamu_umu_loc / dumu_pmu_loc / maxWeight; - if (RandomGenerators::randgenMT.rand() < prob) { - ret[0] = part.px; - ret[1] = part.py; - ret[2] = part.pz; - break; + if (Weight > 1.) { + printf("**WARNING** BoostInvariantHypersurfaceMomentumGenerator::GetMomentum: Weight exceeds unity by %E\n", + Weight - 1.); } + + if (RandomGenerators::randgenMT.rand() < Weight) + break; + } + + if (betar != 0.0 || eta != 0.0) + part = ParticleDecaysMC::LorentzBoost(part, -vx, -vy, -vz); + + + std::vector ret(3, 0.); + ret[0] = part.px; + ret[1] = part.py; + ret[2] = part.pz; + + // Space-time coordinates + double tau = m_FreezeoutModel->taufunc(zetacand); + double r0 = tau * cosheta; + double rz = tau * sinheta; + + double Rperp = m_FreezeoutModel->Rfunc(zetacand); + double rx = Rperp * cosphi; + double ry = Rperp * sinphi; + + ret.push_back(r0); + ret.push_back(rx); + ret.push_back(ry); + ret.push_back(rz); return ret; } @@ -799,11 +844,17 @@ namespace thermalfist { if (GetBeta() != 0.0) part = ParticleDecaysMC::LorentzBoost(part, -vx, -vy, -vz); - std::vector ret(3); + std::vector ret(7); ret[0] = part.px; ret[1] = part.py; ret[2] = part.pz; + // Assume unit sphere at t = 0 + ret[3] = 0.; + ret[4] = sinth * cos(ph); + ret[5] = sinth * sin(ph); + ret[6] = costh; + return ret; } diff --git a/src/library/HRGEventGenerator/SimpleEvent.cpp b/src/library/HRGEventGenerator/SimpleEvent.cpp index 9c9778b..569d03b 100644 --- a/src/library/HRGEventGenerator/SimpleEvent.cpp +++ b/src/library/HRGEventGenerator/SimpleEvent.cpp @@ -1,6 +1,6 @@ /* * Thermal-FIST package - * + * * Copyright (c) 2015-2019 Volodymyr Vovchenko * * GNU General Public License (GPLv3 or later) @@ -11,18 +11,27 @@ namespace thermalfist { - void SimpleEvent::writeToFile(std::ofstream & fout, const EventOutputConfig& config, int eventnumber) + void SimpleEvent::writeToFile(std::ofstream& fout, const EventOutputConfig& config, int eventnumber) { fout << "Event " << eventnumber << std::endl; - fout << "Weight: " << weight << std::endl; - fout << std::setw(20) << "pdgid" - << std::setw(20) << "px[GeV]" - << std::setw(20) << "py[GeV]" - << std::setw(20) << "pz[GeV]"; + if (config.printWeight) + fout << "Weight: " << weight << std::endl; + + fout << std::setw(20) << "pdgid"; + + if (config.printCoordinates) + fout << std::setw(20) << "r0[fm/c]" + << std::setw(20) << "rx[fm]" + << std::setw(20) << "ry[fm]" + << std::setw(20) << "rz[fm]"; if (config.printEnergy) - fout << std::setw(20) << "p0[GeV]"; + fout << std::setw(20) << "p0[GeV/c2]"; + + fout << std::setw(20) << "px[GeV/c]" + << std::setw(20) << "py[GeV/c]" + << std::setw(20) << "pz[GeV/c]"; if (config.printMotherPdg) fout << std::setw(20) << "mother_pdgid"; @@ -32,9 +41,20 @@ namespace thermalfist { fout << std::endl; + fout.precision(10); + fout << std::scientific; + for (size_t i = 0; i < Particles.size(); ++i) { - fout << std::setw(20) << Particles[i].PDGID - << std::setw(20) << Particles[i].px + fout << std::setw(20) << Particles[i].PDGID; + + if (config.printCoordinates) + fout << std::setw(20) << Particles[i].r0 + << std::setw(20) << Particles[i].rx + << std::setw(20) << Particles[i].ry + << std::setw(20) << Particles[i].rz; + + + fout << std::setw(20) << Particles[i].px << std::setw(20) << Particles[i].py << std::setw(20) << Particles[i].pz; @@ -50,12 +70,18 @@ namespace thermalfist { fout << std::endl; } - fout << std::scientific; if (config.printPhotonsLeptons) { for (size_t i = 0; i < PhotonsLeptons.size(); ++i) { - fout << std::setw(20) << PhotonsLeptons[i].PDGID - << std::setw(20) << PhotonsLeptons[i].px + fout << std::setw(20) << PhotonsLeptons[i].PDGID; + + if (config.printCoordinates) + fout << std::setw(20) << PhotonsLeptons[i].r0 + << std::setw(20) << PhotonsLeptons[i].rx + << std::setw(20) << PhotonsLeptons[i].ry + << std::setw(20) << PhotonsLeptons[i].rz; + + fout << std::setw(20) << PhotonsLeptons[i].px << std::setw(20) << PhotonsLeptons[i].py << std::setw(20) << PhotonsLeptons[i].pz; @@ -75,6 +101,35 @@ namespace thermalfist { fout << std::fixed; } + void SimpleEvent::writeToFileForUrqmd(std::ofstream& fout) + { + fout << "# " << Particles.size() << std::endl; + + fout.precision(16); + fout << std::scientific; + + const int tabsize = 23; + + for (size_t i = 0; i < Particles.size(); ++i) { + fout << std::setw(12) << Particles[i].PDGID << " "; + + fout << std::setw(tabsize) << Particles[i].r0 << " " + << std::setw(tabsize) << Particles[i].rx << " " + << std::setw(tabsize) << Particles[i].ry << " " + << std::setw(tabsize) << Particles[i].rz << " "; + + + fout << std::setw(tabsize) << Particles[i].p0 << " "; + fout << std::setw(tabsize) << Particles[i].px << " " + << std::setw(tabsize) << Particles[i].py << " " + << std::setw(tabsize) << Particles[i].pz << " "; + + fout << std::endl; + } + + fout << std::fixed; + } + void SimpleEvent::RapidityBoost(double dY) { for (size_t i = 0; i < Particles.size(); ++i) @@ -83,7 +138,7 @@ namespace thermalfist { AllParticles[i].RapidityBoost(dY); } - SimpleEvent SimpleEvent::MergeEvents(const SimpleEvent & evt1, const SimpleEvent & evt2) + SimpleEvent SimpleEvent::MergeEvents(const SimpleEvent& evt1, const SimpleEvent& evt2) { SimpleEvent ret; ret.Particles.reserve(evt1.Particles.size() + evt2.Particles.size()); diff --git a/src/library/HRGPCE/ThermalModelPCE.cpp b/src/library/HRGPCE/ThermalModelPCE.cpp index ad4a2b5..aa34571 100644 --- a/src/library/HRGPCE/ThermalModelPCE.cpp +++ b/src/library/HRGPCE/ThermalModelPCE.cpp @@ -29,7 +29,7 @@ namespace thermalfist { // A helper ThermalParticleSystem instance to compute the PCE effective charges for all particles ThermalParticleSystem TPShelper = ThermalParticleSystem(*m_model->TPS()); - // Set the nucleon content of the known light nuclei as "decay" products + // Set the nucleon content of the known stable light nuclei as "decay" products PrepareNucleiForPCE(&TPShelper); m_StabilityFlags = StabilityFlags; diff --git a/src/library/HRGVDW/ThermalModelVDW.cpp b/src/library/HRGVDW/ThermalModelVDW.cpp index 5e41e0f..4ad1ac2 100644 --- a/src/library/HRGVDW/ThermalModelVDW.cpp +++ b/src/library/HRGVDW/ThermalModelVDW.cpp @@ -290,6 +290,7 @@ namespace thermalfist { } vector sol = SearchSingleSolution(curmust); + bool fl = true; for(size_t i = 0; i < sol.size(); ++i) if (sol[i] != sol[i]) fl = false; @@ -1135,6 +1136,7 @@ namespace thermalfist { double maxdiff = 0.; for (size_t i = 0; i < xdelta.size(); ++i) { maxdiff = std::max(maxdiff, fabs(xdelta[i])); + maxdiff = std::max(maxdiff, fabs(f[i])); } return (maxdiff < m_MaximumError); }