diff --git a/src/FSSW.cpp b/src/FSSW.cpp index 5f49e22..e049d22 100644 --- a/src/FSSW.cpp +++ b/src/FSSW.cpp @@ -603,22 +603,23 @@ void FSSW::calculate_dN_dxtdy_for_one_particle_species( get22momNEOSBQSCoefficients(surf->Edec, surf->Bn, visCoefficients); } + if (deltaf_kind_ == 31) { + getNewRTACoefficients(temp, visCoefficients); + } + // bulk delta f contribution double bulkPi = 0.0; if (INCLUDE_BULK_DELTAF) { - if (deltaf_kind_ == 21 || deltaf_kind_ == 20) { - bulkPi = surf->bulkPi; // GeV/fm^3 - } else if (deltaf_kind_ == 11) { - bulkPi = surf->bulkPi; // GeV/fm^3 - getbulkvisCoefficients(temp, mu_B, visCoefficients); - } else { - if (deltaf_kind_ == 0) { - bulkPi = surf->bulkPi; // unit in GeV/fm^3 - } else { - bulkPi = surf->bulkPi / hbarC; // unit in fm^-4 - } + bulkPi = surf->bulkPi; // GeV/fm^3 + if (deltaf_kind_ > 0 && deltaf_kind_ < 10) { + bulkPi = surf->bulkPi / hbarC; // unit in fm^-4 + } + if (deltaf_kind_ < 10) { getbulkvisCoefficients(temp, visCoefficients); } + if (deltaf_kind_ == 11) { + getbulkvisCoefficients(temp, mu_B, visCoefficients); + } } // diffusion delta f @@ -695,6 +696,12 @@ void FSSW::calculate_dN_dxtdy_for_one_particle_species( + charge * visCoefficients[5]) + results_ptr[3] * (visCoefficients[1] - visCoefficients[2]))); + } else if (deltaf_kind_ == 31) { + deltaN_bulk = + (unit_factor * prefactor * dsigma_dot_u * bulkPi + / visCoefficients[1] + * (visCoefficients[2] * results_ptr[3] + results_ptr[1] + + visCoefficients[3] * results_ptr[2])); } } @@ -831,6 +838,19 @@ void FSSW::calculate_dN_analytic( deltaN_bulk_term2 = mass * mass / (beta * beta) * deltaN_bulk_term2; deltaN_bulk_term3 = mass * mass * mass / (beta * beta) * deltaN_bulk_term3; + } else if (deltaf_kind_ == 31) { + std::vector thermoCoeffs; + if (sign == -1) { + getNewRTAThermoCoefficients(Temperature, thermoCoeffs, true); + } else { + getNewRTAThermoCoefficients(Temperature, thermoCoeffs, false); + } + deltaN_bulk_term1 = + -Temperature * mass * mass / 3. * thermoCoeffs[0]; + deltaN_bulk_term2 = + Temperature * Temperature * Temperature * thermoCoeffs[1]; + deltaN_bulk_term3 = + Temperature * Temperature * Temperature * thermoCoeffs[2]; } } else { deltaN_bulk_term1 = 0.0; @@ -1478,6 +1498,59 @@ void FSSW::load_deltaf_table_newRTA(std::string filepath) { bulkTable_cs2 >> dummy >> deltaf_coeff_newRTA_bulk_cs2_[j]; } bulkTable_cs2.close(); + + sf_J_dgamma_ = deltaf_newRTA_dgamma_; + sf_J_gammaMin_ = deltaf_newRTA_gamma0_; + sf_J_dm_ = -1.; + sf_J_mMin_ = 0.; + + std::string thermFileName = filepath + "/thermoIntegral.dat/J_gamma_BE.dat"; + loadThermoIntegrals(thermFileName, 4, sf_logJgamma_BE_); + thermFileName = filepath + "/thermoIntegral.dat/J_gamma_FD.dat"; + loadThermoIntegrals(thermFileName, 4, sf_logJgamma_FD_); + thermFileName = filepath + "/thermoIntegral.dat/J_gammaplus2_BE.dat"; + loadThermoIntegrals(thermFileName, 4, sf_logJgammaplus2_BE_); + thermFileName = filepath + "/thermoIntegral.dat/J_gammaplus2_FD.dat"; + loadThermoIntegrals(thermFileName, 4, sf_logJgammaplus2_FD_); + thermFileName = filepath + "/thermoIntegral.dat/J_2_BE.dat"; + loadThermoIntegrals(thermFileName, 1, sf_logJ2_BE_); + thermFileName = filepath + "/thermoIntegral.dat/J_2_FD.dat"; + loadThermoIntegrals(thermFileName, 1, sf_logJ2_FD_); +} + +void FSSW::loadThermoIntegrals( + std::string thermoFileName, const int ncol, + std::vector> &data) { + data.clear(); + // read in thermoIntegral tables + ifstream thermoTable(thermoFileName.c_str()); + if (!thermoTable) { + messager_ << "Can not found file: " << thermoFileName; + messager_.flush("error"); + exit(1); + } + double dummy_d; + std::string dummy; + // read header + std::getline(thermoTable, dummy); + std::vector temp(ncol, 0.); + thermoTable >> sf_J_mMin_; + for (int i = 0; i < ncol; i++) { + thermoTable >> temp[i]; + } + while (!thermoTable.eof()) { + for (unsigned int i = 0; i < temp.size(); i++) { + temp[i] = log(temp[i]); + } + data.push_back(temp); + thermoTable >> dummy_d; + for (int i = 0; i < ncol; i++) { + thermoTable >> temp[i]; + } + if (sf_J_dm_ < 0.) { + sf_J_dm_ = dummy_d - sf_J_mMin_; + } + } } void FSSW::getbulkvisCoefficients( @@ -1700,6 +1773,67 @@ void FSSW::getNewRTACoefficients( + deltaf_coeff_newRTA_bulk_cs2_[idx_T + 1] * frac_T; } +void FSSW::getNewRTAThermoCoefficients( + const double m, std::vector &thermoCoefficients, + const bool bosonFlag) { + thermoCoefficients.resize(3, 0.); // J_gamma, J_gamma+2, J_2 + int idx_m = static_cast((m - sf_J_mMin_) / sf_J_dm_); + const int mtableLength = sf_logJgamma_BE_.size(); + idx_m = std::max(0, std::min(mtableLength - 2, idx_m)); + double frac_m = (m - sf_J_mMin_) / sf_J_dm_ - idx_m; + frac_m = std::max(0., frac_m); // allow frac_m > 1 for extrapolation + + int idx_gamma = static_cast( + (deltaf_newRTA_gamma_ - sf_J_gammaMin_) / sf_J_dgamma_); + const int gammatableLength = sf_logJgamma_BE_[0].size(); + idx_gamma = std::max(0, std::min(gammatableLength - 2, idx_gamma)); + double frac_gamma = + (deltaf_newRTA_gamma_ - sf_J_gammaMin_) / sf_J_dgamma_ - idx_gamma; + frac_gamma = std::max(0., std::min(1., frac_gamma)); + + double temp1, temp2; + // J_gamma + if (bosonFlag) { + temp1 = sf_logJgamma_BE_[idx_m][idx_gamma] * (1. - frac_m) + + sf_logJgamma_BE_[idx_m + 1][idx_gamma] * frac_m; + temp2 = sf_logJgamma_BE_[idx_m][idx_gamma + 1] * (1. - frac_m) + + sf_logJgamma_BE_[idx_m + 1][idx_gamma + 1] * frac_m; + } else { + temp1 = sf_logJgamma_FD_[idx_m][idx_gamma] * (1. - frac_m) + + sf_logJgamma_FD_[idx_m + 1][idx_gamma] * frac_m; + temp2 = sf_logJgamma_FD_[idx_m][idx_gamma + 1] * (1. - frac_m) + + sf_logJgamma_FD_[idx_m + 1][idx_gamma + 1] * frac_m; + } + thermoCoefficients[0] = temp1 * (1. - frac_gamma) + temp2 * frac_gamma; + + // J_gamma + 2 + if (bosonFlag) { + temp1 = sf_logJgammaplus2_BE_[idx_m][idx_gamma] * (1. - frac_m) + + sf_logJgammaplus2_BE_[idx_m + 1][idx_gamma] * frac_m; + temp2 = sf_logJgammaplus2_BE_[idx_m][idx_gamma + 1] * (1. - frac_m) + + sf_logJgammaplus2_BE_[idx_m + 1][idx_gamma + 1] * frac_m; + } else { + temp1 = sf_logJgammaplus2_FD_[idx_m][idx_gamma] * (1. - frac_m) + + sf_logJgammaplus2_FD_[idx_m + 1][idx_gamma] * frac_m; + temp2 = sf_logJgammaplus2_FD_[idx_m][idx_gamma + 1] * (1. - frac_m) + + sf_logJgammaplus2_FD_[idx_m + 1][idx_gamma + 1] * frac_m; + } + thermoCoefficients[1] = temp1 * (1. - frac_gamma) + temp2 * frac_gamma; + + // J_2 + if (bosonFlag) { + temp1 = sf_logJ2_BE_[idx_m][0] * (1. - frac_m) + + sf_logJ2_BE_[idx_m + 1][0] * frac_m; + } else { + temp1 = sf_logJ2_FD_[idx_m][0] * (1. - frac_m) + + sf_logJ2_FD_[idx_m + 1][0] * frac_m; + } + thermoCoefficients[2] = temp1; + for (auto &coeff : thermoCoefficients) { + coeff = std::exp(coeff); + } +} + void FSSW::load_deltaf_qmu_coeff_table(string filename) { ifstream table(filename.c_str()); deltaf_qmu_coeff_table_length_T = 150; diff --git a/src/FSSW.h b/src/FSSW.h index 8eae356..47c09dc 100644 --- a/src/FSSW.h +++ b/src/FSSW.h @@ -118,6 +118,14 @@ class FSSW { int sf_expint_truncate_order; double **sf_bessel_Kn; double **sf_expint_En; + double sf_J_mMin_, sf_J_dm_; + double sf_J_gammaMin_, sf_J_dgamma_; + std::vector> sf_logJgamma_BE_; + std::vector> sf_logJgamma_FD_; + std::vector> sf_logJgammaplus2_BE_; + std::vector> sf_logJgammaplus2_FD_; + std::vector> sf_logJ2_BE_; + std::vector> sf_logJ2_FD_; int flag_output_samples_into_files; int flag_store_samples_in_memory; @@ -178,6 +186,10 @@ class FSSW { void getNewRTACoefficients( const double Tdec, std::vector &visCoefficients); + void getNewRTAThermoCoefficients( + const double m, std::vector &thermoCoefficients, + const bool bosonFlag); + void load_bulk_deltaf_14mom_table(string filepath); void load_CE_deltaf_NEOSBQS_table(string filepath); void load_22mom_deltaf_NEOSBQS_table(string filepath); @@ -185,6 +197,9 @@ class FSSW { double get_deltaf_qmu_coeff(double T, double muB); void load_deltaf_table_newRTA(std::string filepath); + void loadThermoIntegrals( + std::string thermoFileName, const int ncol, + std::vector> &data); void check_samples_in_memory(); int get_number_of_sampled_events() { return (Hadron_list->size()); }; diff --git a/utilities/computeThermodynamicJ.py b/utilities/computeThermodynamicJ.py new file mode 100644 index 0000000..dce1c26 --- /dev/null +++ b/utilities/computeThermodynamicJ.py @@ -0,0 +1,311 @@ +#!/urs/bin/env python3 + +""" + This script compute thermodynamic integrals with special functions +""" + +import mpmath +import numpy as np + +mpmath.mp.dpf = 50 +mpmath.mp.pretty = True + + +def computeThermodynamicJ(nu: float, m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_nu(m) + Eq. (9) in https://www.overleaf.com/project/674f23fdff4ede2888cab91e + """ + res = 0. + if m < 20.: + truncOrder = 5 + else: + truncOrder = 1 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder + 1): + arg = (k * m)**2 / 4. + res += sign**(k-1) * k * ( + 1./(k**(nu + 2)) * mpmath.gamma(nu + 2) + * mpmath.hyp1f2(-0.5, -0.5 - nu/2., -nu/2., arg) + + m**(nu + 2)/4.*np.sqrt(np.pi)*( + mpmath.gamma(-1 - nu/2.) / mpmath.gamma(0.5 - nu/2.) + * mpmath.hyp1f2(0.5 + nu/2., 0.5, 2 + nu/2., arg) + - k * m * mpmath.gamma(-1.5 - nu/2.) / mpmath.gamma(- nu/2.) + * mpmath.hyp1f2(1 + nu/2., 1.5, 2.5 + nu/2., arg)) + ) + return res + + +def thermodynamicJ0(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_0(m) + Eq. (10) in https://www.overleaf.com/project/674f23fdff4ede2888cab91e + """ + res = 0. + truncOrder = 5 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) + res += sign**(k-1) * m * mpmath.besselk(1, arg) + return res + + +def thermodynamicJ1(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_1(m) + Eq. (11) in https://www.overleaf.com/project/674f23fdff4ede2888cab91e + """ + res = 0. + truncOrder = 5 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) + res += sign**(k-1) * m**2 * mpmath.besselk(2, arg) + return res + + +def thermodynamicJ2(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_2(m) + Eq. (12) in https://www.overleaf.com/project/674f23fdff4ede2888cab91e + """ + res = 0. + truncOrder = 5 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) + res += sign**(k-1) * m**2 / k * ( + arg * mpmath.besselk(1, arg) + 3. * mpmath.besselk(2, arg) + ) + return res + + +def thermodynamicJ3(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_3(m) + Eq. (13) in https://www.overleaf.com/project/674f23fdff4ede2888cab91e + """ + res = 0. + truncOrder = 5 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) + res += sign**(k-1) * m**3 / k * ( + arg * mpmath.besselk(2, arg) + 3. * mpmath.besselk(3, arg) + ) + return res + + +def thermodynamicJ0p5(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_0p5(m) + """ + res = 0. + if m < 5.: + truncOrder = 5 + else: + truncOrder = 2 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) / 2. + res += sign**(k-1) * k * np.sqrt(np.pi) / (16. * k**2.5) * ( + np.sqrt(2.) * k**2 * m**2 * np.pi * mpmath.besseli(-0.25, arg)**2 + + np.sqrt(2.) * np.pi * ( + (-3 + k**2 * m**2) * mpmath.besseli(0.25, arg)**2 + + 3 * k * m * mpmath.besseli(0.25, arg) * mpmath.besseli(0.75, arg) + + 3 * k * m * mpmath.besseli(-0.75, arg) * ( + mpmath.besseli(0.25, arg) + - 2 * k * m * mpmath.besseli(0.75, arg)) + + 3 * k**2 * m**2 * ( + mpmath.besseli(0.75, arg)**2 + + mpmath.besseli(1.25, arg)**2) + ) + - 2 * k * m * mpmath.besseli(-0.25, arg) * ( + np.sqrt(2) * k * m * np.pi * mpmath.besseli(0.25, arg) + - 3 * mpmath.besselk(-0.75, arg) + ) + ) + return res + + +def thermodynamicJ1p5(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_1p5(m) + """ + res = 0. + if m < 5.: + truncOrder = 5 + else: + truncOrder = 2 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) / 2. + res += sign**(k-1) * k * np.sqrt(np.pi) / (32. * k**3.5) * ( + 5 * np.sqrt(2.) * k**2 * m**2 * np.pi * mpmath.besseli(-0.25, arg)**2 + + np.sqrt(2.) * np.pi * ( + 5 * (-3 + k**2 * m**2) * mpmath.besseli(0.25, arg)**2 + + k * m * (15 + 8 * k**2 * m**2) * mpmath.besseli(0.25, arg) + * mpmath.besseli(0.75, arg) + + k * m * mpmath.besseli(-0.75, arg) * ( + (15 - 8 * k**2 * m**2) * mpmath.besseli(0.25, arg) + - 30 * k * m * mpmath.besseli(0.75, arg)) + + 15 * k**2 * m**2 * ( + mpmath.besseli(0.75, arg)**2 + + mpmath.besseli(1.25, arg)**2) + ) + + 2 * k * m * mpmath.besseli(-0.25, arg) * ( + - 5 * np.sqrt(2) * k * m * np.pi * mpmath.besseli(0.25, arg) + + (15 + 8 * k**2 * m**2) * mpmath.besselk(-0.75, arg) + ) + ) + return res + + +def thermodynamicJ2p5(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_2p5(m) + """ + res = 0. + if m < 5.: + truncOrder = 5 + else: + truncOrder = 2 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) / 2. + res += sign**(k-1) * k * np.sqrt(np.pi) / (64 * k**4.5) * ( + np.sqrt(2.) * k**2 * m**2 * np.pi * (35 + 8 * k**2 * m**2) + * mpmath.besseli(-0.25, arg)**2 + + np.sqrt(2.) * np.pi * ( + (-105 + 27 * k**2 * m**2 + 8 * k**4 * m**4) + * mpmath.besseli(0.25, arg)**2 + + k * m * (105 + 64 * k**2 * m**2) * mpmath.besseli(0.25, arg) + * mpmath.besseli(0.75, arg) + - k * m * mpmath.besseli(-0.75, arg) * ( + 3 * (-35 + 16 * k**2 * m**2) * mpmath.besseli(0.25, arg) + + 2 * k * m * (105 + 8 * k**2 * m**2) * mpmath.besseli(0.75, arg)) + + k**2 * m**2 * (105 + 8 * k**2 * m**2) * ( + mpmath.besseli(0.75, arg)**2 + + mpmath.besseli(1.25, arg)**2) + ) + - 2 * k * m * mpmath.besseli(-0.25, arg) * ( + np.sqrt(2) * k * m * (35 + 8 * k**2 * m**2) * np.pi * mpmath.besseli(0.25, arg) + - (105 + 64 * k**2 * m**2) * mpmath.besselk(-0.75, arg) + ) + ) + return res + + +def thermodynamicJ3p5(m: float, bosonFlag: bool) -> float: + """ + Compute thermodynamic integral J_3p5(m) + """ + res = 0. + if m < 5.: + truncOrder = 5 + else: + truncOrder = 2 + sign = 1. + if bosonFlag: + sign = -1. + for k in range(1, truncOrder): + arg = (k * m) / 2. + res += sign**(k-1) * k * np.sqrt(np.pi) / (128 * k**5.5) * ( + 5 * np.sqrt(2.) * k**2 * m**2 * np.pi * (63 + 16 * k**2 * m**2) + * mpmath.besseli(-0.25, arg)**2 + + np.sqrt(2.) * np.pi * ( + (-945 + 219 * k**2 * m**2 + 80 * k**4 * m**4) + * mpmath.besseli(0.25, arg)**2 + + k * m * (945 + 600 * k**2 * m**2 + 32 * k**4 * m**4) + * mpmath.besseli(0.25, arg) * mpmath.besseli(0.75, arg) + - k * m * mpmath.besseli(-0.75, arg) * ( + (-945 + 408 * k**2 * m**2 + 32 * k**4 * m**4) * mpmath.besseli(0.25, arg) + + 6 * k * m * (315 + 32 * k**2 * m**2) * mpmath.besseli(0.75, arg)) + + 3 * k**2 * m**2 * (315 + 32 * k**2 * m**2) * ( + mpmath.besseli(0.75, arg)**2 + + mpmath.besseli(1.25, arg)**2) + ) + - 2 * k * m * mpmath.besseli(-0.25, arg) * ( + 5 * np.sqrt(2) * k * m * (63 + 16 * k**2 * m**2) * np.pi * mpmath.besseli(0.25, arg) + - (945 + 600 * k**2 * m**2 + 32 * k**4 * m**4) * mpmath.besselk(-0.75, arg) + ) + ) + return res + +#print(computeThermodynamicJ(0.01, 0.5, False), thermodynamicJ0(0.5, False)) +#print(computeThermodynamicJ(0.01, 0.5, True), thermodynamicJ0(0.5, True)) +#print(computeThermodynamicJ(1.01, 0.5, False), thermodynamicJ1(0.5, False)) +#print(computeThermodynamicJ(1.01, 0.5, True), thermodynamicJ1(0.5, True)) +#print(computeThermodynamicJ(2.01, 0.5, False), thermodynamicJ2(0.5, False)) +#print(computeThermodynamicJ(2.01, 0.5, True), thermodynamicJ2(0.5, True)) +#print(computeThermodynamicJ(3.01, 0.5, False), thermodynamicJ3(0.5, False)) +#print(computeThermodynamicJ(3.01, 0.5, True), thermodynamicJ3(0.5, True)) + +m_min = 0.5 +m_max = 15. +dm = 0.05 +n_points = int((m_max - m_min) / dm) + 1 +mList = np.linspace(m_min, m_max, n_points) + +#for m_i in mList: +# print(f"{m_i:.6e}, {float(thermodynamicJ3p5(m_i, False)):.6e}, {float(computeThermodynamicJ(3.5, m_i, False)):.6e}") +#exit(0) + +gammaList = [0.0, 0.5, 1.0, 1.5] + +for bosonFlag in [True, False]: + if bosonFlag: + f1 = open("J_gamma_BE.dat", "w") + f2 = open("J_gammaplus2_BE.dat", "w") + f3 = open("J_2_BE.dat", "w") + else: + f1 = open("J_gamma_FD.dat", "w") + f2 = open("J_gammaplus2_FD.dat", "w") + f3 = open("J_2_FD.dat", "w") + f1.write(f"# m/T J_gamma (gamma = {gammaList})\n") + f2.write(f"# m/T J_gamma+2 (gamma = {gammaList})\n") + f3.write(f"# m/T J_2\n") + for m_i in mList: + J_2 = thermodynamicJ2(m_i, bosonFlag) + f3.write(f"{m_i:.8e} {float(J_2):.8e}\n") + f1.write(f"{m_i:.8e}") + f2.write(f"{m_i:.8e}") + for gamma_i in gammaList: + if abs(gamma_i - 0) < 1e-10: + J_gamma = thermodynamicJ0(m_i, bosonFlag) + J_gammaplus2 = thermodynamicJ2(m_i, bosonFlag) + elif abs(gamma_i - 0.5) < 1e-10: + J_gamma = thermodynamicJ0p5(m_i, bosonFlag) + J_gammaplus2 = thermodynamicJ2p5(m_i, bosonFlag) + elif abs(gamma_i - 1) < 1e-10: + J_gamma = thermodynamicJ1(m_i, bosonFlag) + J_gammaplus2 = thermodynamicJ3(m_i, bosonFlag) + elif abs(gamma_i - 1.5) < 1e-10: + J_gamma = thermodynamicJ1p5(m_i, bosonFlag) + J_gammaplus2 = thermodynamicJ3p5(m_i, bosonFlag) + else: + J_gamma = computeThermodynamicJ(gamma_i, m_i, bosonFlag) + J_gammaplus2 = computeThermodynamicJ(gamma_i + 2, m_i, bosonFlag) + f1.write(f" {float(J_gamma):.8e}") + f2.write(f" {float(J_gammaplus2):.8e}") + f1.write("\n") + f2.write("\n") + f1.close() + f2.close() + f3.close()