src/afterburner/SmashWrapper.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 // -----------------------------------------
0016 // This is a wrapper for SMASH hadronic afterburner with the JETSCAPE framework
0017 // -----------------------------------------
0018
0019 #include "SmashWrapper.h"
0020
0021 #include "smash/particles.h"
0022 #include "smash/library.h"
0023
0024 #include <string>
0025 #include <filesystem>
0026
0027 using namespace Jetscape;
0028
0029 // Register the module with the base class
0030 RegisterJetScapeModule<SmashWrapper> SmashWrapper::reg("SMASH");
0031
0032 SmashWrapper::SmashWrapper() { SetId("SMASH"); }
0033
0034 void SmashWrapper::InitTask() {
0035   JSINFO << "SMASH: picking SMASH-specific configuration from xml file";
0036   std::string smash_config_file =
0037       GetXMLElementText({"Afterburner", "SMASH", "SMASH_config_file"});
0038   std::string smash_hadron_list =
0039       GetXMLElementText({"Afterburner", "SMASH", "SMASH_particles_file"});
0040   std::string smash_decays_list =
0041       GetXMLElementText({"Afterburner", "SMASH", "SMASH_decaymodes_file"});
0042
0043   // do not store tabulation, which is achieved by an empty tabulations path
0044   std::string tabulations_path("");
0045
0046   auto config = smash::setup_config_and_logging(smash_config_file,
0047                                                 smash_hadron_list,
0048                                                 smash_decays_list);
0049
0050   // Take care of the random seed. This will make SMASH results reproducible.
0051   auto random_seed = (*GetMt19937Generator())();
0052   config.set_value({"General","Randomseed"}, random_seed);
0053
0054   // Read in the rest of configuration
0055   end_time_ = GetXMLElementDouble({"Afterburner", "SMASH", "end_time"});
0056   config.set_value({"General","End_Time"}, end_time_);
0057
0058   JSINFO << "End time until which SMASH propagates is " << end_time_ << " fm/c";
0059   only_final_decays_ =
0060       GetXMLElementInt({"Afterburner", "SMASH", "only_decays"});
0061   if (only_final_decays_) {
0062     JSINFO << "SMASH will only perform resonance decays, no propagation";
0063   }
0064
0065   const std::string smash_version(SMASH_VERSION);
0066   smash::initialize_particles_decays_and_tabulations(config, smash_version,
0067                                                      tabulations_path);
0068
0069   JSINFO << "Seting up SMASH Experiment object";
0070   // output path is just dummy here, because no output from SMASH is foreseen
0071   std::filesystem::path output_path("./smash_output");
0072   smash_experiment_ =
0073       make_shared<smash::Experiment<AfterburnerModus>>(config, output_path);
0074   JSINFO << "Finish initializing SMASH";
0075 }
0076
0077 void SmashWrapper::ExecuteTask() {
0078   AfterburnerModus *modus = smash_experiment_->modus();
0079   // This is necessary to correctly handle indices of particle sets from hydro.
0080   // Every hydro event creates a new structure like jetscape_hadrons_
0081   // with as many events in it as one has samples per hydro
0082   modus->reset_event_numbering();
0083   modus->jetscape_hadrons_ = GatherAfterburnerHadrons();
0084   const int n_events = modus->jetscape_hadrons_.size();
0085   JSINFO << "SMASH: obtained " << n_events << " events from particlization";
0086   // SMASH within JETSCAPE only works with one (the first) ensemble
0087   smash::Particles *smash_particles = smash_experiment_->first_ensemble();
0088   for (unsigned int i = 0; i < n_events; i++) {
0089     JSINFO << "Event " << i << " SMASH starts with "
0090            << modus->jetscape_hadrons_[i].size() << " particles.";
0091     smash_experiment_->initialize_new_event();
0092     if (!only_final_decays_) {
0093       smash_experiment_->run_time_evolution(end_time_);
0094     }
0095     smash_experiment_->do_final_decays();
0096     smash_experiment_->final_output();
0097     smash_particles_to_JS_hadrons(*smash_particles,
0098                                   modus->jetscape_hadrons_[i]);
0099     JSINFO << modus->jetscape_hadrons_[i].size() << " hadrons from SMASH.";
0100     smash_experiment_->increase_event_number();
0101   }
0102 }
0103
0104 void SmashWrapper::WriteTask(weak_ptr<JetScapeWriter> w) {
0105   JSINFO << "SMASH hadronic afterburner printout";
0106   auto f = w.lock();
0107   if (!f) {
0108     return;
0109   }
0110   AfterburnerModus *modus = smash_experiment_->modus();
0111   f->WriteComment("JetScape module: " + GetId());
0112   for (const auto &event : modus->jetscape_hadrons_) {
0113     int i = -1;
0114     for (const auto hadron : event) {
0115       f->WriteWhiteSpace("[" + to_string(++i) + "] H");
0116       f->Write(hadron);
0117     }
0118   }
0119 }
0120
0121 void AfterburnerModus::JS_hadrons_to_smash_particles(
0122     const std::vector<shared_ptr<Hadron>> &JS_hadrons,
0123     smash::Particles &smash_particles) {
0124   smash_particles.reset();
0125   for (const auto JS_hadron : JS_hadrons) {
0126     const FourVector p = JS_hadron->p_in();
0127     const FourVector r = JS_hadron->x_in();
0128     const double mass = JS_hadron->restmass();
0129     smash::PdgCode pdgcode = smash::PdgCode::from_decimal(JS_hadron->pid());
0130     this->try_create_particle(smash_particles, pdgcode, r.t(), r.x(), r.y(),
0131                               r.z(), mass, JS_hadron->e(), JS_hadron->px(),
0132                               JS_hadron->py(), JS_hadron->pz());
0133   }
0134 }
0135
0136 void SmashWrapper::smash_particles_to_JS_hadrons(
0137     const smash::Particles &smash_particles,
0138     std::vector<shared_ptr<Hadron>> &JS_hadrons) {
0139   JS_hadrons.clear();
0140   for (const auto &particle : smash_particles) {
0141     const int hadron_label = 0;
0142     const int hadron_status = 27;
0143     const int hadron_id = particle.pdgcode().get_decimal();
0144     smash::FourVector p = particle.momentum(), r = particle.position();
0145     const FourVector hadron_p(p.x1(), p.x2(), p.x3(), p.x0()),
0146         hadron_r(r.x1(), r.x2(), r.x3(), r.x0());
0147     const double hadron_mass = p.abs();
0148     JS_hadrons.push_back(make_shared<Hadron>(hadron_label, hadron_id,
0149                                              hadron_status, hadron_p, hadron_r,
0150                                              hadron_mass));
0151   }
0152 }