src/framework/LiquefierBase.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015 #include "LiquefierBase.h"
0016 #include <math.h>
0017
0018 namespace Jetscape {
0019
0020 LiquefierBase::LiquefierBase()
0021     : hydro_source_abs_err(1e-10), drop_stat(-11), miss_stat(-13),
0022       neg_stat(-17) {
0023   GetHydroCellSignalConnected = false;
0024 }
0025
0026 void LiquefierBase::get_source(Jetscape::real tau, Jetscape::real x,
0027                                Jetscape::real y, Jetscape::real eta,
0028                                std::array<Jetscape::real, 4> &jmu) const {
0029   jmu = {0.0, 0.0, 0.0, 0.0};
0030   for (const auto &drop_i : dropletlist) {
0031
0032     const auto x_drop = drop_i.get_xmu();
0033     double ds2 = tau * tau + x_drop[0] * x_drop[0] -
0034                  2.0 * tau * x_drop[0] * cosh(eta - x_drop[3]) -
0035                  (x - x_drop[1]) * (x - x_drop[1]) -
0036                  (y - x_drop[2]) * (y - x_drop[2]);
0037
0038     if (tau >= x_drop[0] && ds2 >= 0.0) {
0039       std::array<Jetscape::real, 4> jmu_i = {0.0, 0.0, 0.0, 0.0};
0040       smearing_kernel(tau, x, y, eta, drop_i, jmu_i);
0041       for (int i = 0; i < 4; i++)
0042         jmu[i] += jmu_i[i];
0043     }
0044   }
0045 }
0046
0047 //! This function check the energy momentum conservation at the vertex
0048 //! If vertex does not conserve energy and momentum,
0049 //! a p_missing will be added to the pOut list
0050 void LiquefierBase::check_energy_momentum_conservation(
0051     const std::vector<Parton> &pIn, std::vector<Parton> &pOut) {
0052   FourVector p_init(0., 0., 0., 0.);
0053   for (const auto &iparton : pIn) {
0054     auto temp = iparton.p_in();
0055     if (iparton.pstat() == -1) {
0056       p_init -= temp;
0057     } else {
0058       p_init += temp;
0059     }
0060   }
0061
0062   FourVector p_final(0., 0., 0., 0.);
0063   FourVector x_final(0., 0., 0., 0.);
0064   for (const auto &iparton : pOut) {
0065     auto temp = iparton.p_in();
0066     if (iparton.pstat() == -1) {
0067       p_final -= temp;
0068     } else {
0069       p_final += temp;
0070     }
0071     x_final = iparton.x_in();
0072   }
0073
0074   FourVector p_missing = p_init;
0075   p_missing -= p_final;
0076   if (std::abs(p_missing.t()) > hydro_source_abs_err ||
0077       std::abs(p_missing.x()) > hydro_source_abs_err ||
0078       std::abs(p_missing.y()) > hydro_source_abs_err ||
0079       std::abs(p_missing.z()) > hydro_source_abs_err) {
0080     JSWARN << "A vertex does not conserve energy momentum!";
0081     JSWARN << "E = " << p_missing.t() << " GeV, px = " << p_missing.x()
0082            << " GeV, py = " << p_missing.y() << " GeV, pz = " << p_missing.z()
0083            << " GeV.";
0084     Parton parton_miss(0, 21, miss_stat, p_missing, x_final);
0085     pOut.push_back(parton_miss);
0086   }
0087 }
0088
0089 void LiquefierBase::filter_partons(std::vector<Parton> &pOut) {
0090   // if e_threshold > 0, use e_threshold, else, use |e_threshold|*T
0091   // this should be put into xml later.
0092   auto e_threshold = 2.0; // GeV
0093
0094   for (auto &iparton : pOut) {
0095     if (iparton.pstat() == miss_stat)
0096       continue;
0097
0098     // ignore photons
0099     if (iparton.isPhoton(iparton.pid()))
0100       continue;
0101
0102     // ignore heavy quarks
0103     if (std::abs(iparton.pid()) == 4 || std::abs(iparton.pid()) == 5)
0104       continue;
0105
0106     if (iparton.pstat() == -1) {
0107       // remove negative particles from parton list
0108       //iparton.set_stat(drop_stat);
0109       iparton.set_stat(neg_stat);
0110       continue;
0111     }
0112
0113     // for positive particles, including jet partons and recoil partons
0114     auto tLoc = iparton.x_in().t();
0115     auto xLoc = iparton.x_in().x();
0116     auto yLoc = iparton.x_in().y();
0117     auto zLoc = iparton.x_in().z();
0118     std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0119     GetHydroCellSignal(tLoc, xLoc, yLoc, zLoc, check_fluid_info_ptr);
0120     auto vxLoc = check_fluid_info_ptr->vx;
0121     auto vyLoc = check_fluid_info_ptr->vy;
0122     auto vzLoc = check_fluid_info_ptr->vz;
0123     auto beta2 = vxLoc * vxLoc + vyLoc * vyLoc + vzLoc * vzLoc;
0124     auto gamma = 1.0 / sqrt(1.0 - beta2);
0125     auto E_boosted = gamma * (iparton.e() - iparton.p(1) * vxLoc -
0126                               iparton.p(2) * vyLoc - iparton.p(3) * vzLoc);
0127     if (e_threshold > 0.) {
0128       // drop partons with energy smaller than e_threshold
0129       // (in the local rest frame) from parton list
0130       if (E_boosted < e_threshold) {
0131         iparton.set_stat(drop_stat);
0132         continue;
0133       }
0134     } else {
0135       // drop partons with energy smaller than |e_threshold|*T
0136       // (in the local rest frame) from parton list
0137       auto tempLoc = check_fluid_info_ptr->temperature;
0138       if (E_boosted < std::abs(e_threshold) * tempLoc) {
0139         iparton.set_stat(drop_stat);
0140         continue;
0141       }
0142     }
0143   }
0144 }
0145
0146 void LiquefierBase::add_hydro_sources(std::vector<Parton> &pIn,
0147                                       std::vector<Parton> &pOut) {
0148   if (pOut.size() == 0) {
0149     // the process is freestreaming, ignore
0150     filter_partons(pIn);
0151     return;
0152   }
0153   check_energy_momentum_conservation(pIn, pOut);
0154   filter_partons(pOut);
0155
0156   FourVector p_final;
0157   FourVector p_init;
0158   FourVector x_final;
0159   FourVector x_init;
0160   //cout << "debug, mid ......." << pOut.size() << endl;
0161   // use energy conservation to deterime the source term
0162   const auto weight_init = 1.0;
0163   for (const auto &iparton : pIn) {
0164     auto temp = iparton.p_in();
0165     p_init += temp;
0166     x_init = iparton.x_in();
0167   }
0168
0169   auto weight_final = 0.0;
0170   for (const auto &iparton : pOut) {
0171     if (iparton.pstat() == drop_stat)
0172       continue;
0173     if (iparton.pstat() == miss_stat)
0174       continue;
0175     if (iparton.pstat() == neg_stat)
0176       continue;
0177     auto temp = iparton.p_in();
0178     p_final += temp;
0179     x_final = iparton.x_in();
0180     weight_final = 1.0;
0181   }
0182
0183   if (std::abs(p_init.t() - p_final.t()) / p_init.t() > hydro_source_abs_err ||
0184       std::abs(p_init.x() - p_final.x()) / p_init.t() > hydro_source_abs_err ||
0185       std::abs(p_init.y() - p_final.y()) / p_init.t() > hydro_source_abs_err ||
0186       std::abs(p_init.z() - p_final.z()) / p_init.t() > hydro_source_abs_err) {
0187     auto droplet_t = ((x_final.t() * weight_final + x_init.t() * weight_init) /
0188                       (weight_final + weight_init));
0189     auto droplet_x = ((x_final.x() * weight_final + x_init.x() * weight_init) /
0190                       (weight_final + weight_init));
0191     auto droplet_y = ((x_final.y() * weight_final + x_init.y() * weight_init) /
0192                       (weight_final + weight_init));
0193     auto droplet_z = ((x_final.z() * weight_final + x_init.z() * weight_init) /
0194                       (weight_final + weight_init));
0195
0196     auto droplet_tau = sqrt(droplet_t * droplet_t - droplet_z * droplet_z);
0197     auto droplet_eta =
0198         (0.5 * log((droplet_t + droplet_z) / (droplet_t - droplet_z)));
0199     auto droplet_E = p_init.t() - p_final.t();
0200     auto droplet_px = p_init.x() - p_final.x();
0201     auto droplet_py = p_init.y() - p_final.y();
0202     auto droplet_pz = p_init.z() - p_final.z();
0203
0204     std::array<Jetscape::real, 4> droplet_xmu = {
0205         static_cast<Jetscape::real>(droplet_tau),
0206         static_cast<Jetscape::real>(droplet_x),
0207         static_cast<Jetscape::real>(droplet_y),
0208         static_cast<Jetscape::real>(droplet_eta)};
0209     std::array<Jetscape::real, 4> droplet_pmu = {
0210         static_cast<Jetscape::real>(droplet_E),
0211         static_cast<Jetscape::real>(droplet_px),
0212         static_cast<Jetscape::real>(droplet_py),
0213         static_cast<Jetscape::real>(droplet_pz)};
0214     Droplet drop_i(droplet_xmu, droplet_pmu);
0215     add_a_droplet(drop_i);
0216   }
0217 }
0218
0219 void LiquefierBase::Clear() { dropletlist.clear(); }
0220
0221 Jetscape::real LiquefierBase::get_dropletlist_total_energy() const {
0222   Jetscape::real total_E = 0.0;
0223   for (const auto &drop_i : dropletlist) {
0224     total_E += drop_i.get_pmu()[0];
0225   }
0226   return (total_E);
0227 }
0228
0229 }; // namespace Jetscape