src/initialstate/TrentoInitial.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 #include <cstdlib>
0017 #include <sstream>
0018 #include <fstream>
0019 #include <boost/bind.hpp>
0020 #include <boost/tokenizer.hpp>
0021 #include <algorithm>
0022 #include <functional>
0023 #include <string>
0024 #include "JetScapeLogger.h"
0025
0026 #include "TrentoInitial.h"
0027
0028 namespace Jetscape {
0029
0030 // Register the module with the base class
0031 RegisterJetScapeModule<TrentoInitial> TrentoInitial::reg("TrentoInitial");
0032
0033 namespace {
0034 /// @brief Tokenize a string.  The tokens will be separated by each non-quoted
0035 ///        space or equal character.  Empty tokens are removed.
0036 ///
0037 /// @param input The string to tokenize.
0038 ///
0039 /// @return Vector of tokens.
0040 std::vector<std::string> tokenize(const std::string &input) {
0041   typedef boost::escaped_list_separator<char> separator_type;
0042   separator_type separator("\\",    // The escape characters.
0043                            "= ",    // The separator characters.
0044                            "\"\'"); // The quote characters.
0045
0046   // Tokenize the intput.
0047   boost::tokenizer<separator_type> tokens(input, separator);
0048
0049   // Copy non-empty tokens from the tokenizer into the result.
0050   std::vector<std::string> result;
0051   copy_if(tokens.begin(), tokens.end(), std::back_inserter(result),
0052           !boost::bind(&std::string::empty, _1));
0053   return result;
0054 }
0055 } // end namespace
0056
0057 // See header for explanation.
0058 TrentoInitial::TrentoInitial() : InitialState() { SetId("Trento"); }
0059
0060 TrentoInitial::~TrentoInitial() = default;
0061
0062 void TrentoInitial::InitTask() {
0063   JSINFO << " Initialzie TRENTo initial condition ";
0064
0065   // TRENTO OPTION DESK
0066   using namespace trento;
0067   using OptDesc = po::options_description;
0068   using VecStr = std::vector<std::string>;
0069   OptDesc main_opts{};
0070   main_opts.add_options()(
0071       "projectile",
0072       po::value<VecStr>()
0073           ->required()
0074           ->notifier( // use a lambda to verify there are exactly two projectiles
0075               [](const VecStr &projectiles) {
0076                 if (projectiles.size() != 2)
0077                   throw po::required_option{"projectile"};
0078               }),
0079       "projectile symbols")("number-events", po::value<int>()->default_value(1),
0080                             "number of events");
0081
0082   // Make all main arguments positional.
0083   po::positional_options_description positional_opts{};
0084   positional_opts.add("projectile", 2).add("number-events", 1);
0085
0086   using VecPath = std::vector<fs::path>;
0087   OptDesc general_opts{"general options"};
0088   general_opts.add_options()("help,h", "show this help message and exit")(
0089       "version", "print version information and exit")(
0090       "bibtex", "print bibtex entry and exit")
0091       // ("default-config", "print a config file with default settings and exit")
0092       ("config-file,c", po::value<VecPath>()->value_name("FILE"),
0093        "configuration file\n(can be passed multiple times)");
0094
0095   OptDesc output_opts{"output options"};
0096   output_opts.add_options()("quiet,q", po::bool_switch(),
0097                             "do not print event properties to stdout")(
0098       "output,o", po::value<fs::path>()->value_name("PATH"),
0099       "HDF5 file or directory for text files")(
0100       "no-header", po::bool_switch(), "do not write headers to text files");
0101
0102   OptDesc phys_opts{"physical options"};
0103   phys_opts.add_options()(
0104       "reduced-thickness,p",
0105       po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0106       "reduced thickness parameter")(
0107       "fluctuation,k",
0108       po::value<double>()->value_name("FLOAT")->default_value(1., "1"),
0109       "gamma fluctuation shape parameter")(
0110       "nucleon-width,w",
0111       po::value<double>()->value_name("FLOAT")->default_value(.5, "0.5"),
0112       "Gaussian nucleon width [fm]")(
0113       "nucleon-min-dist,d",
0114       po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0115       "minimum nucleon-nucleon distance [fm]")(
0116       "mean-coeff,m",
0117       po::value<double>()->value_name("FLOAT")->default_value(1., "1."),
0118       "rapidity mean coefficient")(
0119       "std-coeff,s",
0120       po::value<double>()->value_name("FLOAT")->default_value(3., "3."),
0121       "rapidity std coefficient")(
0122       "skew-coeff,t",
0123       po::value<double>()->value_name("FLOAT")->default_value(0., "0."),
0124       "rapidity skew coefficient")(
0125       "skew-type,r", po::value<int>()->value_name("INT")->default_value(1, "1"),
0126       "rapidity skew type: 1: relative, 2: absolute, other: no skew")(
0127       "jacobian,j",
0128       po::value<double>()->value_name("FLOAT")->default_value(0.8, "0.8"),
0129       "<pt>/<mt> used in Jacobian")(
0130       "normalization,n",
0131       po::value<double>()->value_name("FLOAT")->default_value(1., "1"),
0132       "normalization factor");
0133
0134   OptDesc coll_opts{"collision options"};
0135   coll_opts.add_options()(
0136       "beam-energy,e",
0137       po::value<double>()->value_name("FLOAT")->default_value(2760, "2760"),
0138       "collision beam energy sqrt(s) [GeV], initializes cross section")(
0139       "cross-section,x",
0140       po::value<double>()->value_name("FLOAT")->default_value(-1, "off"),
0141       "manual inelastic nucleon-nucleon cross section sigma_NN [fm^2]")(
0142       "b-min", po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0143       "minimum impact parameter [fm]")(
0144       "b-max",
0145       po::value<double>()->value_name("FLOAT")->default_value(-1., "auto"),
0146       "maximum impact parameter [fm]")(
0147       "npart-min", po::value<int>()->value_name("INT")->default_value(0, "0"),
0148       "minimum Npart cut")("npart-max",
0149                            po::value<int>()->value_name("INT")->default_value(
0150                                std::numeric_limits<int>::max(), "INT_MAX"),
0151                            "maximum Npart cut")(
0152       "s-min", po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0153       "minimum entropy cut")(
0154       "s-max",
0155       po::value<double>()->value_name("FLOAT")->default_value(
0156           std::numeric_limits<double>::max(), "DOUBLE_MAX"),
0157       "maxmimum entropy cut")(
0158       "random-seed",
0159       po::value<int64_t>()->value_name("INT")->default_value(-1, "auto"),
0160       "random seed")(
0161       "ncoll,b", po::bool_switch(),
0162       "calculate # of binary collision and binary collision density");
0163
0164   OptDesc grid_opts{"grid options"};
0165   grid_opts.add_options()(
0166       "xy-max",
0167       po::value<double>()->value_name("FLOAT")->default_value(10., "10.0"),
0168       "xy max [fm]\n(transverse grid from -max to +max)")(
0169       "xy-step",
0170       po::value<double>()->value_name("FLOAT")->default_value(0.2, "0.2"),
0171       "transverse step size [fm]")(
0172       "eta-max",
0173       po::value<double>()->value_name("FLOAT")->default_value(0.0, "0.0"),
0174       "pseudorapidity max \n(eta grid from -max to +max)")(
0175       "eta-step",
0176       po::value<double>()->value_name("FLOAT")->default_value(0.5, "0.5"),
0177       "pseudorapidity step size");
0178
0179   // Make a meta-group containing all the option groups except the main
0180   // positional options (don't want the auto-generated usage info for those).
0181   OptDesc usage_opts{};
0182   usage_opts.add(general_opts)
0183       .add(output_opts)
0184       .add(phys_opts)
0185       .add(coll_opts)
0186       .add(grid_opts);
0187
0188   // Now a meta-group containing _all_ options.
0189   OptDesc all_opts{};
0190   all_opts.add(usage_opts).add(main_opts);
0191
0192   // Will be used several times.
0193   const std::string usage_str{
0194       "usage: trento [options] projectile projectile [number-events = 1]\n"};
0195   const std::string usage_str3d{
0196       "To operate in 3D mode, make sure --eta-max is nonzero.\n"};
0197
0198   // NOW LETS FILL IN THE OPTION DESK
0199   auto phy_opts = GetXMLElement({"IS", "Trento", "PhysicsInputs"});
0200   auto cut_opts = GetXMLElement({"IS", "Trento", "CutInputs"});
0201   auto trans_opts = GetXMLElement({"IS", "Trento", "TransInputs"});
0202   auto longi_opts = GetXMLElement({"IS", "Trento", "LongiInputs"});
0203
0204   double xymax = GetXMax(), dxy = GetXStep();
0205   double etamax = GetZMax(), deta = GetZStep();
0206
0207   auto random_seed = (*GetMt19937Generator())();
0208   //TEMPORARY FOR TESTING
0209   //auto random_seed = 1;
0210   //TEMPORARY
0211   JSINFO << "Random seed used for Trento " << random_seed;
0212
0213   std::string proj(phy_opts->Attribute("projectile"));
0214   std::string targ(phy_opts->Attribute("target"));
0215   double sqrts = std::atof(phy_opts->Attribute("sqrts"));
0216   double cross_section = std::atof(phy_opts->Attribute("cross-section"));
0217   double normalization = std::atof(phy_opts->Attribute("normalization"));
0218
0219   int cen_low = std::atoi(cut_opts->Attribute("centrality-low"));
0220   int cen_high = std::atoi(cut_opts->Attribute("centrality-high"));
0221
0222   double p = std::atof(trans_opts->Attribute("reduced-thickness"));
0223   double k = std::atof(trans_opts->Attribute("fluctuation"));
0224   double w = std::atof(trans_opts->Attribute("nucleon-width"));
0225   double d = std::atof(trans_opts->Attribute("nucleon-min-dist"));
0226
0227   double mean = std::atof(longi_opts->Attribute("mean-coeff"));
0228   double var = std::atof(longi_opts->Attribute("std-coeff"));
0229   double skew = std::atof(longi_opts->Attribute("skew-coeff"));
0230   int skew_type = std::atof(longi_opts->Attribute("skew-type"));
0231   double J = std::atof(longi_opts->Attribute("jacobian"));
0232
0233   std::string options1 =
0234       +" --random-seed " + std::to_string(random_seed) + " --cross-section " +
0235       std::to_string(cross_section) + " --beam-energy " +
0236       std::to_string(sqrts) + " --reduced-thickness " + std::to_string(p) +
0237       " --fluctuation " + std::to_string(k) + " --nucleon-width " +
0238       std::to_string(w) + " --nucleon-min-dist " + std::to_string(d) +
0239       " --mean-coeff " + std::to_string(mean) + " --std-coeff " +
0240       std::to_string(var) + " --skew-coeff " + std::to_string(skew) +
0241       " --skew-type " + std::to_string(skew_type) + " --jacobian " +
0242       std::to_string(J) + " --quiet ";
0243   std::string options2 = " --normalization " + std::to_string(normalization) +
0244                          " --ncoll " // calcualte # of binary collision
0245                          + " --xy-max " + std::to_string(xymax) +
0246                          " --xy-step " + std::to_string(dxy) + " --eta-max " +
0247                          std::to_string(etamax) + " --eta-step " +
0248                          std::to_string(deta);
0249   // Handle centrality table, not normzlized, default grid, 2D (fast) !!!
0250   std::string cmd_basic = proj + " " + targ + " 10000 " + options1;
0251   VarMap var_map_basic{};
0252   po::store(po::command_line_parser(tokenize(cmd_basic))
0253                 .options(all_opts)
0254                 .positional(positional_opts)
0255                 .run(),
0256             var_map_basic);
0257
0258   std::string options_cut = "";
0259   if (cen_low == 0 && cen_high == 100) {
0260     JSINFO << "TRENTo Minimum Biased Mode Generates 0-100(%) of nuclear "
0261               "inelastic cross-section";
0262   } else {
0263     auto Ecut = GenCenTab(proj, targ, var_map_basic, cen_low, cen_high);
0264     double Ehigh = Ecut.first * normalization; // rescale the cut
0265     double Elow = Ecut.second * normalization; // rescale the cut
0266
0267     JSINFO << "The total energy density cut for centrality = [" << cen_low
0268            << ", " << cen_high << "] (%) is:";
0269     JSINFO << Elow << "<dE/deta(eta=0)<" << Ehigh;
0270     options_cut = " --s-max " + std::to_string(Ehigh) + " --s-min " +
0271                   std::to_string(Elow);
0272     // Set trento configuration
0273   }
0274   std::string cmd =
0275       proj + " " + targ + " 1 " + options1 + options2 + options_cut;
0276   JSINFO << cmd;
0277   VarMap var_map{};
0278   po::store(po::command_line_parser(tokenize(cmd))
0279                 .options(all_opts)
0280                 .positional(positional_opts)
0281                 .run(),
0282             var_map);
0283   TrentoGen_ = std::make_shared<trento::Collider>(var_map);
0284   SetRanges(xymax, xymax, etamax);
0285   SetSteps(dxy, dxy, deta);
0286   JSINFO << "TRENTo set";
0287 }
0288
0289 bool compare_E(trento::records r1, trento::records r2) {
0290   return r1.mult > r2.mult;
0291 }
0292
0293 std::pair<double, double> TrentoInitial::GenCenTab(std::string proj,
0294                                                    std::string targ,
0295                                                    VarMap var_map, int cL,
0296                                                    int cH) {
0297   // Terminate for nonsense
0298   if (cL < 0 || cL > 100 || cH < 0 || cH > 100 || cH < cL) {
0299     JSWARN << "Wrong centrality cuts! To be terminated.";
0300     exit(-1);
0301   }
0302   // These are all the parameters that could change the shape of centrality tables
0303   // Normalization prefactor parameter is factorized
0304   // They form a table header
0305   trento::Collider another_collider(var_map);
0306   double beamE = var_map["beam-energy"].as<double>();
0307   double xsection = var_map["cross-section"].as<double>();
0308   double pvalue = var_map["reduced-thickness"].as<double>();
0309   double fluct = var_map["fluctuation"].as<double>();
0310   double nuclw = var_map["nucleon-width"].as<double>();
0311   double dmin = var_map["nucleon-min-dist"].as<double>();
0312   char buffer[512];
0313   std::sprintf(buffer, "%s-%s-E-%1.0f-X-%1.2f-p-%1.2f-k-%1.2f-w-%1.2f-d-%1.2f",
0314                proj.c_str(), targ.c_str(), beamE, xsection, pvalue, fluct,
0315                nuclw, dmin);
0316   std::string header(buffer);
0317   JSINFO << "TRENTO centrality table header: " << header;
0318   // Create headering string hash tage for these parameter combination
0319   // Use this tag as a unique table filename for this specific parameter set
0320   std::hash<std::string> hash_function;
0321   size_t header_hash = hash_function(header);
0322   JSINFO << "Hash tag for this header: " << header_hash;
0323   // create dir incase it does not exist
0324   std::system("mkdir -p ./trento_data");
0325   char filename[512];
0326   std::sprintf(filename, "./trento_data/%zu", header_hash);
0327   // Step1: check it a table exist
0328   std::ifstream infile(filename);
0329   double Etab[101];
0330   double buff1, buff2;
0331   std::string line;
0332   if (infile.good()) {
0333     JSINFO << "The required centrality table exists. Load the table.";
0334     int i = 0;
0335     while (std::getline(infile, line)) {
0336       if (line[0] != '#') {
0337         std::istringstream iss(line);
0338         iss >> buff1 >> buff2 >> Etab[i];
0339         i++;
0340       }
0341     }
0342     infile.close();
0343   } else {
0344     JSINFO << "TRENTo is generating new centrality table for this new "
0345               "parameter set";
0346     JSINFO << "It may take 10(s) to 1(min).";
0347
0348     another_collider.run_events();
0349     // Get all records and sort according to totoal energy
0350     auto event_records = another_collider.all_records();
0351     std::sort(event_records.begin(), event_records.end(), compare_E);
0352     // write centrality table
0353     int nstep = std::ceil(event_records.size() / 100);
0354     std::ofstream fout(filename);
0355     fout << "#\tproj\ttarj\tsqrts\tx\tp\tk\tw\td\n"
0356          << "#\t" << proj << "\t" << targ << "\t" << beamE << "\t" << xsection
0357          << "\t" << pvalue << "\t" << fluct << "\t" << nuclw << "\t" << dmin
0358          << "\n"
0359          << "#\tcen_L\tcen_H\tun-normalized total density\n";
0360     Etab[0] = 1e10;
0361     for (int i = 1; i < 100; i += 1) {
0362       auto ee = event_records[i * nstep];
0363       fout << i - 1 << "\t" << i << "\t" << ee.mult << std::endl;
0364       Etab[i] = ee.mult;
0365     }
0366     auto ee = event_records.back();
0367     fout << 99 << "\t" << 100 << "\t" << ee.mult << std::endl;
0368     Etab[100] = ee.mult;
0369     fout.close();
0370   }
0371   JSINFO << "#########" << Etab[cL] << " " << Etab[cH];
0372   return std::make_pair(Etab[cL], Etab[cH]);
0373 }
0374
0375 void TrentoInitial::Exec() {
0376   JSINFO << " Exec TRENTo initial condition ";
0377   TrentoGen_->run_events();
0378
0379   JSINFO << " TRENTo event info: ";
0380   auto tmp_event = TrentoGen_->expose_event();
0381   info_.impact_parameter = TrentoGen_->all_records().back().b;
0382   info_.num_participant = tmp_event.npart();
0383   info_.num_binary_collisions = tmp_event.ncoll();
0384   info_.total_entropy = tmp_event.multiplicity();
0385   info_.ecc = tmp_event.eccentricity();
0386   info_.psi = tmp_event.participant_plane();
0387   info_.xmid =
0388       -GetXMax() + tmp_event.mass_center_index().first * tmp_event.dxy();
0389   info_.ymid =
0390       -GetYMax() + tmp_event.mass_center_index().second * tmp_event.dxy();
0391   JSINFO << "b\tnpart\tncoll\tET\t(x-com, y-com) (fm)";
0392   JSINFO << info_.impact_parameter << "\t" << info_.num_participant << "\t"
0393          << info_.num_binary_collisions << "\t" << info_.total_entropy << "\t"
0394          << "(" << info_.xmid << ", " << info_.ymid << ")";
0395
0396   JSINFO << " Load TRENTo density and ncoll density to JETSCAPE memory ";
0397   auto density_field = tmp_event.density_grid();
0398   auto ncoll_field = tmp_event.TAB_grid();
0399   JSINFO << density_field.num_elements() << " density elements";
0400   for (int i = 0; i < density_field.num_elements(); i++) {
0401     entropy_density_distribution_.push_back(density_field.data()[i]);
0402   }
0403   JSINFO << ncoll_field.num_elements() << " ncoll elements";
0404   for (int i = 0; i < ncoll_field.num_elements(); i++) {
0405     num_of_binary_collisions_.push_back(ncoll_field.data()[i]);
0406   }
0407   JSINFO << " TRENTO event generated and loaded ";
0408 }
0409
0410 void TrentoInitial::Clear() {
0411   VERBOSE(2) << " : Finish creating initial condition ";
0412   entropy_density_distribution_.clear();
0413   num_of_binary_collisions_.clear();
0414 }
0415
0416 } // end namespace Jetscape