src/hadronization/iSpectraSamplerWrapper.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 iSpectraSampler (iSS) with the JETSCAPE framework
0017 // -----------------------------------------
0018
0019 #include "JetScapeLogger.h"
0020 #include "iSpectraSamplerWrapper.h"
0021
0022 #include <memory>
0023 #include <string>
0024 #include <fstream>
0025
0026 using namespace Jetscape;
0027
0028 // Register the module with the base class
0029 RegisterJetScapeModule<iSpectraSamplerWrapper>
0030     iSpectraSamplerWrapper::reg("iSS");
0031
0032 iSpectraSamplerWrapper::iSpectraSamplerWrapper() { SetId("iSS"); }
0033
0034 iSpectraSamplerWrapper::~iSpectraSamplerWrapper() {}
0035
0036 void iSpectraSamplerWrapper::InitTask() {
0037
0038   JSINFO << "Initialize a particle sampler (iSS)";
0039
0040   std::string input_file =
0041       GetXMLElementText({"SoftParticlization", "iSS", "iSS_input_file"});
0042   std::string table_path =
0043       GetXMLElementText({"SoftParticlization", "iSS", "iSS_table_path"});
0044   std::string particle_table_path =
0045       GetXMLElementText({"SoftParticlization", "iSS",
0046                          "iSS_particle_table_path"});
0047   std::string working_path =
0048       GetXMLElementText({"SoftParticlization", "iSS", "iSS_working_path"});
0049   int hydro_mode =
0050       GetXMLElementInt({"SoftParticlization", "iSS", "hydro_mode"});
0051   int number_of_repeated_sampling = GetXMLElementInt(
0052       {"SoftParticlization", "iSS", "number_of_repeated_sampling"});
0053   int flag_perform_decays = GetXMLElementInt(
0054       {"SoftParticlization", "iSS", "Perform_resonance_decays"});
0055   int afterburner_type = (
0056       GetXMLElementInt({"SoftParticlization", "iSS", "afterburner_type"}));
0057
0058   if (!boost_invariance) {
0059     hydro_mode = 2;
0060   }
0061
0062   iSpectraSampler_ptr_ = std::unique_ptr<iSS>(
0063           new iSS(working_path, table_path, particle_table_path, input_file));
0064   iSpectraSampler_ptr_->paraRdr_ptr->readFromFile(input_file);
0065
0066   // overwrite some parameters
0067   int echoLevel = GetXMLElementInt({"vlevel"});
0068   iSpectraSampler_ptr_->paraRdr_ptr->setVal("JSechoLevel", echoLevel);
0069
0070   iSpectraSampler_ptr_->paraRdr_ptr->setVal("hydro_mode", hydro_mode);
0071   iSpectraSampler_ptr_->paraRdr_ptr->setVal("afterburner_type",
0072                                             afterburner_type);
0073   iSpectraSampler_ptr_->paraRdr_ptr->setVal("output_samples_into_files", 0);
0074   iSpectraSampler_ptr_->paraRdr_ptr->setVal("use_OSCAR_format", 0);
0075   iSpectraSampler_ptr_->paraRdr_ptr->setVal("use_gzip_format", 0);
0076   iSpectraSampler_ptr_->paraRdr_ptr->setVal("use_binary_format", 0);
0077   iSpectraSampler_ptr_->paraRdr_ptr->setVal("store_samples_in_memory", 1);
0078   iSpectraSampler_ptr_->paraRdr_ptr->setVal("number_of_repeated_sampling",
0079                                             number_of_repeated_sampling);
0080   iSpectraSampler_ptr_->paraRdr_ptr->setVal("perform_decays",
0081                                             flag_perform_decays);
0082
0083   // set default parameters
0084   iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_shear", 1);
0085   iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_bulk", 0);
0086   iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_rhob", 0);
0087   iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_diff", 0);
0088
0089   iSpectraSampler_ptr_->paraRdr_ptr->setVal("include_deltaf_shear", 1);
0090   iSpectraSampler_ptr_->paraRdr_ptr->setVal("include_deltaf_bulk", 0);
0091   iSpectraSampler_ptr_->paraRdr_ptr->setVal("bulk_deltaf_kind", 1);
0092   iSpectraSampler_ptr_->paraRdr_ptr->setVal("include_deltaf_diffusion", 0);
0093
0094   iSpectraSampler_ptr_->paraRdr_ptr->setVal("restrict_deltaf", 0);
0095   iSpectraSampler_ptr_->paraRdr_ptr->setVal("deltaf_max_ratio", 1.0);
0096   iSpectraSampler_ptr_->paraRdr_ptr->setVal("f0_is_not_small", 1);
0097
0098   iSpectraSampler_ptr_->paraRdr_ptr->setVal("calculate_vn", 0);
0099   iSpectraSampler_ptr_->paraRdr_ptr->setVal("MC_sampling", 4);
0100
0101   iSpectraSampler_ptr_->paraRdr_ptr->setVal(
0102       "sample_upto_desired_particle_number", 0);
0103   iSpectraSampler_ptr_->paraRdr_ptr->echo();
0104 }
0105
0106 void iSpectraSamplerWrapper::Exec() {
0107   JSINFO << "running iSS ...";
0108
0109   // generate symbolic links with music_input_file
0110   std::string music_input_file_path = GetXMLElementText(
0111           {"Hydro", "MUSIC", "MUSIC_input_file"});
0112   std::string working_path =
0113       GetXMLElementText({"SoftParticlization", "iSS", "iSS_working_path"});
0114   std::string music_input = working_path + "/music_input";
0115   std::ifstream inputfile(music_input.c_str());
0116   if (!inputfile.good()) {
0117     std::ostringstream system_command;
0118     system_command << "ln -s " << music_input_file_path << " "
0119                    << music_input;
0120     system(system_command.str().c_str());
0121   }
0122   inputfile.close();
0123
0124   int status = iSpectraSampler_ptr_->read_in_FO_surface();
0125   if (status != 0) {
0126     JSWARN << "Some errors happened in reading in the hyper-surface";
0127     exit(-1);
0128   }
0129
0130   auto random_seed = (*GetMt19937Generator())(); // get random seed
0131   iSpectraSampler_ptr_->set_random_seed(random_seed);
0132   VERBOSE(2) << "Random seed used for the iSS module" << random_seed;
0133
0134   status = iSpectraSampler_ptr_->generate_samples();
0135   if (status != 0) {
0136     JSWARN << "Some errors happened in generating particle samples";
0137     exit(-1);
0138   }
0139   PassHadronListToJetscape();
0140   JSINFO << "iSS finished.";
0141 }
0142
0143 void iSpectraSamplerWrapper::Clear() {
0144   VERBOSE(2) << "Finish the particle sampling";
0145   iSpectraSampler_ptr_->clear();
0146   for (unsigned i = 0; i < Hadron_list_.size(); i++) {
0147     Hadron_list_.at(i).clear();
0148   }
0149   Hadron_list_.clear();
0150 }
0151
0152 void iSpectraSamplerWrapper::PassHadronListToJetscape() {
0153   unsigned int nev = iSpectraSampler_ptr_->get_number_of_sampled_events();
0154   VERBOSE(2) << "Passing all sampled hadrons to the JETSCAPE framework";
0155   VERBOSE(4) << "number of events to pass : " << nev;
0156   for (unsigned int iev = 0; iev < nev; iev++) {
0157     std::vector<shared_ptr<Hadron>> hadrons;
0158     unsigned int nparticles =
0159         (iSpectraSampler_ptr_->get_number_of_particles(iev));
0160     VERBOSE(4) << "event " << iev << ": number of particles = " << nparticles;
0161     for (unsigned int ipart = 0; ipart < nparticles; ipart++) {
0162       iSS_Hadron current_hadron =
0163           (iSpectraSampler_ptr_->get_hadron(iev, ipart));
0164       int hadron_label = 0;
0165       int hadron_status = 11;
0166       int hadron_id = current_hadron.pid;
0167       //int hadron_id = 1;   // just for testing need to be changed to the line above
0168       double hadron_mass = current_hadron.mass;
0169       FourVector hadron_p(current_hadron.px, current_hadron.py,
0170                           current_hadron.pz, current_hadron.E);
0171       FourVector hadron_x(current_hadron.x, current_hadron.y, current_hadron.z,
0172                           current_hadron.t);
0173
0174       // create a JETSCAPE Hadron
0175       hadrons.push_back(make_shared<Hadron>(hadron_label, hadron_id,
0176                                             hadron_status, hadron_p, hadron_x,
0177                                             hadron_mass));
0178       //Hadron* jetscape_hadron = new Hadron(hadron_label, hadron_id, hadron_status, hadron_p, hadron_x, hadron_mass);
0179       //(*Hadron_list_)[iev]->push_back(*jetscape_hadron);
0180     }
0181     Hadron_list_.push_back(hadrons);
0182   }
0183   VERBOSE(4) << "JETSCAPE received " << Hadron_list_.size() << " events.";
0184   for (unsigned int iev = 0; iev < Hadron_list_.size(); iev++) {
0185     VERBOSE(4) << "In event " << iev << " JETSCAPE received "
0186                << Hadron_list_.at(iev).size() << " particles.";
0187   }
0188 }
0189
0190 void iSpectraSamplerWrapper::WriteTask(weak_ptr<JetScapeWriter> w) {
0191   VERBOSE(4) << "In iSpectraSamplerWrapper::WriteTask";
0192   auto f = w.lock();
0193   if (!f)
0194     return;
0195
0196   f->WriteComment("JetScape module: " + GetId());
0197   if (Hadron_list_.size() > 0) {
0198     f->WriteComment("Final State Bulk Hadrons");
0199     for (unsigned int j = 0; j < Hadron_list_.size(); j++) {
0200       vector<shared_ptr<Hadron>> hadVec = Hadron_list_.at(j);
0201       for (unsigned int i = 0; i < hadVec.size(); i++) {
0202         f->WriteWhiteSpace("[" + to_string(i) + "] H");
0203         f->Write(hadVec.at(i));
0204       }
0205     }
0206   } else {
0207     f->WriteComment("There are no bulk Hadrons");
0208   }
0209 }