diff --git a/CHANGELOG.md b/CHANGELOG.md index 610b3d4..45b9664 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,13 @@ 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.4.1] + +Date: 2022-10-18 + +- GUI: Added more clarity concerning the loading of the decays list +- GUI: Number of known decays channels shown in the table for unstable hadrons +- Added explicit calculation of Fermi-Dirac integrals at T = 0 ## [Version 1.4] diff --git a/CMakeLists.txt b/CMakeLists.txt index 77be144..fa7b0bb 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 4) -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/IdealGasFunctions.h b/include/HRGBase/IdealGasFunctions.h index 681b43f..01de971 100644 --- a/include/HRGBase/IdealGasFunctions.h +++ b/include/HRGBase/IdealGasFunctions.h @@ -28,7 +28,9 @@ namespace thermalfist { * \brief Identifies the thermodynamic function. * */ - enum Quantity { ParticleDensity, EnergyDensity, EntropyDensity, Pressure, chi2, chi3, chi4, ScalarDensity }; + enum Quantity { ParticleDensity, EnergyDensity, EntropyDensity, Pressure, chi2, chi3, chi4, ScalarDensity, + chi2difull, chi3difull, chi4difull + }; /** * \brief Identifies whether quantum statistics * are to be computed using the cluster expansion @@ -109,11 +111,12 @@ namespace thermalfist { double BoltzmannTdndmu(int N, double T, double mu, double m, double deg); /** - * \brief Computes the n-th order susceptibility for a Maxwell-Boltzmann gas. + * \brief Computes the n-th order dimensionless susceptibility for a Maxwell-Boltzmann gas. * * Computes \f$ \chi_n \equiv \frac{\partial^n p/T^4}{\partial (mu/T)^n} \f$ * for a Maxwell-Boltzmann gas. * + * \param N Susceptibility order. * \param T Temperature [GeV]. * \param mu Chemical potential [GeV]. * \param m Particle's mass [GeV]. @@ -122,6 +125,21 @@ namespace thermalfist { */ double BoltzmannChiN(int N, double T, double mu, double m, double deg); + /** + * \brief Computes the n-th order dimensionfull susceptibility for a Maxwell-Boltzmann gas. + * + * Computes \f$ \chi_n \equiv \frac{\partial^n p}{\partial mu^n} \f$ + * for a Maxwell-Boltzmann gas. + * + * \param N Susceptibility order. + * \param T Temperature [GeV]. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return \f$ \chi_n \f$. + */ + double BoltzmannChiNDimensionfull(int N, double T, double mu, double m, double deg); + /** * \brief Computes the particle number density of a quantum ideal gas using cluster expansion. * @@ -191,11 +209,12 @@ namespace thermalfist { double QuantumClusterExpansionTdndmu(int N, int statistics, double T, double mu, double m, double deg, int order = 1); /** - * \brief Computes the n-th order susceptibility for a quantum ideal gas using cluster expansion. + * \brief Computes the n-th order dimensionless susceptibility for a quantum ideal gas using cluster expansion. * * Computes \f$ \chi_n \equiv \frac{\partial^n p/T^4}{\partial (mu/T)^n} \f$ * for a quantum ideal gas using cluster expansion. * + * \param N Susceptibility order. * \param T Temperature [GeV]. * \param mu Chemical potential [GeV]. * \param m Particle's mass [GeV]. @@ -204,6 +223,21 @@ namespace thermalfist { */ double QuantumClusterExpansionChiN(int N, int statistics, double T, double mu, double m, double deg, int order = 1); + /** + * \brief Computes the n-th order dimensionfull susceptibility for a quantum ideal gas using cluster expansion. + * + * Computes \f$ \chi_n \equiv \frac{\partial^n p}{\partial mu^n} \f$ + * for a quantum ideal gas using cluster expansion. + * + * \param N Susceptibility order. + * \param T Temperature [GeV]. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return \f$ \chi_n \f$ [GeV^{4-N}]. + */ + double QuantumClusterExpansionChiNDimensionfull(int N, int statistics, double T, double mu, double m, double deg, int order = 1); + /** * \brief Computes the particle number density of a quantum ideal gas using 32-point Gauss-Laguerre quadratures. * @@ -265,11 +299,12 @@ namespace thermalfist { double QuantumNumericalIntegrationTdndmu(int N, int statistics, double T, double mu, double m, double deg); /** - * \brief Computes the n-th order susceptibility for a quantum ideal gas using 32-point Gauss-Laguerre quadratures. + * \brief Computes the n-th order dimensionless susceptibility for a quantum ideal gas using 32-point Gauss-Laguerre quadratures. * * Computes \f$ \chi_n \equiv \frac{\partial^n p/T^4}{\partial (mu/T)^n} \f$ * for a quantum ideal gas using 32-point Gauss-Laguerre quadratures. * + * \param N Susceptibility order. * \param T Temperature [GeV]. * \param mu Chemical potential [GeV]. * \param m Particle's mass [GeV]. @@ -278,6 +313,21 @@ namespace thermalfist { */ double QuantumNumericalIntegrationChiN(int N, int statistics, double T, double mu, double m, double deg); + /** + * \brief Computes the n-th order dimensionfull susceptibility for a quantum ideal gas using 32-point Gauss-Laguerre quadratures. + * + * Computes \f$ \chi_n \equiv \frac{\partial^n p}{\partial mu^n} \f$ + * for a quantum ideal gas using 32-point Gauss-Laguerre quadratures. + * + * \param N Susceptibility order. + * \param T Temperature [GeV]. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return \f$ \chi_n \f$ [GeV^{n-4}. + */ + double QuantumNumericalIntegrationChiNDimensionfull(int N, int statistics, double T, double mu, double m, double deg); + /** * \brief Auxiliary function. * @@ -365,12 +415,13 @@ namespace thermalfist { double FermiNumericalIntegrationLargeMuTdndmu(int N, double T, double mu, double m, double deg); /** - * \brief Computes the n-th order susceptibility for a Fermi-Dirac ideal gas + * \brief Computes the n-th order dimensionless susceptibility for a Fermi-Dirac ideal gas * at mu > m. * * Computes \f$ \chi_n \equiv \frac{\partial^n p/T^4}{\partial (mu/T)^n} \f$ * for a Fermi-Dirac ideal gas at mu > m. * + * \param N Susceptibility order. * \param T Temperature [GeV]. * \param mu Chemical potential [GeV]. * \param m Particle's mass [GeV]. @@ -379,6 +430,108 @@ namespace thermalfist { */ double FermiNumericalIntegrationLargeMuChiN(int N, double T, double mu, double m, double deg); + /** + * \brief Computes the n-th order dimensionfull susceptibility for a Fermi-Dirac ideal gas + * at mu > m. + * + * Computes \f$ \chi_n \equiv \frac{\partial^n p}{\partial mu^n} \f$ + * for a Fermi-Dirac ideal gas at mu > m. + * + * \param N Susceptibility order. + * \param T Temperature [GeV]. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return \f$ \chi_n \f$. + */ + double FermiNumericalIntegrationLargeMuChiNDimensionfull(int N, double T, double mu, double m, double deg); + + + /** + * \brief Computes the particle number density of a Fermi-Dirac ideal gas at zero temperature. + * + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return Particle number density [fm-3]. + */ + double FermiZeroTDensity(double mu, double m, double deg); + + /** + * \brief Computes the pressure of a Fermi-Dirac ideal gas at zero temperature. + * + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return Pressure [GeV fm-3]. + */ + double FermiZeroTPressure(double mu, double m, double deg); + + /** + * \brief Computes the energy density of a Fermi-Dirac ideal gas at zero temperature. + * + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return Energy density [GeV fm-3]. + */ + double FermiZeroTEnergyDensity(double mu, double m, double deg); + + /** + * \brief Computes the entropy density of a Fermi-Dirac ideal gas at zero temperature. + * + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return Entropy density [GeV fm-3]. + */ + double FermiZeroTEntropyDensity(double mu, double m, double deg); + + /** + * \brief Computes the scalar density of a Fermi-Dirac ideal gas at zero temperature. + * + * \param T Temperature [GeV]. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return Scalar density [fm-3]. + */ + double FermiZeroTScalarDensity(double mu, double m, double deg); // TODO: Check for correctness + + double FermiZeroTdn1dmu1(double mu, double m, double deg); + double FermiZeroTdn2dmu2(double mu, double m, double deg); + double FermiZeroTdn3dmu3(double mu, double m, double deg); + double FermiZeroTdndmu(int N, double mu, double m, double deg); + + /** + * \brief Computes the n-th order susceptibility for a Fermi-Dirac ideal gas at zero temperature. + * + * Computes \f$ \chi_n \equiv \frac{\partial^n p/T^4}{\partial (mu/T)^n} \f$ + * for a Fermi-Dirac ideal gas at mu > m. + * + * \param N Susceptibility order. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return \f$ \chi_n \f$. + */ + double FermiZeroTChiN(int N, double mu, double m, double deg); + + /** + * \brief Computes the n-th order susceptibility scaled by T^n for a Fermi-Dirac ideal gas at zero temperature. + * + * Computes \f$ \chi_n \equiv \frac{\partial^n p}{\partial (mu)^n} \f$ + * for a Fermi-Dirac ideal gas at mu > m. + * Useful for zero temperature where the dimensionless susceptibility would be singular. + * + * \param N Susceptibility order. + * \param mu Chemical potential [GeV]. + * \param m Particle's mass [GeV]. + * \param deg Internal degeneracy factor. + * \return \f$ \chi_n \f$ [GeV^{4-N}]. + */ + double FermiZeroTChiNDimensionfull(int N, double mu, double m, double deg); + /** * \brief Calculation of a generic ideal gas function. * diff --git a/include/HRGBase/ThermalModelBase.h b/include/HRGBase/ThermalModelBase.h index c5485ef..f10cc6f 100644 --- a/include/HRGBase/ThermalModelBase.h +++ b/include/HRGBase/ThermalModelBase.h @@ -422,6 +422,15 @@ namespace thermalfist { */ double ChemicalPotential(int i) const; + /** + * \brief Sets the chemical potential of particle species i. + * + * + * \param i 0-based index of particle species + * \param chem value of the chemical potential + */ + virtual void SetChemicalPotential(int i, double chem); + /** * \brief Chemical potential entering the ideal gas expressions of particle species i. * diff --git a/include/HRGBase/ThermalParticle.h b/include/HRGBase/ThermalParticle.h index c047b08..3c32747 100644 --- a/include/HRGBase/ThermalParticle.h +++ b/include/HRGBase/ThermalParticle.h @@ -216,6 +216,22 @@ namespace thermalfist { */ double chi(int index, const ThermalModelParameters ¶ms, bool useWidth = 0, double mu = 0.) const; + /** + * \brief Computes the ideal gas dimensionfull susceptibility \f$ \chi_n \equiv \frac{\partial^n p}{\partial mu^n} \f$. + * + * Computes the dimensionfull susceptibility \f$ \chi_n \equiv \frac{\partial^n p}{\partial mu^n} \f$ + * of the corresponding ideal gas. + * Takes into account chemical non-equilibrium fugacity factors + * and finite resonance widths. + * + * \param index Order of the susceptibility. + * \param params Structure containing the temperature value and the chemical factors. + * \param useWidth Whether finite widths are taken into account. + * \param mu Chemical potential. + * \return Value of the computed susceptility [GeV^{4-n]. + */ + double chiDimensionfull(int index, const ThermalModelParameters& params, bool useWidth = 0, double mu = 0.) const; + /** * \brief Computes the scaled variance of particle number fluctuations * in the ideal gas. diff --git a/include/HRGPCE/ThermalModelPCE.h b/include/HRGPCE/ThermalModelPCE.h index 09ccaf0..d42f0c0 100644 --- a/include/HRGPCE/ThermalModelPCE.h +++ b/include/HRGPCE/ThermalModelPCE.h @@ -154,6 +154,7 @@ namespace thermalfist { * \return Vector of chemical potentials of all particles species, as resulted from the last CalculatePCE() call */ const std::vector& ChemicalPotentials() const { return m_ChemCurrent; } + std::vector& ChemicalPotentials() { return m_ChemCurrent; } /** * \return The system volume, as resulted from the last CalculatePCE() call @@ -187,6 +188,9 @@ namespace thermalfist { * This is only necessary if energy-dependent Breit-Wigner widths are used. */ void ApplyFixForBoseCondensation(); + + const std::vector< std::vector >& EffectiveCharges() const { return m_EffectiveCharges; } + const std::vector& StableMapTo() const { return m_StableMapTo; } protected: diff --git a/input/list/PDG2020_modular/list-charged-leptons.dat b/input/list/PDG2020_modular/list-charged-leptons.dat new file mode 100644 index 0000000..2cec40a --- /dev/null +++ b/input/list/PDG2020_modular/list-charged-leptons.dat @@ -0,0 +1,5 @@ +# leptons +# pdgid name stable mass[GeV] degeneracy statistics B Q S C |S| |C| width[GeV] threshold[GeV] + 11 e- 1 5.109989461E-04 2 1 0 -1 0 0 0 0 0 0 + 13 mu- 1 1.056583745E-01 2 1 0 -1 0 0 0 0 0 0 + 15 tau- 1 1.77686E+00 2 1 0 -1 0 0 0 0 0 0 diff --git a/src/gui/QtThermalFIST/mainwindow.cpp b/src/gui/QtThermalFIST/mainwindow.cpp index 2d78b45..6bdbe97 100644 --- a/src/gui/QtThermalFIST/mainwindow.cpp +++ b/src/gui/QtThermalFIST/mainwindow.cpp @@ -48,7 +48,7 @@ MainWindow::MainWindow(QWidget *parent) leList = new QLineEdit("");//QApplication::applicationDirPath()); leList->setReadOnly(true); if (TPS->Particles().size() > 0) - leList->setText(listpath); + leList->setText(listpath + " + decays.dat"); buttonLoad = new QPushButton(tr("Load particle list...")); connect(buttonLoad, SIGNAL(clicked()), this, SLOT(loadList())); @@ -157,15 +157,30 @@ void MainWindow::loadDecays() QString listpathprefix = QString(ThermalFIST_INPUT_FOLDER) + "/list"; if (leList->text().size() != 0) listpathprefix = leList->text(); - QString path = QFileDialog::getOpenFileName(this, tr("Open file with decays"), listpathprefix); - if (path.length() > 0) + QStringList pathdecays = QFileDialog::getOpenFileNames(this, tr("Open file with decays"), listpathprefix); + if (pathdecays.length() > 0) { - TPS->LoadDecays(path.toStdString()); + std::vector decayspaths; + for (int i = 0; i < pathdecays.length(); ++i) + decayspaths.push_back(pathdecays[i].toStdString()); + TPS->LoadDecays(decayspaths); model->ChangeTPS(TPS); tab1->resetTPS(); tab2->resetTPS(); tab5->resetTPS(); tabEditor->resetTPS(); + + leList->setText(clists); + + if (QFileInfo(QString::fromStdString(decayspaths[0])).dir().absolutePath() != + QFileInfo(clists).dir().absolutePath()) + leList->setText(leList->text() + " + " + QString::fromStdString(decayspaths[0])); + else + leList->setText(leList->text() + " + " + QFileInfo(QString::fromStdString(decayspaths[0])).fileName()); + + for (int idec = 1; idec < decayspaths.size(); ++idec) { + leList->setText(leList->text() + " + " + QFileInfo(QString::fromStdString(decayspaths[idec])).fileName()); + } } } @@ -243,22 +258,54 @@ void MainWindow::loadList() // { ThermalParticleSystem::flag_noexcitednuclei, ThermalParticleSystem::flag_nonuclei } //); - QString decpath = QFileInfo(pathlist[0]).absolutePath() + "/decays.dat"; + std::vector decays; + QStringList decpath; + decpath.push_back(QFileInfo(pathlist[0]).absolutePath() + "/decays.dat"); //if (!TPS->CheckDecayChannelsAreSpecified() && !QFileInfo(decpath).exists()) { - if (!QFileInfo(decpath).exists()) { - decpath = QFileDialog::getOpenFileName(this, tr("Open file with decays"), decpath); - if (decpath.length() > 0) - { - TPS->LoadDecays(decpath.toStdString()); + if (!QFileInfo(decpath[0]).exists()) { + + QMessageBox::StandardButton reply; + reply = QMessageBox::question(this, + "Decays", + "Decays file was not found at `decays.dat`. Would you like to load decays from another file?", + QMessageBox::Yes | QMessageBox::No); + if (reply == QMessageBox::Yes) { + decpath = QFileDialog::getOpenFileNames(this, tr("Open file with decays"), decpath[0]); + if (decpath.length() > 0) + { + for (int i = 0; i < decpath.length(); ++i) + decays.push_back(decpath[i].toStdString()); + } } } - std::vector decays; - decays.push_back(decpath.toStdString()); + else { + decays.push_back(decpath[0].toStdString()); + } *TPS = ThermalParticleSystem(paths, decays); //TPS->SetSortMode(ThermalParticleSystem::SortByBaryonAndMassAndPDG); model->ChangeTPS(TPS); leList->setText(pathlist[0]); + if (pathlist.size() > 1) { + for (int il = 1; il < pathlist.size(); ++il) { + leList->setText(leList->text() + " + " + QFileInfo(pathlist[il]).fileName()); + } + } + clists = leList->text(); + + if (decays.size() > 0) { + + if (QFileInfo(QString::fromStdString(decays[0])).dir().absolutePath() != + QFileInfo(clists).dir().absolutePath()) + leList->setText(leList->text() + " + " + QString::fromStdString(decays[0])); + else + leList->setText(leList->text() + " + " + QFileInfo(QString::fromStdString(decays[0])).fileName()); + + for (int idec = 1; idec < decays.size(); ++idec) { + leList->setText(leList->text() + " + " + QFileInfo(QString::fromStdString(decays[idec])).fileName()); + } + } + tab1->resetTPS(); tab2->resetTPS(); diff --git a/src/gui/QtThermalFIST/mainwindow.h b/src/gui/QtThermalFIST/mainwindow.h index ce41892..5904bd6 100644 --- a/src/gui/QtThermalFIST/mainwindow.h +++ b/src/gui/QtThermalFIST/mainwindow.h @@ -50,7 +50,8 @@ class MainWindow : public QMainWindow thermalfist::ThermalParticleSystem *TPS; thermalfist::ThermalModelBase *model; - QString cpath = ""; + QString cpath = ""; + QString clists = ""; public: MainWindow(QWidget *parent = 0); diff --git a/src/gui/QtThermalFIST/modeltab.cpp b/src/gui/QtThermalFIST/modeltab.cpp index 0ad8bb1..0506a5a 100644 --- a/src/gui/QtThermalFIST/modeltab.cpp +++ b/src/gui/QtThermalFIST/modeltab.cpp @@ -123,14 +123,14 @@ ModelTab::ModelTab(QWidget *parent, ThermalModelBase *modelop) layParameters->setAlignment(Qt::AlignLeft); QLabel *labelTemperature = new QLabel(tr("T (MeV):")); spinTemperature = new QDoubleSpinBox(); - spinTemperature->setMinimum(1.); + spinTemperature->setMinimum(0.); spinTemperature->setMaximum(10000.); spinTemperature->setValue(model->Parameters().T * 1e3); spinTemperature->setToolTip(tr("Temperature")); QLabel *labelmuB = new QLabel(tr("μB (MeV):")); spinmuB = new QDoubleSpinBox(); - spinmuB->setMinimum(-1000.); - spinmuB->setMaximum(2000.); + spinmuB->setMinimum(-10000.); + spinmuB->setMaximum(10000.); spinmuB->setValue(model->Parameters().muB * 1e3); spinmuB->setToolTip(tr("Baryochemical potential")); QLabel *labelgammaq = new QLabel(tr("γq:")); diff --git a/src/gui/QtThermalFIST/tablemodel.cpp b/src/gui/QtThermalFIST/tablemodel.cpp index bb10837..2374df0 100644 --- a/src/gui/QtThermalFIST/tablemodel.cpp +++ b/src/gui/QtThermalFIST/tablemodel.cpp @@ -49,6 +49,8 @@ QVariant TableModel::data(const QModelIndex &index, int role) const if (col==2) return model->TPS()->Particles()[RowToParticle[row]].Mass(); if ((col==3 && model->TPS()->Particles()[RowToParticle[row]].IsStable()) || (col==4 && model->TPS()->Particles()[RowToParticle[row]].IsNeutral()!=0)) return "*"; + if (col == 3 && !model->TPS()->Particles()[RowToParticle[row]].IsStable()) + return QString("%1 decays").arg(model->TPS()->Particles()[RowToParticle[row]].Decays().size()); if (col==5) { if (model->TPS()->Particles()[RowToParticle[row]].BaryonCharge()>0) return "+" + QString::number(model->TPS()->Particles()[RowToParticle[row]].BaryonCharge()); else if (model->TPS()->Particles()[RowToParticle[row]].BaryonCharge()<0) return QString::number(model->TPS()->Particles()[RowToParticle[row]].BaryonCharge()); diff --git a/src/library/HRGBase/IdealGasFunctions.cpp b/src/library/HRGBase/IdealGasFunctions.cpp index f357be1..197d11b 100644 --- a/src/library/HRGBase/IdealGasFunctions.cpp +++ b/src/library/HRGBase/IdealGasFunctions.cpp @@ -61,6 +61,11 @@ namespace thermalfist { return BoltzmannTdndmu(N - 1, T, mu, m, deg) / pow(T, 3) / xMath::GeVtoifm3(); } + double BoltzmannChiNDimensionfull(int N, double T, double mu, double m, double deg) + { + return BoltzmannTdndmu(N - 1, T, mu, m, deg) / pow(T, N - 1) / xMath::GeVtoifm3(); + } + double QuantumClusterExpansionDensity(int statistics, double T, double mu, double m, double deg, int order) { double sign = 1.; @@ -192,22 +197,29 @@ namespace thermalfist { } + double QuantumClusterExpansionChiNDimensionfull(int N, int statistics, double T, double mu, double m, double deg, int order) + { + return QuantumClusterExpansionTdndmu(N - 1, statistics, T, mu, m, deg, order) / pow(T, N - 1) / xMath::GeVtoifm3(); + } + // Gauss-Legendre 32-point quadrature for [0,1] interval - const double *legx32 = NumericalIntegration::coefficients_xleg32_zeroone; - const double *legw32 = NumericalIntegration::coefficients_wleg32_zeroone; + const double* legx32 = NumericalIntegration::coefficients_xleg32_zeroone; + const double* legw32 = NumericalIntegration::coefficients_wleg32_zeroone; // Gauss-Laguerre 32-point quadrature for [0,\infty] interval - const double *lagx32 = NumericalIntegration::coefficients_xlag32; - const double *lagw32 = NumericalIntegration::coefficients_wlag32; + const double* lagx32 = NumericalIntegration::coefficients_xlag32; + const double* lagw32 = NumericalIntegration::coefficients_wlag32; double QuantumNumericalIntegrationDensity(int statistics, double T, double mu, double m, double deg) { if (statistics == 0) return BoltzmannDensity(T, mu, m, deg); + if (statistics == 1 && T == 0.) return FermiZeroTDensity(mu, m, deg); if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuDensity(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationDensity: Bose-Einstein condensation, mass = %lf, mu = %lf\n", m, mu); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; double ret = 0.; double moverT = m / T; @@ -225,12 +237,14 @@ namespace thermalfist { double QuantumNumericalIntegrationPressure(int statistics, double T, double mu, double m, double deg) { if (statistics == 0) return BoltzmannPressure(T, mu, m, deg); + if (statistics == 1 && T == 0.) return FermiZeroTPressure(mu, m, deg); if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuPressure(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationPressure: Bose-Einstein condensation\n"); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; double ret = 0.; double moverT = m / T; @@ -250,12 +264,14 @@ namespace thermalfist { double QuantumNumericalIntegrationEnergyDensity(int statistics, double T, double mu, double m, double deg) { if (statistics == 0) return BoltzmannEnergyDensity(T, mu, m, deg); + if (statistics == 1 && T == 0.) return FermiZeroTEnergyDensity(mu, m, deg); if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuEnergyDensity(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationEnergyDensity: Bose-Einstein condensation\n"); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; double ret = 0.; double moverT = m / T; @@ -272,18 +288,22 @@ namespace thermalfist { double QuantumNumericalIntegrationEntropyDensity(int statistics, double T, double mu, double m, double deg) { + if (T == 0.) + return 0.; return (QuantumNumericalIntegrationPressure(statistics, T, mu, m, deg) + QuantumNumericalIntegrationEnergyDensity(statistics, T, mu, m, deg) - mu * QuantumNumericalIntegrationDensity(statistics, T, mu, m, deg)) / T; } double QuantumNumericalIntegrationScalarDensity(int statistics, double T, double mu, double m, double deg) { if (statistics == 0) return BoltzmannScalarDensity(T, mu, m, deg); + if (statistics == 1 && T == 0.) return FermiZeroTScalarDensity(mu, m, deg); if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuScalarDensity(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationScalarDensity: Bose-Einstein condensation\n"); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; double ret = 0.; double moverT = m / T; @@ -300,13 +320,15 @@ namespace thermalfist { double QuantumNumericalIntegrationT1dn1dmu1(int statistics, double T, double mu, double m, double deg) { - if (statistics == 0) return BoltzmannTdndmu(1, T, mu, m, deg); - if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuT1dn1dmu1(T, mu, m, deg); + if (statistics == 0) return BoltzmannTdndmu(1, T, mu, m, deg); + if (statistics == 1 && T == 0.) return 0.; + if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuT1dn1dmu1(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationT1dn1dmu1: Bose-Einstein condensation\n"); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; double ret = 0.; double moverT = m / T; @@ -324,13 +346,15 @@ namespace thermalfist { double QuantumNumericalIntegrationT2dn2dmu2(int statistics, double T, double mu, double m, double deg) { - if (statistics == 0) return BoltzmannTdndmu(2, T, mu, m, deg); - if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuT2dn2dmu2(T, mu, m, deg); + if (statistics == 0) return BoltzmannTdndmu(2, T, mu, m, deg); + if (statistics == 1 && T == 0.) return 0.; + if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuT2dn2dmu2(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationT2dn2dmu2: Bose-Einstein condensation\n"); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; double ret = 0.; double moverT = m / T; @@ -349,13 +373,15 @@ namespace thermalfist { double QuantumNumericalIntegrationT3dn3dmu3(int statistics, double T, double mu, double m, double deg) { - if (statistics == 0) return BoltzmannTdndmu(3, T, mu, m, deg); - if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuT3dn3dmu3(T, mu, m, deg); + if (statistics == 0) return BoltzmannTdndmu(3, T, mu, m, deg); + if (statistics == 1 && T == 0.) return 0.; + if (statistics == 1 && mu > m) return FermiNumericalIntegrationLargeMuT3dn3dmu3(T, mu, m, deg); if (statistics == -1 && mu > m) { printf("**WARNING** QuantumNumericalIntegrationT3dn3dmu3: Bose-Einstein condensation\n"); calculationHadBECIssue = true; return 0.; } + if (statistics == -1 && T == 0.) return 0.; //printf("\n"); @@ -400,6 +426,20 @@ namespace thermalfist { return QuantumNumericalIntegrationTdndmu(N - 1, statistics, T, mu, m, deg) / pow(T, 3) / xMath::GeVtoifm3(); } + double QuantumNumericalIntegrationChiNDimensionfull(int N, int statistics, double T, double mu, double m, double deg) + { + if (statistics == 1 && T == 0.0) + return FermiZeroTChiNDimensionfull(N, mu, m, deg); + if (statistics == -1 && T == 0.0) { + if (mu >= m) { + printf("**WARNING** QuantumNumericalIntegrationChiNDimensionfull: Bose-Einstein condensation\n"); + calculationHadBECIssue = true; + } + return 0.; + } + return QuantumNumericalIntegrationTdndmu(N - 1, statistics, T, mu, m, deg) / pow(T, N-1) / xMath::GeVtoifm3(); + } + double psi(double x) { if (x == 0.0) @@ -637,6 +677,127 @@ namespace thermalfist { return FermiNumericalIntegrationLargeMuTdndmu(N - 1, T, mu, m, deg) / pow(T, 3) / xMath::GeVtoifm3(); } + double FermiNumericalIntegrationLargeMuChiNDimensionfull(int N, double T, double mu, double m, double deg) + { + return FermiNumericalIntegrationLargeMuTdndmu(N - 1, T, mu, m, deg) / pow(T, N-1) / xMath::GeVtoifm3(); + } + + double FermiZeroTDensity(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + return deg / 6. / xMath::Pi() / xMath::Pi() * pf * pf * pf * xMath::GeVtoifm3(); + } + + double FermiZeroTPressure(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + if (m == 0.0) { + return deg / 24. / xMath::Pi() / xMath::Pi() * pf * pf * pf * pf * xMath::GeVtoifm3(); + } + double m2 = m * m; + return deg / 48. / xMath::Pi() / xMath::Pi() * + ( + mu * pf * (2. * mu * mu - 5. * m2) + - 3. * m2 * m2 * log(m / (mu + pf)) + ) * xMath::GeVtoifm3(); + } + + double FermiZeroTEnergyDensity(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + if (m == 0.0) { + return deg / 8. / xMath::Pi() / xMath::Pi() * pf * pf * pf * pf * xMath::GeVtoifm3(); + } + double m2 = m * m; + return deg / 16. / xMath::Pi() / xMath::Pi() * + ( + mu * pf * (2. * mu * mu - m2) + + m2 * m2 * log(m / (mu + pf)) + ) * xMath::GeVtoifm3(); + } + + double FermiZeroTEntropyDensity(double mu, double m, double deg) + { + return 0.0; + } + + double FermiZeroTScalarDensity(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + if (m == 0.0) { + return 0.; + } + double m2 = m * m; + return deg * m / 4. / xMath::Pi() / xMath::Pi() * + ( + mu * pf + + m2 * log(m / (mu + pf)) + ) * xMath::GeVtoifm3(); + } + + double FermiZeroTdn1dmu1(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + return deg / 2. / xMath::Pi() / xMath::Pi() * mu * pf * xMath::GeVtoifm3(); + } + + double FermiZeroTdn2dmu2(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + return deg / 2. / xMath::Pi() / xMath::Pi() * (mu * mu + pf * pf) / pf * xMath::GeVtoifm3(); + } + + double FermiZeroTdn3dmu3(double mu, double m, double deg) + { + if (m >= mu) + return 0.0; + double pf = sqrt(mu * mu - m * m); + return deg / 2. / xMath::Pi() / xMath::Pi() * mu * (3. * pf * pf - mu * mu) / pf / pf / pf * xMath::GeVtoifm3(); + } + + double FermiZeroTdndmu(int N, double mu, double m, double deg) + { + if (N < 0 || N>3) { + printf("**ERROR** FermiNumericalIntegrationLargeMuTdndmu: N < 0 or N > 3\n"); + exit(1); + } + if (N == 0) + return FermiZeroTDensity(mu, m, deg); + + if (N == 1) + return FermiZeroTdn1dmu1(mu, m, deg); + + if (N == 2) + return FermiZeroTdn2dmu2(mu, m, deg); + + return FermiZeroTdn3dmu3(mu, m, deg); + } + + double FermiZeroTChiN(int N, double mu, double m, double deg) + { + printf("**ERROR** FermiZeroTChiN: This quantity is infinite by definition at T = 0!\n"); + exit(1); + return 0.0; + //return FermiNumericalIntegrationLargeMuTdndmu(N - 1, mu, m, deg) / pow(T, 3) / xMath::GeVtoifm3(); + } + + double FermiZeroTChiNDimensionfull(int N, double mu, double m, double deg) + { + return FermiZeroTdndmu(N - 1, mu, m, deg) / xMath::GeVtoifm3(); + } + double IdealGasQuantity(Quantity quantity, QStatsCalculationType calctype, int statistics, double T, double mu, double m, double deg, int order) { if (statistics == 0) { @@ -656,6 +817,12 @@ namespace thermalfist { return BoltzmannChiN(3, T, mu, m, deg); if (quantity == chi4) return BoltzmannChiN(4, T, mu, m, deg); + if (quantity == chi2difull) + return BoltzmannChiNDimensionfull(2, T, mu, m, deg); + if (quantity == chi3difull) + return BoltzmannChiNDimensionfull(3, T, mu, m, deg); + if (quantity == chi4difull) + return BoltzmannChiNDimensionfull(4, T, mu, m, deg); } else { if (calctype == ClusterExpansion) { @@ -675,6 +842,12 @@ namespace thermalfist { return QuantumClusterExpansionChiN(3, statistics, T, mu, m, deg, order); if (quantity == chi4) return QuantumClusterExpansionChiN(4, statistics, T, mu, m, deg, order); + if (quantity == chi2difull) + return QuantumClusterExpansionChiNDimensionfull(2, statistics, T, mu, m, deg, order); + if (quantity == chi3difull) + return QuantumClusterExpansionChiNDimensionfull(3, statistics, T, mu, m, deg, order); + if (quantity == chi4difull) + return QuantumClusterExpansionChiNDimensionfull(4, statistics, T, mu, m, deg, order); } else { if (quantity == ParticleDensity) @@ -693,6 +866,12 @@ namespace thermalfist { return QuantumNumericalIntegrationChiN(3, statistics, T, mu, m, deg); if (quantity == chi4) return QuantumNumericalIntegrationChiN(4, statistics, T, mu, m, deg); + if (quantity == chi2difull) + return QuantumNumericalIntegrationChiNDimensionfull(2, statistics, T, mu, m, deg); + if (quantity == chi3difull) + return QuantumNumericalIntegrationChiNDimensionfull(3, statistics, T, mu, m, deg); + if (quantity == chi4difull) + return QuantumNumericalIntegrationChiNDimensionfull(4, statistics, T, mu, m, deg); } } printf("**WARNING** IdealGasFunctions::IdealGasQuantity: Unknown quantity\n"); diff --git a/src/library/HRGBase/ThermalModelBase.cpp b/src/library/HRGBase/ThermalModelBase.cpp index e0fcd64..3f28ad8 100644 --- a/src/library/HRGBase/ThermalModelBase.cpp +++ b/src/library/HRGBase/ThermalModelBase.cpp @@ -374,6 +374,15 @@ namespace thermalfist { return m_Chem[i]; } + void ThermalModelBase::SetChemicalPotential(int i, double chem) + { + if (i < 0 || i >= static_cast(m_Chem.size())) { + printf("**ERROR** ThermalModelBase::SetChemicalPotential(int i): i is out of bounds!"); + exit(1); + } + m_Chem[i] = chem; + } + double ThermalModelBase::FullIdealChemicalPotential(int i) const { if (i < 0 || i >= static_cast(m_Chem.size())) { @@ -1405,10 +1414,12 @@ namespace thermalfist { fACd = m_THM->CalculateAbsoluteCharmDensityModulo(); } - vector wprim; - wprim.resize(m_THM->Densities().size()); - for (size_t i = 0; i < wprim.size(); ++i) - wprim[i] = m_THM->ParticleScaledVariance(i); + vector chi2s; + chi2s.resize(m_THM->Densities().size()); + for (size_t i = 0; i < chi2s.size(); ++i) { + chi2s[i] = m_THM->TPS()->Particle(i).chiDimensionfull(2, m_THM->Parameters(), m_THM->UseWidth(), m_THM->ChemicalPotential(i) + m_THM->MuShift(i)) + * xMath::GeVtoifm3(); + } int NNN = 0; if (m_THM->ConstrainMuQ()) NNN++; @@ -1425,14 +1436,12 @@ namespace thermalfist { if (m_THM->ConstrainMuQ()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i]; // Update: remove division by Q/B to allow for charge neutrality ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2); @@ -1444,14 +1453,12 @@ namespace thermalfist { if (m_THM->ConstrainMuS()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).Strangeness() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).Strangeness() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i]; // Update: remove division by Q/B to allow for charge neutrality ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2); @@ -1463,14 +1470,12 @@ namespace thermalfist { if (m_THM->ConstrainMuC()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).Charm() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).Charm() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).Charm() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).Charm() * chi2s[i]; // Update: remove division by Q/B to allow for charge neutrality ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2); @@ -1491,14 +1496,12 @@ namespace thermalfist { if (m_THM->ConstrainMuQ()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i]; ret(i1, i2) = d1 / fASd - fSd / fASd / fASd * d2; @@ -1508,14 +1511,12 @@ namespace thermalfist { if (m_THM->ConstrainMuS()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).Strangeness() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).Strangeness() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i]; ret(i1, i2) = d1 / fASd - fSd / fASd / fASd * d2; @@ -1525,14 +1526,12 @@ namespace thermalfist { if (m_THM->ConstrainMuC()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).Charm() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).Charm() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).Charm() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).Charm() * chi2s[i]; ret(i1, i2) = d1 / fASd - fSd / fASd / fASd * d2; @@ -1551,14 +1550,12 @@ namespace thermalfist { if (m_THM->ConstrainMuQ()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).ElectricCharge() *chi2s[i]; ret(i1, i2) = d1 / fACd - fCd / fACd / fACd * d2; @@ -1568,14 +1565,12 @@ namespace thermalfist { if (m_THM->ConstrainMuS()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).Strangeness() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).Strangeness() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i]; ret(i1, i2) = d1 / fACd - fCd / fACd / fACd * d2; @@ -1585,14 +1580,12 @@ namespace thermalfist { if (m_THM->ConstrainMuC()) { d1 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).Charm() * m_THM->Densities()[i] * wprim[i]; - d1 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).Charm() * chi2s[i]; d2 = 0.; - for (size_t i = 0; i < wprim.size(); ++i) - d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).Charm() * m_THM->Densities()[i] * wprim[i]; - d2 /= m_THM->Parameters().T; + for (size_t i = 0; i < chi2s.size(); ++i) + d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).Charm() * chi2s[i]; ret(i1, i2) = d1 / fACd - fCd / fACd / fACd * d2; @@ -1859,12 +1852,6 @@ namespace thermalfist { } } - vector tfug(4, 0.); - tfug[0] = exp(m_THM->Parameters().muB / m_THM->Parameters().T); - tfug[1] = exp(m_THM->Parameters().muQ / m_THM->Parameters().T); - tfug[2] = exp(m_THM->Parameters().muS / m_THM->Parameters().T); - tfug[3] = exp(m_THM->Parameters().muC / m_THM->Parameters().T); - m_THM->FillChemicalPotentials(); m_THM->CalculatePrimordialDensities(); @@ -1892,22 +1879,23 @@ namespace thermalfist { } } - vector wprim; - wprim.resize(m_THM->Densities().size()); - for (size_t i = 0; i < wprim.size(); ++i) wprim[i] = m_THM->ParticleScaledVariance(i); + vector chi2s; + chi2s.resize(m_THM->Densities().size()); + for (size_t i = 0; i < chi2s.size(); ++i) { + chi2s[i] = m_THM->TPS()->Particle(i).chiDimensionfull(2, m_THM->Parameters(), m_THM->UseWidth(), m_THM->ChemicalPotential(i) + m_THM->MuShift(i)) + * xMath::GeVtoifm3(); + } vector< vector > deriv(4, vector(4)), derivabs(4, vector(4)); for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) { deriv[i][j] = 0.; - for (size_t part = 0; part < wprim.size(); ++part) - deriv[i][j] += m_THM->TPS()->Particles()[part].GetCharge(i) * m_THM->TPS()->Particles()[part].GetCharge(j) * m_THM->Densities()[part] * wprim[part]; - deriv[i][j] /= m_THM->Parameters().T; + for (size_t part = 0; part < chi2s.size(); ++part) + deriv[i][j] += m_THM->TPS()->Particles()[part].GetCharge(i) * m_THM->TPS()->Particles()[part].GetCharge(j) * chi2s[part]; derivabs[i][j] = 0.; - for (size_t part = 0; part < wprim.size(); ++part) - derivabs[i][j] += m_THM->TPS()->Particles()[part].GetAbsCharge(i) * m_THM->TPS()->Particles()[part].GetCharge(j) * m_THM->Densities()[part] * wprim[part]; - derivabs[i][j] /= m_THM->Parameters().T; + for (size_t part = 0; part < chi2s.size(); ++part) + derivabs[i][j] += m_THM->TPS()->Particles()[part].GetAbsCharge(i) * m_THM->TPS()->Particles()[part].GetCharge(j) * chi2s[part]; } diff --git a/src/library/HRGBase/ThermalModelIdeal.cpp b/src/library/HRGBase/ThermalModelIdeal.cpp index 38b8c10..00d6b3e 100644 --- a/src/library/HRGBase/ThermalModelIdeal.cpp +++ b/src/library/HRGBase/ThermalModelIdeal.cpp @@ -45,10 +45,13 @@ namespace thermalfist { void ThermalModelIdeal::CalculateTwoParticleCorrelations() { int NN = m_densities.size(); - vector tN(NN), tW(NN); - for (int i = 0; i < NN; ++i) tN[i] = m_densities[i]; + vector tN(NN); + for (int i = 0; i < NN; ++i) + tN[i] = m_densities[i]; - for (int i = 0; i < NN; ++i) tW[i] = ParticleScaledVariance(i); + vector chi2s(NN); + for (int i = 0; i < NN; ++i) + chi2s[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_Chem[i]); m_PrimCorrel.resize(NN); for (int i = 0; i < NN; ++i) m_PrimCorrel[i].resize(NN); @@ -57,8 +60,7 @@ namespace thermalfist { for (int i = 0; i < NN; ++i) for (int j = 0; j < NN; ++j) { m_PrimCorrel[i][j] = 0.; - if (i == j) m_PrimCorrel[i][j] += m_densities[i] * tW[i]; - m_PrimCorrel[i][j] /= m_Parameters.T; + if (i == j) m_PrimCorrel[i][j] += m_densities[i] * chi2s[i]; } for (int i = 0; i < NN; ++i) { diff --git a/src/library/HRGBase/ThermalParticle.cpp b/src/library/HRGBase/ThermalParticle.cpp index b453c70..48e2226 100644 --- a/src/library/HRGBase/ThermalParticle.cpp +++ b/src/library/HRGBase/ThermalParticle.cpp @@ -707,12 +707,22 @@ namespace thermalfist { } + double ThermalParticle::chiDimensionfull(int index, const ThermalModelParameters& params, bool useWidth, double mu) const + { + if (index == 0) return Density(params, IdealGasFunctions::Pressure, useWidth, mu) / pow(xMath::GeVtoifm(), 3); + if (index == 1) return Density(params, IdealGasFunctions::ParticleDensity, useWidth, mu) / pow(xMath::GeVtoifm(), 3); + if (index == 2) return Density(params, IdealGasFunctions::chi2difull, useWidth, mu); + if (index == 3) return Density(params, IdealGasFunctions::chi3difull, useWidth, mu); + if (index == 4) return Density(params, IdealGasFunctions::chi4difull, useWidth, mu); + return 1.; + } + double ThermalParticle::ScaledVariance(const ThermalModelParameters ¶ms, bool useWidth, double mu) const { if (m_Degeneracy == 0.0) return 1.; if (m_Statistics == 0) return 1.; double dens = Density(params, IdealGasFunctions::ParticleDensity, useWidth, mu); if (dens == 0.) return 1.; - double ret = chi(2, params, useWidth, mu) / chi(1, params, useWidth, mu); + double ret = params.T * chiDimensionfull(2, params, useWidth, mu) / chiDimensionfull(1, params, useWidth, mu); if (ret != ret) ret = 1.; return ret; } @@ -723,7 +733,7 @@ namespace thermalfist { if (m_Statistics == 0) return 1.; double dens = Density(params, IdealGasFunctions::ParticleDensity, useWidth, mu); if (dens == 0.) return 1.; - double ret = chi(3, params, useWidth, mu) / chi(2, params, useWidth, mu); + double ret = params.T * chiDimensionfull(3, params, useWidth, mu) / chiDimensionfull(2, params, useWidth, mu); if (ret != ret) ret = 1.; return ret; } @@ -734,7 +744,7 @@ namespace thermalfist { if (m_Statistics == 0) return 1.; double dens = Density(params, IdealGasFunctions::ParticleDensity, useWidth, mu); if (dens == 0.) return 1.; - double ret = chi(4, params, useWidth, mu) / chi(2, params, useWidth, mu); + double ret = params.T * params.T * chiDimensionfull(4, params, useWidth, mu) / chiDimensionfull(2, params, useWidth, mu); if (ret != ret) ret = 1.; return ret; } diff --git a/src/library/HRGEV/ThermalModelEVCrosstermsLegacy.cpp b/src/library/HRGEV/ThermalModelEVCrosstermsLegacy.cpp index 4e9441d..a783ffb 100644 --- a/src/library/HRGEV/ThermalModelEVCrosstermsLegacy.cpp +++ b/src/library/HRGEV/ThermalModelEVCrosstermsLegacy.cpp @@ -489,9 +489,16 @@ namespace thermalfist { void ThermalModelEVCrosstermsLegacy::CalculateTwoParticleCorrelations() { int NN = m_densities.size(); - vector tN(NN), tW(NN); + vector tN(NN);// , tW(NN); for (int i = 0; i < NN; ++i) tN[i] = DensityId(i); - for (int i = 0; i < NN; ++i) tW[i] = ScaledVarianceId(i); + //for (int i = 0; i < NN; ++i) tW[i] = ScaledVarianceId(i); + vector chi2id(NN); + for (int i = 0; i < NN; ++i) { + double dMu = 0.; + for (int j = 0; j < m_TPS->ComponentsNumber(); ++j) dMu += -m_Virial[i][j] * m_Ps[j]; + chi2id[i] = m_densities[i] / tN[i] * m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_Chem[i] + dMu) * xMath::GeVtoifm3(); + } + MatrixXd densMatrix(NN, NN); VectorXd solVector(NN), xVector(NN), xVector2(NN); @@ -545,7 +552,8 @@ namespace thermalfist { for (int i = 0; i < NN; ++i) for (int j = i; j < NN; ++j) { for (int l = 0; l < NN; ++l) - xVector[l] = tN[l] / m_Parameters.T * tW[l] * coefs[l][i] * coefs[l][j]; + //xVector[l] = tN[l] / m_Parameters.T * tW[l] * coefs[l][i] * coefs[l][j]; + xVector[l] = chi2id[l] * coefs[l][i] * coefs[l][j]; solVector = decomp.solve(xVector); if (1) { //if (decomp.info()==Eigen::Success) { diff --git a/src/library/HRGEV/ThermalModelEVDiagonal.cpp b/src/library/HRGEV/ThermalModelEVDiagonal.cpp index f50e0c1..d4be9d2 100644 --- a/src/library/HRGEV/ThermalModelEVDiagonal.cpp +++ b/src/library/HRGEV/ThermalModelEVDiagonal.cpp @@ -170,6 +170,8 @@ namespace thermalfist { m_Pressure = Pressure(0.); double mnc = pow(m_Parameters.T, 4.) * pow(xMath::GeVtoifm(), 3.); + if (m_Parameters.T == 0.0) + mnc = pow(xMath::mnucleon(), 4.) * pow(xMath::GeVtoifm(), 3.); eqs.SetMnc(mnc); jac.SetMnc(mnc); double x0 = log(m_Pressure / mnc); @@ -241,26 +243,32 @@ namespace thermalfist { void ThermalModelEVDiagonal::CalculateTwoParticleCorrelations() { int NN = m_densities.size(); - vector tN(NN), tW(NN); + vector tN(NN); for (int i = 0; i < NN; ++i) tN[i] = DensityId(i, m_Pressure); - for (int i = 0; i < NN; ++i) tW[i] = ScaledVarianceId(i, m_Pressure); + + vector chi2id(NN); + for (int i = 0; i < NN; ++i) { + chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure) * xMath::GeVtoifm3(); + if (tN[i] > 0.0) + chi2id[i] *= m_densities[i] / tN[i]; + } m_PrimCorrel.resize(NN); for (int i = 0; i < NN; ++i) m_PrimCorrel[i].resize(NN); m_TotalCorrel = m_PrimCorrel; - for (int i = 0; i < NN; ++i) + for (int i = 0; i < NN; ++i) { for (int j = 0; j < NN; ++j) { m_PrimCorrel[i][j] = 0.; - if (i == j) m_PrimCorrel[i][j] += m_densities[i] * tW[i]; - m_PrimCorrel[i][j] += -m_v[i] * m_densities[i] * m_densities[j] * tW[i]; - m_PrimCorrel[i][j] += -m_v[j] * m_densities[i] * m_densities[j] * tW[j]; + if (i == j) m_PrimCorrel[i][j] += chi2id[i]; + m_PrimCorrel[i][j] += -m_v[i] * m_densities[j] * chi2id[i]; + m_PrimCorrel[i][j] += -m_v[j] * m_densities[i] * chi2id[j]; double tmp = 0.; for (size_t k = 0; k < m_densities.size(); ++k) - tmp += m_v[k] * m_v[k] * tW[k] * m_densities[k]; + tmp += m_v[k] * m_v[k] * chi2id[k]; m_PrimCorrel[i][j] += m_densities[i] * m_densities[j] * tmp; - m_PrimCorrel[i][j] /= m_Parameters.T; } + } for (int i = 0; i < NN; ++i) { m_wprim[i] = m_PrimCorrel[i][i]; diff --git a/src/library/HRGVDW/ThermalModelVDW.cpp b/src/library/HRGVDW/ThermalModelVDW.cpp index 448125e..4c6060f 100644 --- a/src/library/HRGVDW/ThermalModelVDW.cpp +++ b/src/library/HRGVDW/ThermalModelVDW.cpp @@ -795,7 +795,8 @@ namespace thermalfist { vector chi2id(m_densities.size()); for (int i = 0; iParticles()[i].chi(2, m_Parameters, m_UseWidth, m_MuStar[i]); + //chi2id[i] = m_TPS->Particles()[i].chi(2, m_Parameters, m_UseWidth, m_MuStar[i]); + chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_MuStar[i]); for (int i = 0; i chi2s(NN, 0.); for (int i = 0; iTPS()->Particles()[i].chi(2, m_THM->Parameters(), + chi2s[i] = m_THM->TPS()->Particles()[i].chiDimensionfull(2, m_THM->Parameters(), m_THM->UseWidth(), m_THM->ChemicalPotential(i) + x[m_THM->m_MapTodMuStar[i]] ); @@ -1102,7 +1104,7 @@ namespace thermalfist { double tmps = 1.; if (ns[ti] != 0.) tmps = np[ti] / ns[ti]; - xVector[l] += chi2s[ti] * m_THM->m_Parameters.T * m_THM->m_Parameters.T * pow(xMath::GeVtoifm(), 3) * tmps; + xVector[l] += chi2s[ti] * pow(xMath::GeVtoifm(), 3) * tmps; } } @@ -1122,7 +1124,7 @@ namespace thermalfist { double tmps = 1.; if (ns[j] != 0.) tmps = np[j] / ns[j]; - dnjdmukp[j] += tmps * chi2s[j] * m_THM->m_Parameters.T * m_THM->m_Parameters.T * pow(xMath::GeVtoifm(), 3); + dnjdmukp[j] += tmps * chi2s[j] * pow(xMath::GeVtoifm(), 3); } } }