trento/src/trento.cxx

0001 // TRENTO: Reduced Thickness Event-by-event Nuclear Topology
0002 // Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
0003 // TRENTO3D: Three-dimensional extension of TRENTO by Weiyao Ke
0004 // MIT License
0005
0006 #include <iostream>
0007 #include <stdexcept>
0008 #include <string>
0009 #include <vector>
0010 #include <climits>
0011
0012 #include <boost/filesystem.hpp>
0013 #include <boost/filesystem/fstream.hpp>
0014 #include <boost/program_options.hpp>
0015
0016 #include "collider.h"
0017 #include "fwd_decl.h"
0018
0019 // CMake sets this definition.
0020 // Fall back to a sane default.
0021 #ifndef TRENTO_VERSION_STRING
0022 #define TRENTO_VERSION_STRING "dev"
0023 #endif
0024
0025 namespace trento {
0026
0027 namespace {
0028
0029 void print_version() {
0030   std::cout << "trento " << TRENTO_VERSION_STRING << '\n';
0031 }
0032
0033 void print_bibtex() {
0034   std::cout <<
0035     "TRENTO:\n"
0036     "@article{Moreland:2014oya,\n"
0037     "      author         = \"Moreland, J. Scott and Bernhard, Jonah E. and Bass,\n"
0038     "                        Steffen A.\",\n"
0039     "      title          = \"{Alternative ansatz to wounded nucleon and binary\n"
0040     "                        collision scaling in high-energy nuclear collisions}\",\n"
0041     "      journal        = \"Phys.Rev.\",\n"
0042     "      number         = \"1\",\n"
0043     "      volume         = \"C92\",\n"
0044     "      pages          = \"011901\",\n"
0045     "      doi            = \"10.1103/PhysRevC.92.011901\",\n"
0046     "      year           = \"2015\",\n"
0047     "      eprint         = \"1412.4708\",\n"
0048     "      archivePrefix  = \"arXiv\",\n"
0049     "      primaryClass   = \"nucl-th\",\n"
0050     "      SLACcitation   = \"%%CITATION = ARXIV:1412.4708;%%\",\n"
0051     "}\n\n";
0052
0053   std::cout <<
0054     "TRENTO3D:\n"
0055     "@article{Ke:2016jrd,\n"
0056     "      author         = \"Ke, Weiyao and Moreland, J. Scott and Bernhard, \n"
0057     "                         Jonah E. and Bass, Steffen A.\",\n"
0058     "      title          = \"{Constraints on rapidity-dependent initial conditions\n"
0059     "                        from charged particle pseudorapidity densities and\n"
0060     "                        two-particle correlations}\",\n"
0061     "      journal        = \"Phys. Rev.\",\n"
0062     "      volume         = \"C96\",\n"
0063     "      year           = \"2017\",\n"
0064     "      number         = \"4\",\n"
0065     "      pages          = \"044912\",\n"
0066     "      doi            = \"10.1103/PhysRevC.96.044912\",\n"
0067     "      eprint         = \"1610.08490\",\n"
0068     "      archivePrefix  = \"arXiv\",\n"
0069     "      primaryClass   = \"nucl-th\",\n"
0070     "      SLACcitation   = \"%%CITATION = ARXIV:1610.08490;%%\"\n";
0071 }
0072
0073 // TODO
0074 // void print_default_config() {
0075 //   std::cout << "to do\n";
0076 // }
0077
0078 }  // unnamed namespace
0079
0080 }  // namespace trento
0081
0082 int main(int argc, char* argv[]) {
0083   using namespace trento;
0084
0085   // Parse options with boost::program_options.
0086   // There are quite a few options, so let's separate them into logical groups.
0087   using OptDesc = po::options_description;
0088
0089   using VecStr = std::vector<std::string>;
0090   OptDesc main_opts{};
0091   main_opts.add_options()
0092     ("projectile", po::value<VecStr>()->required()->
0093      notifier(  // use a lambda to verify there are exactly two projectiles
0094          [](const VecStr& projectiles) {
0095            if (projectiles.size() != 2)
0096             throw po::required_option{"projectile"};
0097            }),
0098      "projectile symbols")
0099     ("number-events", po::value<int>()->default_value(1),
0100      "number of events");
0101
0102   // Make all main arguments positional.
0103   po::positional_options_description positional_opts{};
0104   positional_opts
0105     .add("projectile", 2)
0106     .add("number-events", 1);
0107
0108   using VecPath = std::vector<fs::path>;
0109   OptDesc general_opts{"general options"};
0110   general_opts.add_options()
0111     ("help,h", "show this help message and exit")
0112     ("version", "print version information and exit")
0113     ("bibtex", "print bibtex entry and exit")
0114     // ("default-config", "print a config file with default settings and exit")
0115     ("config-file,c", po::value<VecPath>()->value_name("FILE"),
0116      "configuration file\n(can be passed multiple times)");
0117
0118   OptDesc output_opts{"output options"};
0119   output_opts.add_options()
0120     ("quiet,q", po::bool_switch(),
0121      "do not print event properties to stdout")
0122     ("output,o", po::value<fs::path>()->value_name("PATH"),
0123      "HDF5 file or directory for text files")
0124     ("no-header", po::bool_switch(),
0125      "do not write headers to text files");
0126
0127   OptDesc phys_opts{"physical options"};
0128   phys_opts.add_options()
0129     ("reduced-thickness,p",
0130      po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0131      "reduced thickness parameter")
0132     ("fluctuation,k",
0133      po::value<double>()->value_name("FLOAT")->default_value(1., "1"),
0134      "gamma fluctuation shape parameter")
0135     ("nucleon-width,w",
0136      po::value<double>()->value_name("FLOAT")->default_value(.5, "0.5"),
0137      "Gaussian nucleon width [fm]")
0138     ("nucleon-min-dist,d",
0139      po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0140      "minimum nucleon-nucleon distance [fm]")
0141     ("mean-coeff,m",
0142      po::value<double>()->value_name("FLOAT")->default_value(1., "1."),
0143      "rapidity mean coefficient")
0144     ("std-coeff,s",
0145      po::value<double>()->value_name("FLOAT")->default_value(3., "3."),
0146      "rapidity std coefficient")
0147     ("skew-coeff,t",
0148      po::value<double>()->value_name("FLOAT")->default_value(0., "0."),
0149      "rapidity skew coefficient")
0150      ("skew-type,r",
0151       po::value<int>()->value_name("INT")->default_value(1, "1"),
0152       "rapidity skew type: 1: relative, 2: absolute, other: no skew")
0153     ("jacobian,j",
0154      po::value<double>()->value_name("FLOAT")->default_value(0.8, "0.8"),
0155      "<pt>/<mt> used in Jacobian")
0156     ("normalization,n",
0157      po::value<double>()->value_name("FLOAT")->default_value(1., "1"),
0158      "normalization factor");
0159
0160   OptDesc coll_opts{"collision options"};
0161   coll_opts.add_options()
0162     ("beam-energy,e",
0163      po::value<double>()->value_name("FLOAT")->default_value(2760, "2760"),
0164      "collision beam energy sqrt(s) [GeV], initializes cross section")
0165     ("cross-section,x",
0166      po::value<double>()->value_name("FLOAT")->default_value(-1, "off"),
0167      "manual inelastic nucleon-nucleon cross section sigma_NN [fm^2]")
0168     ("b-min",
0169      po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0170      "minimum impact parameter [fm]")
0171     ("b-max",
0172      po::value<double>()->value_name("FLOAT")->default_value(-1., "auto"),
0173      "maximum impact parameter [fm]")
0174     ("npart-min",
0175      po::value<int>()->value_name("INT")->default_value(0, "0"),
0176      "minimum Npart cut")
0177     ("npart-max",
0178      po::value<int>()->value_name("INT")->default_value(
0179      std::numeric_limits<int>::max(), "INT_MAX"), "maximum Npart cut")
0180     ("s-min",
0181      po::value<double>()->value_name("FLOAT")->default_value(0., "0"),
0182      "minimum entropy cut")
0183     ("s-max",
0184      po::value<double>()->value_name("FLOAT")->default_value(
0185      std::numeric_limits<double>::max(), "DOUBLE_MAX"), "maxmimum entropy cut")
0186     ("random-seed",
0187      po::value<int64_t>()->value_name("INT")->default_value(-1, "auto"),
0188      "random seed")
0189     ("ncoll,b", po::bool_switch(),
0190      "calculate # of binary collision and binary collision density");
0191
0192   OptDesc grid_opts{"grid options"};
0193   grid_opts.add_options()
0194     ("xy-max",
0195      po::value<double>()->value_name("FLOAT")->default_value(10., "10.0"),
0196      "xy max [fm]\n(transverse grid from -max to +max)")
0197     ("xy-step",
0198      po::value<double>()->value_name("FLOAT")->default_value(0.2, "0.2"),
0199      "transverse step size [fm]")
0200     ("eta-max",
0201      po::value<double>()->value_name("FLOAT")->default_value(0.0, "0.0"),
0202      "pseudorapidity max \n(eta grid from -max to +max)")
0203     ("eta-step",
0204      po::value<double>()->value_name("FLOAT")->default_value(0.5, "0.5"),
0205      "pseudorapidity step size");
0206
0207   // Make a meta-group containing all the option groups except the main
0208   // positional options (don't want the auto-generated usage info for those).
0209   OptDesc usage_opts{};
0210   usage_opts
0211     .add(general_opts)
0212     .add(output_opts)
0213     .add(phys_opts)
0214     .add(coll_opts)
0215     .add(grid_opts);
0216
0217   // Now a meta-group containing _all_ options.
0218   OptDesc all_opts{};
0219   all_opts
0220     .add(usage_opts)
0221     .add(main_opts);
0222
0223   // Will be used several times.
0224   const std::string usage_str{
0225     "usage: trento [options] projectile projectile [number-events = 1]\n"};
0226   const std::string usage_str3d{
0227     "To operate in 3D mode, make sure --eta-max is nonzero.\n"};
0228
0229   try {
0230     // Initialize a VarMap (boost::program_options::variables_map).
0231     // It will contain all configuration values.
0232     VarMap var_map{};
0233
0234     // Parse command line options.
0235     po::store(po::command_line_parser(argc, argv)
0236         .options(all_opts).positional(positional_opts).run(), var_map);
0237
0238     // Handle options that imply immediate exit.
0239     // Must do this _before_ po::notify() since that can throw exceptions.
0240     if (var_map.count("help")) {
0241       std::cout
0242         << usage_str << usage_str3d
0243         << "\n"
0244            "projectile = { p | d | Cu | Cu2 | Xe | Au | Au2 | Pb | U | U2 | U3 }\n"
0245         << usage_opts
0246         << "\n"
0247            "see the online documentation for complete usage information\n";
0248       return 0;
0249     }
0250     if (var_map.count("version")) {
0251       print_version();
0252       return 0;
0253     }
0254     if (var_map.count("bibtex")) {
0255       print_bibtex();
0256       return 0;
0257     }
0258     // if (var_map.count("default-config")) {
0259     //   print_default_config();
0260     //   return 0;
0261     // }
0262
0263     // Merge any config files.
0264     if (var_map.count("config-file")) {
0265       // Everything except general_opts.
0266       OptDesc config_file_opts{};
0267       config_file_opts
0268         .add(main_opts)
0269         .add(output_opts)
0270         .add(phys_opts)
0271         .add(grid_opts);
0272
0273       for (const auto& path : var_map["config-file"].as<VecPath>()) {
0274         if (!fs::exists(path)) {
0275           throw po::error{
0276             "configuration file '" + path.string() + "' not found"};
0277         }
0278         fs::ifstream ifs{path};
0279         po::store(po::parse_config_file(ifs, config_file_opts), var_map);
0280       }
0281     }
0282
0283     // Save all the final values into var_map.
0284     // Exceptions may occur here.
0285     po::notify(var_map);
0286
0287     // Go!
0288     Collider collider{var_map};
0289     collider.run_events();
0290   }
0291   catch (const po::required_option&) {
0292     // Handle this exception separately from others.
0293     // This occurs e.g. when the program is excuted with no arguments.
0294     std::cerr << usage_str << usage_str3d
0295               << "run 'trento --help' for more information\n";
0296     return 1;
0297   }
0298   catch (const std::exception& e) {
0299     // For all other exceptions just output the error message.
0300     std::cerr << e.what() << '\n';
0301     return 1;
0302   }
0303
0304   return 0;
0305 }