From e4512e2c4dfcc542c9d852dfdce7976bdac01791 Mon Sep 17 00:00:00 2001 From: Chun Shen Date: Tue, 16 Apr 2024 13:27:11 -0400 Subject: [PATCH 1/3] avoid nan or inf when E and p^z are too close for particle samples --- src/FSSW.cpp | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/src/FSSW.cpp b/src/FSSW.cpp index 34a7526..f1626d8 100644 --- a/src/FSSW.cpp +++ b/src/FSSW.cpp @@ -1954,12 +1954,30 @@ int FSSW::sample_momemtum_from_a_fluid_cell( // assigned to the return variables pT = sqrt(pLab[1]*pLab[1] + pLab[2]*pLab[2]); - if (pT < 5.) { - phi = atan2(pLab[2], pLab[1]); - double y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3])); - y_minus_eta_s = y - surf->eta; - return(1); + phi = atan2(pLab[2], pLab[1]); + double y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3])); + if (std::isnan(y) || std::isinf(y)) { + // pLab[0] and pLab[3] are too close + // recompute pLab[3] with E, pT, and mass + pLab[3] = sqrt(pLab[0]*pLab[0] - pT*pT - mass*mass); + y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3])); + if (std::isnan(y) || std::isinf(y)) { + std::cout << "[Error]: sampled y is " << y << "!" + << std::endl; + std::cout << "umu = " + << umu[0] << ", " << umu[1] << ", " << umu[2] + << ", " << umu[3] << std::endl; + std::cout << "pLRF = " + << pLRF[0] << ", " << pLRF[1] << ", " + << pLRF[2] << ", " << pLRF[3] << std::endl; + std::cout << "pLab = " + << pLab[0] << ", " << pLab[1] << ", " + << pLab[2] << ", " << pLab[3] << std::endl; + exit(1); + } } + y_minus_eta_s = y - surf->eta; + return(1); } tries++; } From 464627d531e03af5ead0002db730087404fbd732 Mon Sep 17 00:00:00 2001 From: Chun Shen Date: Tue, 16 Apr 2024 13:47:20 -0400 Subject: [PATCH 2/3] suppress the compiler warnings --- src/FSSW.cpp | 22 ++++----- src/FermionMomentumSampler.cpp | 2 - src/MomentumSamplerBase.h | 2 +- src/emissionfunction.cpp | 84 +++++++++++++++++----------------- src/iSS.h | 1 - src/readindata.cpp | 2 - 6 files changed, 54 insertions(+), 59 deletions(-) diff --git a/src/FSSW.cpp b/src/FSSW.cpp index f1626d8..a61cccf 100644 --- a/src/FSSW.cpp +++ b/src/FSSW.cpp @@ -469,17 +469,17 @@ void FSSW::combine_samples_to_OSCAR() { oscar << setw(10) << ipart + 1 << " " << setw(10) << (*(*Hadron_list)[iev])[ipart].pid << " "; - sprintf(line_buffer, - "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e", - (*(*Hadron_list)[iev])[ipart].px, - (*(*Hadron_list)[iev])[ipart].py, - (*(*Hadron_list)[iev])[ipart].pz, - (*(*Hadron_list)[iev])[ipart].E, - (*(*Hadron_list)[iev])[ipart].mass, - (*(*Hadron_list)[iev])[ipart].x, - (*(*Hadron_list)[iev])[ipart].y, - (*(*Hadron_list)[iev])[ipart].z, - (*(*Hadron_list)[iev])[ipart].t); + snprintf(line_buffer, 500, + "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e", + (*(*Hadron_list)[iev])[ipart].px, + (*(*Hadron_list)[iev])[ipart].py, + (*(*Hadron_list)[iev])[ipart].pz, + (*(*Hadron_list)[iev])[ipart].E, + (*(*Hadron_list)[iev])[ipart].mass, + (*(*Hadron_list)[iev])[ipart].x, + (*(*Hadron_list)[iev])[ipart].y, + (*(*Hadron_list)[iev])[ipart].z, + (*(*Hadron_list)[iev])[ipart].t); oscar << line_buffer << endl; } } diff --git a/src/FermionMomentumSampler.cpp b/src/FermionMomentumSampler.cpp index 4ee3c7c..ac54f5c 100644 --- a/src/FermionMomentumSampler.cpp +++ b/src/FermionMomentumSampler.cpp @@ -80,14 +80,12 @@ double FermionMomentumSampler::CDF_2(const double Etilde) const { double FermionMomentumSampler::CDF_2(const double Etilde, const double mtilde) const { double res = 0.; - double sign = 1.; for (int n = 0; n < trunc_order_; n++) { int n1 = n + 1; res += (1./(n1*n1*n1)*exp(-mtilde*n)*( (mtilde*n1*(mtilde*n1 + 2) + 2) - exp((mtilde - Etilde)*n1)*(Etilde*n1*(Etilde*n1 + 2) + 2)) ); - sign *= -1.; } return(res); } diff --git a/src/MomentumSamplerBase.h b/src/MomentumSamplerBase.h index ade9423..b9d2cf8 100644 --- a/src/MomentumSamplerBase.h +++ b/src/MomentumSamplerBase.h @@ -29,7 +29,7 @@ class MomentumSamplerBase { std::vector CDF_2_; MomentumSamplerBase(); - ~MomentumSamplerBase() {}; + virtual ~MomentumSamplerBase() {}; virtual double CDF_0(const double Etilde) const {return(0.0*Etilde);} virtual double CDF_1(const double Etilde) const {return(0.0*Etilde);} diff --git a/src/emissionfunction.cpp b/src/emissionfunction.cpp index 3ccb70f..f14e65b 100644 --- a/src/emissionfunction.cpp +++ b/src/emissionfunction.cpp @@ -1344,7 +1344,7 @@ void EmissionFunctionArray::sample_using_dN_dxtdetady_smooth_pT_phi() { long number_to_sample = determine_number_to_sample( dN, sampling_model, sampling_para1); // write to control file - sprintf(line_buffer, "%lu\n", number_to_sample); + snprintf(line_buffer, 500, "%lu\n", number_to_sample); control_str_buffer << line_buffer; control_writing_signal++; if (control_writing_signal == iSS_data::NUMBER_OF_LINES_TO_WRITE) { @@ -1539,18 +1539,18 @@ void EmissionFunctionArray::sample_using_dN_dxtdetady_smooth_pT_phi() { // write to sample file if (!USE_OSCAR_FORMAT) { - sprintf(line_buffer, - "%lu %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", - FO_idx, surf->tau, surf->xpt, surf->ypt, y_minus_eta_s, - pT, phi, surf->da0, surf->da1, surf->da2, vx, vy, y, - eta_s, E, p_z, t, z); + snprintf(line_buffer, 500, + "%lu %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", + FO_idx, surf->tau, surf->xpt, surf->ypt, y_minus_eta_s, + pT, phi, surf->da0, surf->da1, surf->da2, vx, vy, y, + eta_s, E, p_z, t, z); } // To be combined to OSCAR if (USE_OSCAR_FORMAT) { - sprintf(line_buffer, - "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e\n", - px, py, p_z, E, mass, surf->xpt, surf->ypt, z, t); + snprintf(line_buffer, 500, + "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e\n", + px, py, p_z, E, mass, surf->xpt, surf->ypt, z, t); } sample_str_buffer << line_buffer; sample_writing_signal++; @@ -1787,7 +1787,7 @@ void EmissionFunctionArray::sample_using_dN_pTdpTdphidy() { long number_to_sample = determine_number_to_sample( dN, sampling_model, sampling_para1); // write to control file - sprintf(line_buffer, "%lu\n", number_to_sample); + snprintf(line_buffer, 500, "%lu\n", number_to_sample); control_str_buffer << line_buffer; control_writing_signal++; if (control_writing_signal == iSS_data::NUMBER_OF_LINES_TO_WRITE) { @@ -2050,17 +2050,17 @@ void EmissionFunctionArray::sample_using_dN_pTdpTdphidy() { // write to sample file if (!USE_OSCAR_FORMAT) { - sprintf(line_buffer, - "%lu %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", - FO_idx, surf->tau, surf->xpt, surf->ypt, y_minus_eta_s, - pT, phi, surf->da0, surf->da1, surf->da2, - surf->u1/surf->u0, surf->u2/surf->u0, y, eta_s, E, - p_z, t, z); + snprintf(line_buffer, 500, + "%lu %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", + FO_idx, surf->tau, surf->xpt, surf->ypt, y_minus_eta_s, + pT, phi, surf->da0, surf->da1, surf->da2, + surf->u1/surf->u0, surf->u2/surf->u0, y, eta_s, E, + p_z, t, z); } else { // To be combined to OSCAR - sprintf(line_buffer, - "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e\n", - px, py, p_z, E, mass, surf->xpt, surf->ypt, z, t); + snprintf(line_buffer, 500, + "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e\n", + px, py, p_z, E, mass, surf->xpt, surf->ypt, z, t); } sample_str_buffer << line_buffer; sample_writing_signal++; @@ -2699,17 +2699,17 @@ void EmissionFunctionArray::combine_samples_to_OSCAR() { oscar << setw(10) << ipart + 1 << " " << setw(10) << (*(*Hadron_list)[iev])[ipart].pid << " "; - sprintf(line_buffer, - "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e", - (*(*Hadron_list)[iev])[ipart].px, - (*(*Hadron_list)[iev])[ipart].py, - (*(*Hadron_list)[iev])[ipart].pz, - (*(*Hadron_list)[iev])[ipart].E, - (*(*Hadron_list)[iev])[ipart].mass, - (*(*Hadron_list)[iev])[ipart].x, - (*(*Hadron_list)[iev])[ipart].y, - (*(*Hadron_list)[iev])[ipart].z, - (*(*Hadron_list)[iev])[ipart].t); + snprintf(line_buffer, 500, + "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e", + (*(*Hadron_list)[iev])[ipart].px, + (*(*Hadron_list)[iev])[ipart].py, + (*(*Hadron_list)[iev])[ipart].pz, + (*(*Hadron_list)[iev])[ipart].E, + (*(*Hadron_list)[iev])[ipart].mass, + (*(*Hadron_list)[iev])[ipart].x, + (*(*Hadron_list)[iev])[ipart].y, + (*(*Hadron_list)[iev])[ipart].z, + (*(*Hadron_list)[iev])[ipart].t); oscar << line_buffer << endl; } } @@ -3442,7 +3442,7 @@ void EmissionFunctionArray::sample_using_dN_dxtdy_4all_particles_conventional() control_writing_signal++; if (flag_output_samples_into_files == 1) { // write to control file - sprintf(line_buffer, "%lu\n", number_to_sample); + snprintf(line_buffer, 500, "%lu\n", number_to_sample); control_str_buffer << line_buffer; if (control_writing_signal == iSS_data::NUMBER_OF_LINES_TO_WRITE) { @@ -4440,19 +4440,19 @@ std::string EmissionFunctionArray::add_one_sampled_particle( // write to sample file if (!USE_OSCAR_FORMAT) { - sprintf(line_buffer, - "%lu %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", - FO_idx, surf->tau, surf->xpt, surf->ypt, - y_minus_eta_s, pT, phi, surf->da0, surf->da1, - surf->da2, surf->u1/surf->u0, - surf->u2/surf->u0, rapidity_y, eta_s, E, p_z, - t, z); + snprintf(line_buffer, 500, + "%lu %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", + FO_idx, surf->tau, surf->xpt, surf->ypt, + y_minus_eta_s, pT, phi, surf->da0, surf->da1, + surf->da2, surf->u1/surf->u0, + surf->u2/surf->u0, rapidity_y, eta_s, E, p_z, + t, z); } else { // To be combined to OSCAR - sprintf(line_buffer, - "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e\n", - px, py, p_z, E, mass, surf->xpt, - surf->ypt, z, t); + snprintf(line_buffer, 500, + "%24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e\n", + px, py, p_z, E, mass, surf->xpt, + surf->ypt, z, t); } text_string = line_buffer; diff --git a/src/iSS.h b/src/iSS.h index 8701207..dc6d78c 100644 --- a/src/iSS.h +++ b/src/iSS.h @@ -27,7 +27,6 @@ class iSS { std::vector FOsurf_Q_; std::shared_ptr ran_gen_ptr_; - int Nparticle; int flag_PCE_; AfterburnerType afterburner_type_; diff --git a/src/readindata.cpp b/src/readindata.cpp index 406b3d7..01098ca 100644 --- a/src/readindata.cpp +++ b/src/readindata.cpp @@ -608,7 +608,6 @@ void read_FOdata::read_FOsurfdat_hydro_analysis_boost_invariant( string input; double temp_tau, temp_xpt, temp_ypt; double temp_vx, temp_vy; - int idx = 0; surfdat_stream << path_ << "/hyper_surface_2+1d.dat"; std::ifstream surfdat(surfdat_stream.str().c_str()); getline(surfdat, input, '\n' ); @@ -670,7 +669,6 @@ void read_FOdata::read_FOsurfdat_hydro_analysis_boost_invariant( surf_ptr.push_back(surf_elem); - idx++; getline(surfdat, input, '\n' ); } surfdat.close(); From 294a234029ada9ba3954bef73592037670435700 Mon Sep 17 00:00:00 2001 From: Chun Shen Date: Tue, 16 Apr 2024 14:05:51 -0400 Subject: [PATCH 3/3] incrase precision for particle momentum vector --- src/FSSW.cpp | 42 ++++++++++++------------------------------ src/FSSW.h | 6 +++--- src/data_struct.h | 1 + 3 files changed, 16 insertions(+), 33 deletions(-) diff --git a/src/FSSW.cpp b/src/FSSW.cpp index a61cccf..2bd7843 100644 --- a/src/FSSW.cpp +++ b/src/FSSW.cpp @@ -27,7 +27,7 @@ #include "Stopwatch.h" using iSS_data::AMOUNT_OF_OUTPUT; -using iSS_data::Vec4; +using iSS_data::VecD4; using std::cout; using std::endl; using std::string; @@ -1945,39 +1945,21 @@ int FSSW::sample_momemtum_from_a_fluid_cell( if (ran_gen_ptr->rand_uniform() < accept_prob) { // accept the sample // now we need to boost the momentum to the lab frame - Vec4 pLRF = {static_cast(p0), static_cast(px), - static_cast(py), static_cast(pz)}; - Vec4 umu = {surf->u_tz[0], surf->u_tz[1], - surf->u_tz[2], surf->u_tz[3]}; - Vec4 pLab = {0., 0., 0., 0.}; + VecD4 pLRF = {p0, px, py, pz}; + VecD4 umu = {surf->u_tz[0], surf->u_tz[1], + surf->u_tz[2], surf->u_tz[3]}; + VecD4 pLab = {0., 0., 0., 0.}; boost_vector_back_to_lab_frame(pLRF, pLab, umu); // assigned to the return variables - pT = sqrt(pLab[1]*pLab[1] + pLab[2]*pLab[2]); - phi = atan2(pLab[2], pLab[1]); + // if E and p^z are too close, resample one double y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3])); - if (std::isnan(y) || std::isinf(y)) { - // pLab[0] and pLab[3] are too close - // recompute pLab[3] with E, pT, and mass - pLab[3] = sqrt(pLab[0]*pLab[0] - pT*pT - mass*mass); - y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3])); - if (std::isnan(y) || std::isinf(y)) { - std::cout << "[Error]: sampled y is " << y << "!" - << std::endl; - std::cout << "umu = " - << umu[0] << ", " << umu[1] << ", " << umu[2] - << ", " << umu[3] << std::endl; - std::cout << "pLRF = " - << pLRF[0] << ", " << pLRF[1] << ", " - << pLRF[2] << ", " << pLRF[3] << std::endl; - std::cout << "pLab = " - << pLab[0] << ", " << pLab[1] << ", " - << pLab[2] << ", " << pLab[3] << std::endl; - exit(1); - } + if (!std::isnan(y) && !std::isinf(y)) { + pT = sqrt(pLab[1]*pLab[1] + pLab[2]*pLab[2]); + phi = atan2(pLab[2], pLab[1]); + y_minus_eta_s = y - surf->eta; + return(1); } - y_minus_eta_s = y - surf->eta; - return(1); } tries++; } @@ -2018,7 +2000,7 @@ void FSSW::add_one_sampled_particle( // this function boost a 4-vector from the fluid local rest frame to // the lab frame. The fluid velocity is u^\mu. void FSSW::boost_vector_back_to_lab_frame( - Vec4 &p_LRF, Vec4 &p_lab, Vec4 &umu) const { + VecD4 &p_LRF, VecD4 &p_lab, VecD4 &umu) const { double p_dot_u = 0.; for (int i = 1; i < 4; i++) p_dot_u += p_LRF[i]*umu[i]; diff --git a/src/FSSW.h b/src/FSSW.h index 6676060..c26c236 100644 --- a/src/FSSW.h +++ b/src/FSSW.h @@ -207,9 +207,9 @@ class FSSW { const int particle_monval, const double mass, const double pT, const double phi, const double y_minus_eta_s, const double eta_s); - void boost_vector_back_to_lab_frame(iSS_data::Vec4 &p_LRF, - iSS_data::Vec4 &p_lab, - iSS_data::Vec4 &umu) const; + void boost_vector_back_to_lab_frame(iSS_data::VecD4 &p_LRF, + iSS_data::VecD4 &p_lab, + iSS_data::VecD4 &umu) const; double bilinearInterp(std::vector>&mat, int idx_e1, int idx_e2, int idx_nB1, int idx_nB2, double x_fraction, double y_fraction); diff --git a/src/data_struct.h b/src/data_struct.h index 529da97..36c0583 100644 --- a/src/data_struct.h +++ b/src/data_struct.h @@ -9,6 +9,7 @@ namespace iSS_data { const double hbarC=0.197327053; //GeV*fm typedef std::array Vec4; + typedef std::array VecD4; typedef std::array ViscousVec; const int NUMBER_OF_LINES_TO_WRITE = 100000;