Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
iSpectraSamplerWrapper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file iSpectraSamplerWrapper.cc
1 /*******************************************************************************
2  * Copyright (c) The JETSCAPE Collaboration, 2018
3  *
4  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
5  *
6  * For the list of contributors see AUTHORS.
7  *
8  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
9  *
10  * or via email to bugs.jetscape@gmail.com
11  *
12  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
13  * See COPYING for details.
14  ******************************************************************************/
15 // -----------------------------------------
16 // This is a wrapper for iSpectraSampler (iSS) with the JETSCAPE framework
17 // -----------------------------------------
18 
19 #include "JetScapeLogger.h"
20 #include "iSpectraSamplerWrapper.h"
21 
22 #include <memory>
23 #include <string>
24 #include <fstream>
25 
26 using namespace Jetscape;
27 
28 // Register the module with the base class
29 RegisterJetScapeModule<iSpectraSamplerWrapper>
30  iSpectraSamplerWrapper::reg("iSS");
31 
32 iSpectraSamplerWrapper::iSpectraSamplerWrapper() { SetId("iSS"); }
33 
34 iSpectraSamplerWrapper::~iSpectraSamplerWrapper() {}
35 
36 void iSpectraSamplerWrapper::InitTask() {
37 
38  JSINFO << "Initialize a particle sampler (iSS)";
39 
40  std::string input_file =
41  GetXMLElementText({"SoftParticlization", "iSS", "iSS_input_file"});
42  std::string table_path =
43  GetXMLElementText({"SoftParticlization", "iSS", "iSS_table_path"});
44  std::string particle_table_path =
45  GetXMLElementText({"SoftParticlization", "iSS",
46  "iSS_particle_table_path"});
47  std::string working_path =
48  GetXMLElementText({"SoftParticlization", "iSS", "iSS_working_path"});
49  int hydro_mode =
50  GetXMLElementInt({"SoftParticlization", "iSS", "hydro_mode"});
51  int number_of_repeated_sampling = GetXMLElementInt(
52  {"SoftParticlization", "iSS", "number_of_repeated_sampling"});
53  int flag_perform_decays = GetXMLElementInt(
54  {"SoftParticlization", "iSS", "Perform_resonance_decays"});
55  int afterburner_type = (
56  GetXMLElementInt({"SoftParticlization", "iSS", "afterburner_type"}));
57 
58  if (!boost_invariance) {
59  hydro_mode = 2;
60  }
61 
62  iSpectraSampler_ptr_ = std::unique_ptr<iSS>(
63  new iSS(working_path, table_path, particle_table_path, input_file));
64  iSpectraSampler_ptr_->paraRdr_ptr->readFromFile(input_file);
65 
66  // overwrite some parameters
67  int echoLevel = GetXMLElementInt({"vlevel"});
68  iSpectraSampler_ptr_->paraRdr_ptr->setVal("JSechoLevel", echoLevel);
69 
70  iSpectraSampler_ptr_->paraRdr_ptr->setVal("hydro_mode", hydro_mode);
71  iSpectraSampler_ptr_->paraRdr_ptr->setVal("afterburner_type",
72  afterburner_type);
73  iSpectraSampler_ptr_->paraRdr_ptr->setVal("output_samples_into_files", 0);
74  iSpectraSampler_ptr_->paraRdr_ptr->setVal("use_OSCAR_format", 0);
75  iSpectraSampler_ptr_->paraRdr_ptr->setVal("use_gzip_format", 0);
76  iSpectraSampler_ptr_->paraRdr_ptr->setVal("use_binary_format", 0);
77  iSpectraSampler_ptr_->paraRdr_ptr->setVal("store_samples_in_memory", 1);
78  iSpectraSampler_ptr_->paraRdr_ptr->setVal("number_of_repeated_sampling",
79  number_of_repeated_sampling);
80  iSpectraSampler_ptr_->paraRdr_ptr->setVal("perform_decays",
81  flag_perform_decays);
82 
83  // set default parameters
84  iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_shear", 1);
85  iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_bulk", 0);
86  iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_rhob", 0);
87  iSpectraSampler_ptr_->paraRdr_ptr->setVal("turn_on_diff", 0);
88 
89  iSpectraSampler_ptr_->paraRdr_ptr->setVal("include_deltaf_shear", 1);
90  iSpectraSampler_ptr_->paraRdr_ptr->setVal("include_deltaf_bulk", 0);
91  iSpectraSampler_ptr_->paraRdr_ptr->setVal("bulk_deltaf_kind", 1);
92  iSpectraSampler_ptr_->paraRdr_ptr->setVal("include_deltaf_diffusion", 0);
93 
94  iSpectraSampler_ptr_->paraRdr_ptr->setVal("restrict_deltaf", 0);
95  iSpectraSampler_ptr_->paraRdr_ptr->setVal("deltaf_max_ratio", 1.0);
96  iSpectraSampler_ptr_->paraRdr_ptr->setVal("f0_is_not_small", 1);
97 
98  iSpectraSampler_ptr_->paraRdr_ptr->setVal("calculate_vn", 0);
99  iSpectraSampler_ptr_->paraRdr_ptr->setVal("MC_sampling", 4);
100 
101  iSpectraSampler_ptr_->paraRdr_ptr->setVal(
102  "sample_upto_desired_particle_number", 0);
103  iSpectraSampler_ptr_->paraRdr_ptr->echo();
104 }
105 
106 void iSpectraSamplerWrapper::Exec() {
107  JSINFO << "running iSS ...";
108 
109  // generate symbolic links with music_input_file
110  std::string music_input_file_path = GetXMLElementText(
111  {"Hydro", "MUSIC", "MUSIC_input_file"});
112  std::string working_path =
113  GetXMLElementText({"SoftParticlization", "iSS", "iSS_working_path"});
114  std::string music_input = working_path + "/music_input";
115  std::ifstream inputfile(music_input.c_str());
116  if (!inputfile.good()) {
117  std::ostringstream system_command;
118  system_command << "ln -s " << music_input_file_path << " "
119  << music_input;
120  system(system_command.str().c_str());
121  }
122  inputfile.close();
123 
124  int status = iSpectraSampler_ptr_->read_in_FO_surface();
125  if (status != 0) {
126  JSWARN << "Some errors happened in reading in the hyper-surface";
127  exit(-1);
128  }
129 
130  auto random_seed = (*GetMt19937Generator())(); // get random seed
131  iSpectraSampler_ptr_->set_random_seed(random_seed);
132  VERBOSE(2) << "Random seed used for the iSS module" << random_seed;
133 
134  status = iSpectraSampler_ptr_->generate_samples();
135  if (status != 0) {
136  JSWARN << "Some errors happened in generating particle samples";
137  exit(-1);
138  }
139  PassHadronListToJetscape();
140  JSINFO << "iSS finished.";
141 }
142 
143 void iSpectraSamplerWrapper::Clear() {
144  VERBOSE(2) << "Finish the particle sampling";
145  iSpectraSampler_ptr_->clear();
146  for (unsigned i = 0; i < Hadron_list_.size(); i++) {
147  Hadron_list_.at(i).clear();
148  }
149  Hadron_list_.clear();
150 }
151 
152 void iSpectraSamplerWrapper::PassHadronListToJetscape() {
153  unsigned int nev = iSpectraSampler_ptr_->get_number_of_sampled_events();
154  VERBOSE(2) << "Passing all sampled hadrons to the JETSCAPE framework";
155  VERBOSE(4) << "number of events to pass : " << nev;
156  for (unsigned int iev = 0; iev < nev; iev++) {
157  std::vector<shared_ptr<Hadron>> hadrons;
158  unsigned int nparticles =
159  (iSpectraSampler_ptr_->get_number_of_particles(iev));
160  VERBOSE(4) << "event " << iev << ": number of particles = " << nparticles;
161  for (unsigned int ipart = 0; ipart < nparticles; ipart++) {
162  iSS_Hadron current_hadron =
163  (iSpectraSampler_ptr_->get_hadron(iev, ipart));
164  int hadron_label = 0;
165  int hadron_status = 11;
166  int hadron_id = current_hadron.pid;
167  //int hadron_id = 1; // just for testing need to be changed to the line above
168  double hadron_mass = current_hadron.mass;
169  FourVector hadron_p(current_hadron.px, current_hadron.py,
170  current_hadron.pz, current_hadron.E);
171  FourVector hadron_x(current_hadron.x, current_hadron.y, current_hadron.z,
172  current_hadron.t);
173 
174  // create a JETSCAPE Hadron
175  hadrons.push_back(make_shared<Hadron>(hadron_label, hadron_id,
176  hadron_status, hadron_p, hadron_x,
177  hadron_mass));
178  //Hadron* jetscape_hadron = new Hadron(hadron_label, hadron_id, hadron_status, hadron_p, hadron_x, hadron_mass);
179  //(*Hadron_list_)[iev]->push_back(*jetscape_hadron);
180  }
181  Hadron_list_.push_back(hadrons);
182  }
183  VERBOSE(4) << "JETSCAPE received " << Hadron_list_.size() << " events.";
184  for (unsigned int iev = 0; iev < Hadron_list_.size(); iev++) {
185  VERBOSE(4) << "In event " << iev << " JETSCAPE received "
186  << Hadron_list_.at(iev).size() << " particles.";
187  }
188 }
189 
190 void iSpectraSamplerWrapper::WriteTask(weak_ptr<JetScapeWriter> w) {
191  VERBOSE(4) << "In iSpectraSamplerWrapper::WriteTask";
192  auto f = w.lock();
193  if (!f)
194  return;
195 
196  f->WriteComment("JetScape module: " + GetId());
197  if (Hadron_list_.size() > 0) {
198  f->WriteComment("Final State Bulk Hadrons");
199  for (unsigned int j = 0; j < Hadron_list_.size(); j++) {
200  vector<shared_ptr<Hadron>> hadVec = Hadron_list_.at(j);
201  for (unsigned int i = 0; i < hadVec.size(); i++) {
202  f->WriteWhiteSpace("[" + to_string(i) + "] H");
203  f->Write(hadVec.at(i));
204  }
205  }
206  } else {
207  f->WriteComment("There are no bulk Hadrons");
208  }
209 }