rivetoutput.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2021 Christian Holm Christensen
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#include "smash/rivetoutput.h"
 
#include <Rivet/Rivet.hh>
#include <Rivet/Tools/Logging.hh>
 
#include "smash/logging.h"
 
namespace smash {
RivetOutput::RivetOutput(const bf::path& path, std::string name,
                         const bool full_event, const int total_N,
                         const int proj_N, const OutputParameters& out_par)
    : HepMcInterface(name, full_event, total_N, proj_N),
      handler_(),
      filename_(path / (name + ".yoda")),
      need_init_(true),
      rivet_confs_(out_par.subcon_for_rivet) {
  handler_ = std::make_shared<Rivet::AnalysisHandler>();
  setup();
}
 
RivetOutput::~RivetOutput() {
  logg[LOutput].debug() << "Writing Rivet results to " << filename_
                        << std::endl;
  analysis_handler_proxy()->finalize();
  analysis_handler_proxy()->writeData(filename_.string());
}
 
void RivetOutput::at_eventend(const Particles& particles,
                              const int32_t event_number,
                              const EventInfo& event) {
  HepMcInterface::at_eventend(particles, event_number, event);
 
  // Initialize Rivet on first event
  if (need_init_) {
    logg[LOutput].debug() << "Initialising Rivet" << std::endl;
    need_init_ = false;
    analysis_handler_proxy()->init(event_);
  }
 
  logg[LOutput].debug() << "Analysing event " << event_number << std::endl;
  // Let Rivet analyse the event
  analysis_handler_proxy()->analyze(event_);
}
 
void RivetOutput::add_analysis(const std::string& name) {
  analysis_handler_proxy()->addAnalysis(name);
}
 
void RivetOutput::add_path(const std::string& path) {
  Rivet::addAnalysisLibPath(path);
  Rivet::addAnalysisDataPath(path);
}
 
void RivetOutput::add_preload(const std::string& file) {
  analysis_handler_proxy()->readData(file);
}
 
void RivetOutput::set_ignore_beams(bool ignore) {
  logg[LOutput].info() << "Ignore beams? " << (ignore ? "yes" : "no")
                       << std::endl;
  analysis_handler_proxy()->setIgnoreBeams(ignore);
}
 
void RivetOutput::set_log_level(const std::string& name,
                                const std::string& level) {
  std::string fname(name);
  if (fname.rfind("Rivet", 0) != 0) {
    fname = "Rivet." + fname;
  }
 
  auto upcase = [](const std::string& s) {
    std::string out(s);
    std::transform(out.begin(), out.end(), out.begin(),
                   [](char c) { return std::toupper(c); });
    return out;
  };
 
  try {
    Rivet::Log::setLevel(fname, Rivet::Log::getLevelFromName(upcase(level)));
  } catch (...) {
  }
}
 
void RivetOutput::set_cross_section(double xs, double xserr) {
  analysis_handler_proxy()->setCrossSection(xs, xserr, true);
}
 
void RivetOutput::setup() {
  logg[LOutput].debug() << "Setting up from configuration:\n"
                        << rivet_confs_.to_string() << std::endl;
 
  // Paths to analyses libraries and data
  if (rivet_confs_.has_value({"Paths"})) {
    logg[LOutput].info() << "Processing paths" << std::endl;
    std::vector<std::string> path = rivet_confs_.take({"Paths"});
    for (auto p : path)
      add_path(p);
  }
 
  // Data files to pre-load e.g., for centrality configurations
  if (rivet_confs_.has_value({"Preloads"})) {
    logg[LOutput].info() << "Processing preloads" << std::endl;
    std::vector<std::string> prel = rivet_confs_.take({"Preloads"});
    for (auto p : prel)
      add_preload(p);
  }
 
  // Analyses (including options) to add to run
  if (rivet_confs_.has_value({"Analyses"})) {
    logg[LOutput].info() << "Processing analyses" << std::endl;
    std::vector<std::string> anas = rivet_confs_.take({"Analyses"});
    for (auto p : anas)
      add_analysis(p);
  }
 
  // Whether Rivet should ignore beams
  if (rivet_confs_.has_value({"Ignore_Beams"})) {
    set_ignore_beams(rivet_confs_.take({"Ignore_Beams"}));
  } else {
    // we must explicity tell Rivet, through the handler, to ignore beam checks
    set_ignore_beams(true);
  }
 
  // Cross sections
  if (rivet_confs_.has_value({"Cross_Section"})) {
    std::array<double, 2> xs = rivet_confs_.take({"Cross_Section"});
    set_cross_section(xs[0], xs[1]);
  }
 
  // Logging in Rivet
  if (rivet_confs_.has_value({"Logging"})) {
    std::map<std::string, std::string> logs = rivet_confs_.take({"Logging"});
    for (auto nl : logs)
      set_log_level(nl.first, nl.second);
  }
 
  // Treatment of event weights in Rivet
  if (rivet_confs_.has_value({"Weights"})) {
    auto wconf = rivet_confs_["Weights"];
 
    // Do not care about multi weights - bool
    if (wconf.has_value({"No_Multi"})) {
      analysis_handler_proxy()->skipMultiWeights(wconf.take({"No_Multi"}));
    }
 
    // Set nominal weight name
    if (wconf.has_value({"Nominal"})) {
      analysis_handler_proxy()->setNominalWeightName(wconf.take({"Nominal"}));
    }
 
    // Set cap (maximum) on weights
    if (wconf.has_value({"Cap"})) {
      analysis_handler_proxy()->setWeightCap(wconf.take({"Cap"}));
    }
 
    // Whether to smear for NLO calculations
    if (wconf.has_value({"NLO_Smearing"})) {
      analysis_handler_proxy()->setNLOSmearing(wconf.take({"NLO_Smearing"}));
    }
 
    // Select which weights to enable
    if (wconf.has_value({"Select"})) {
      std::vector<std::string> sel = wconf.take({"Select"});
      std::stringstream s;
      int comma = 0;
      for (auto w : sel)
        s << (comma++ ? "," : "") << w;
      analysis_handler_proxy()->selectMultiWeights(s.str());
    }
 
    // Select weights to disable
    if (wconf.has_value({"Deselect"})) {
      std::vector<std::string> sel = wconf.take({"Deselect"});
      std::stringstream s;
      int comma = 0;
      for (auto w : sel)
        s << (comma++ ? "," : "") << w;
      analysis_handler_proxy()->deselectMultiWeights(s.str());
    }
  }
  logg[LOutput].debug() << "After processing configuration:\n"
                        << rivet_confs_.to_string() << std::endl;
}
}  // namespace smash