MARLEYHelper.cxx

Go to the documentation of this file.

// framework includes
#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft includes
#include "larsim/EventGenerator/MARLEY/MARLEYHelper.h"
#include "larsim/EventGenerator/MARLEY/MarleyParameterSetWalker.h"
#include "nurandom/RandomUtils/NuRandomService.h"
#include "nusimdata/SimulationBase/MCParticle.h"

// ROOT includes
#include "TInterpreter.h"
#include "TROOT.h"

// MARLEY includes
#include "marley/Event.hh"
#include "marley/Particle.hh"
#include "marley/RootJSONConfig.hh"

namespace {
  // We need to convert from MARLEY's energy units (MeV) to LArSoft's
  // (GeV) using this conversion factor
  constexpr double MeV_to_GeV = 1e-3;
}

//------------------------------------------------------------------------------
evgen::MARLEYHelper::MARLEYHelper(const fhicl::ParameterSet& pset,
                                  rndm::NuRandomService& rand_service,
                                  const std::string& helper_name)
  : fHelperName(helper_name)
{
  // Configure MARLEY using the FHiCL parameters
  this->reconfigure(pset);

  // Register this MARLEY generator with the NuRandomService. For simplicity,
  // we use a lambda as the seeder function (see NuRandomService.h for
  // details). This allows the SeedService to automatically re-seed MARLEY
  // whenever necessary. The user can set an explicit seed for MARLEY in the
  // FHiCL configuration using the "seed" parameter. If you need to get the
  // seed for MARLEY from the SeedService, note that we're using use the value
  // of the input variable helper_name as its generator instance name.
  rndm::NuRandomService::seed_t marley_seed = rand_service.registerEngine(
    [this](rndm::NuRandomService::EngineId const& /* unused */,
           rndm::NuRandomService::seed_t lar_seed) -> void {
      if (fMarleyGenerator && fMarleyGenerator.get()) {
        auto seed = static_cast<uint_fast64_t>(lar_seed);
        fMarleyGenerator->reseed(seed);
      }
    },
    fHelperName,
    pset,
    {"seed"});

  // Unless I'm mistaken, the call to registerEngine should seed the generator
  // with the seed from the FHiCL configuration file if one is included, but it
  // doesn't appear to do so (as of 16 Aug 2016, larsoft v06_03_00). As a
  // workaround, I manually reseed the generator (if needed) here using the
  // result of the call to registerEngine, which will be the seed from the
  // FHiCL file if one was given.
  // TODO: figure out what's going on here, and remove this workaround as
  // needed
  uint_fast64_t marley_cast_seed = static_cast<uint_fast64_t>(marley_seed);
  if (marley_cast_seed != fMarleyGenerator->get_seed()) {
    fMarleyGenerator->reseed(marley_cast_seed);
  }

  // Log initialization information from the MARLEY generator
  MF_LOG_INFO(fHelperName) << fMarleyLogStream.str();
  fMarleyLogStream = std::stringstream();

  // Do any needed setup of the MARLEY class dictionaries
  load_marley_dictionaries();
}

//------------------------------------------------------------------------------
void evgen::MARLEYHelper::add_marley_particles(simb::MCTruth& truth,
                                               const std::vector<marley::Particle*>& particles,
                                               const TLorentzVector& vtx_pos,
                                               bool track)
{
  // Loop over the vector of MARLEY particles and add simb::MCParticle
  // versions of each of them to the MCTruth object.
  for (const marley::Particle* p : particles) {
    // Treat all of these particles as primaries, which have negative
    // track IDs by convention
    int trackID = -1 * (truth.NParticles() + 1);

    int pdg = p->pdg_code();
    double mass = p->mass() * MeV_to_GeV;
    double px = p->px() * MeV_to_GeV;
    double py = p->py() * MeV_to_GeV;
    double pz = p->pz() * MeV_to_GeV;
    double E = p->total_energy() * MeV_to_GeV;
    TLorentzVector mom(px, py, pz, E);

    int status = 0; // don't track the particles in LArG4 by default
    if (track) status = 1;

    simb::MCParticle part(trackID /* trackID to use in Geant4 */,
                          pdg,
                          "MARLEY",
                          -1 /* primary particle */,
                          mass,
                          status);

    part.AddTrajectoryPoint(vtx_pos, mom);
    truth.Add(part);
  }
}

//------------------------------------------------------------------------------
simb::MCTruth evgen::MARLEYHelper::create_MCTruth(const TLorentzVector& vtx_pos,
                                                  marley::Event* marley_event)
{
  simb::MCTruth truth;

  truth.SetOrigin(simb::kSuperNovaNeutrino);

  marley::Event event = fMarleyGenerator->create_event();

  // Add the initial and final state particles to the MCTruth object.
  add_marley_particles(truth, event.get_initial_particles(), vtx_pos, false);
  add_marley_particles(truth, event.get_final_particles(), vtx_pos, true);

  // calculate a few parameters for the call to SetNeutrino
  const marley::Particle& nu = event.projectile();
  const marley::Particle& lep = event.ejectile();
  double qx = nu.px() - lep.px();
  double qy = nu.py() - lep.py();
  double qz = nu.pz() - lep.pz();
  double Enu = nu.total_energy();
  double Elep = lep.total_energy();
  double Q2 = qx * qx + qy * qy + qz * qz - std::pow(Enu - Elep, 2);

  // For definitions of Bjorken x, etc., a good reference is Mark Thomson's
  // set of slides on deep inelastic scattering (http://tinyurl.com/hcn5n6l)
  double bjorken_x = Q2 / (2 * event.target().mass() * (Enu - Elep));
  double inelasticity_y = 1. - Elep / Enu;

  // Include the initial excitation energy of the final-state nucleus when
  // calculating W (the final-state invariant mass of the hadronic system)
  // since the other parameters (X, Y) also take into account the 2-2
  // scattering reaction only.
  const marley::Particle& res = event.residue();
  double hadronic_mass_W = res.mass() + event.Ex();

  // TODO: do a more careful job of setting the parameters here
  truth.SetNeutrino(simb::kCC,                 // change when MARLEY can handle NC
                    simb::kUnknownInteraction, // not sure what the mode should be
                    simb::kUnknownInteraction, // not sure what the interaction type should be
                    marley_utils::get_nucleus_pid(18, 40), // Ar-40 PDG code
                    marley_utils::NEUTRON,                 // nucleon PDG
                    0, // MARLEY handles low enough energies that we shouldn't need HitQuark
                    hadronic_mass_W * MeV_to_GeV,
                    bjorken_x,      // dimensionless
                    inelasticity_y, // dimensionless
                    Q2 * std::pow(MeV_to_GeV, 2));

  if (marley_event) *marley_event = event;

  // Process the MARLEY logging messages (if any) captured by our
  // stringstream and forward them to the messagefacility logger
  std::string line;
  while (std::getline(fMarleyLogStream, line)) {
    MF_LOG_INFO(fHelperName) << line;
  }

  // Reset the MARLEY log stream
  fMarleyLogStream = std::stringstream();

  return truth;
}

//------------------------------------------------------------------------------
std::string evgen::MARLEYHelper::find_file(const std::string& fileName, const std::string& fileType)
{
  cet::search_path searchPath("FW_SEARCH_PATH");

  std::string fullName;
  searchPath.find_file(fileName, fullName);

  if (fullName.empty())
    throw cet::exception("MARLEYHelper")
      << "Cannot find MARLEY " << fileType << " data file '" << fileName << '\'';

  return fullName;
}

//------------------------------------------------------------------------------
void evgen::MARLEYHelper::load_full_paths_into_json(marley::JSON& json,
                                                    const std::string& key,
                                                    bool missing_ok)
{
  if (json.has_key(key)) {

    marley::JSON& value = json.at(key);

    if (value.is_array()) {
      // Replace each file name (which may appear in the FHiCL configuration
      // without a full path) with the full path found using cetlib
      for (auto& element : value.array_range()) {
        element = find_file(element.to_string(), key);
      }
    }

    else
      value = find_file(value.to_string(), key);
  }
  else if (!missing_ok)
    throw cet::exception("MARLEYHelper")
      << "Missing \"" << key << "\" key in the MARLEY parameters.";
}

//------------------------------------------------------------------------------
void evgen::MARLEYHelper::reconfigure(const fhicl::ParameterSet& pset)
{
  // Convert the FHiCL parameters into a JSON object that MARLEY can understand
  evgen::MarleyParameterSetWalker mpsw;
  pset.walk(mpsw);

  marley::JSON& json = mpsw.get_json();

  // Update the reaction and structure data file names to the full paths
  // using cetlib to search for them
  load_full_paths_into_json(json, "reactions", false);
  load_full_paths_into_json(json, "structure", true);

  // Also update the path for a neutrino source spectrum given in a ROOT
  // TFile
  if (json.has_key("source")) {
    marley::JSON& source_object = json.at("source");

    if (source_object.has_key("tfile")) { load_full_paths_into_json(source_object, "tfile"); }
  }

  // Create a new MARLEY configuration based on the JSON parameters
  MF_LOG_INFO("MARLEYHelper " + fHelperName) << "MARLEY will now use"
                                                " the JSON configuration\n"
                                             << json.dump_string() << '\n';
  marley::RootJSONConfig config(json);

  // Create a new marley::Generator object based on the current configuration
  fMarleyGenerator = std::make_unique<marley::Generator>(config.create_generator());
}

//------------------------------------------------------------------------------
void evgen::MARLEYHelper::load_marley_dictionaries()
{
  static bool already_loaded_marley_dict = false;

  if (already_loaded_marley_dict) return;

  // Current (24 July 2016) versions of ROOT 6 require runtime
  // loading of headers for custom classes in order to use
  // dictionaries correctly. If we're running ROOT 6+, do the
  // loading here, and give the user guidance if there are any
  // problems.
  //
  // This is the same technique used in the MARLEY source code
  // for the executable (src/marley.cc). If you change how this
  // code works, please sync changes with the executable as well.
  if (gROOT->GetVersionInt() >= 60000) {
    MF_LOG_INFO("MARLEYHelper " + fHelperName)
      << "ROOT 6 or greater"
      << " detected. Loading class information\nfrom headers"
      << " \"marley/Particle.hh\" and \"marley/Event.hh\"";
    TInterpreter::EErrorCode* ec = new TInterpreter::EErrorCode();
    gInterpreter->ProcessLine("#include \"marley/Particle.hh\"", ec);
    if (*ec != 0)
      throw cet::exception("MARLEYHelper " + fHelperName)
        << "Error loading MARLEY header Particle.hh. For MARLEY headers stored"
        << " in /path/to/include/marley/, please add /path/to/include"
        << " to your ROOT_INCLUDE_PATH environment variable and"
        << " try again.";
    gInterpreter->ProcessLine("#include \"marley/Event.hh\"");
    if (*ec != 0)
      throw cet::exception("MARLEYHelper")
        << "Error loading"
        << " MARLEY header Event.hh. For MARLEY headers stored in"
        << " /path/to/include/marley/, please add /path/to/include"
        << " to your ROOT_INCLUDE_PATH environment variable and"
        << " try again.";
  }

  // No further action is required for ROOT 5 because the compiled
  // dictionaries (which are linked to this algorithm) contain all of
  // the needed information
  already_loaded_marley_dict = true;
}