evgen::MARLEYHelper Class Reference

#include "MARLEYHelper.h"

Classes
struct	Config
struct	Source_Config

Public Types
using	Name = fhicl::Name
using	Comment = fhicl::Comment

Public Member Functions
	MARLEYHelper (const fhicl::Table< Config > &conf, rndm::NuRandomService &rand_service, const std::string &generator_name)
	MARLEYHelper (const fhicl::ParameterSet &pset, rndm::NuRandomService &rand_service, const std::string &generator_name)
void	reconfigure (const fhicl::Table< Config > &conf)
simb::MCTruth	create_MCTruth (const TLorentzVector &vtx_pos, marley::Event *marley_event=nullptr)
marley::Generator &	get_generator ()
const marley::Generator &	get_generator () const
std::string	find_file (const std::string &fileName, const std::string &fileType)

Protected Member Functions
template<typename T >
marley::JSON	fhicl_atom_to_json (const fhicl::Atom< T > *atom)
template<typename T >
marley::JSON	fhicl_optional_atom_to_json (const fhicl::OptionalAtom< T > *opt_atom)
template<typename S >
marley::JSON	fhicl_sequence_to_json (const S *sequence)
marley::JSON	fhicl_parameter_to_json (const fhicl::detail::ParameterBase *par)
void	add_marley_particles (simb::MCTruth &truth, const std::vector< marley::Particle * > &particles, const TLorentzVector &vtx_pos, bool track)
void	load_full_paths_into_json (marley::JSON &json, const std::string &array_name)
void	load_marley_dictionaries ()

Protected Attributes
std::unique_ptr< marley::Generator >	fMarleyGenerator
std::string	fHelperName
std::stringstream	fMarleyLogStream

Detailed Description

Definition at line 48 of file MARLEYHelper.h.

Member Typedef Documentation

using evgen::MARLEYHelper::Comment = fhicl::Comment

Definition at line 53 of file MARLEYHelper.h.

using evgen::MARLEYHelper::Name = fhicl::Name

Definition at line 52 of file MARLEYHelper.h.

Constructor & Destructor Documentation

evgen::MARLEYHelper::MARLEYHelper	(	const fhicl::Table< Config > &	conf,
		rndm::NuRandomService &	rand_service,
		const std::string &	generator_name
	)

Definition at line 31 of file MARLEYHelper.cxx.

References fHelperName, fMarleyGenerator, fMarleyLogStream, fhicl::Table< T, KeysToIgnore >::get_PSet(), load_marley_dictionaries(), LOG_INFO, reconfigure(), rndm::NuRandomService::registerEngine(), and seed.

  : fHelperName(helper_name)
{
  // Configure MARLEY using the FHiCL parameters
  this->reconfigure(conf);

  // 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.get()->reseed(seed);
      }
    },
    fHelperName, conf.get_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
  LOG_INFO(fHelperName) << fMarleyLogStream.str();
  fMarleyLogStream = std::stringstream();

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

evgen::MARLEYHelper::MARLEYHelper ( const fhicl::ParameterSet &  pset, rndm::NuRandomService &  rand_service, const std::string &  generator_name  )

inline

Definition at line 247 of file MARLEYHelper.h.

References create_MCTruth(), and reconfigure().

        : MARLEYHelper(fhicl::Table<Config>(pset, {}),
        rand_service, generator_name) {}

Member Function Documentation

void evgen::MARLEYHelper::add_marley_particles ( simb::MCTruth &  truth, const std::vector< marley::Particle * > &  particles, const TLorentzVector &  vtx_pos, bool  track  )

protected

Definition at line 80 of file MARLEYHelper.cxx.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), E, simb::MCTruth::NParticles(), and part.

Referenced by create_MCTruth(), and fhicl_sequence_to_json().

{
  // 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 = nullptr
	)

Definition at line 111 of file MARLEYHelper.cxx.

References add_marley_particles(), fHelperName, fMarleyGenerator, fMarleyLogStream, simb::kCC, simb::kSuperNovaNeutrino, simb::kUnknownInteraction, LOG_INFO, simb::MCTruth::SetNeutrino(), and simb::MCTruth::SetOrigin().

Referenced by MARLEYHelper().

{
  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)) {
    LOG_INFO(fHelperName) << line;
  }

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

  return truth;
}

template<typename T >

marley::JSON evgen::MARLEYHelper::fhicl_atom_to_json ( const fhicl::Atom< T > *  atom )

inlineprotected

Definition at line 269 of file MARLEYHelper.h.

      {
        return marley::JSON(atom->operator()());
      }

template<typename T >

marley::JSON evgen::MARLEYHelper::fhicl_optional_atom_to_json ( const fhicl::OptionalAtom< T > *  opt_atom )

inlineprotected

Definition at line 275 of file MARLEYHelper.h.

References fhicl::detail::atom::value().

      {
        T value;
        if (opt_atom->operator()(value)) return marley::JSON(value);
        // Return a null JSON object if the fhicl::OptionalAtom wasn't used
        else return marley::JSON();
      }

marley::JSON evgen::MARLEYHelper::fhicl_parameter_to_json ( const fhicl::detail::ParameterBase *  par )

protected

Definition at line 289 of file MARLEYHelper.cxx.

References fhicl::is_atom(), fhicl::detail::ParameterBase::is_optional(), fhicl::is_sequence(), fhicl::is_table(), fhicl::detail::ParameterBase::key(), fhicl::detail::ParameterBase::parameter_type(), and fhicl::detail::ParameterBase::should_use().

Referenced by fhicl_sequence_to_json(), and reconfigure().

{
  // Don't bother with the rest of the analysis if we've been handed a
  // nullptr or if the parameter is disabled in the current configuration
  if (par && par->should_use()) {

    auto type = par->parameter_type();
    if (fhicl::is_atom(type)) {

      // This is an ugly hack to load FHiCL atoms into JSON objects. We have
      // to do it type-by-type using the current (12/2016) implementation of
      // the fhicl::Atom and fhicl::OptionalAtom template classes.
      if (!par->is_optional()) {
        if (auto atm = dynamic_cast<const fhicl::Atom<double>* >(par))
          return fhicl_atom_to_json<>(atm);
        else if (auto atm = dynamic_cast<const
          fhicl::Atom<std::string>* >(par))
        {
          return fhicl_atom_to_json<>(atm);
        }
        else throw cet::exception("MarleyGen") << "Failed to deduce"
          << " the type of the FHiCL atom " << par->key();
      }
      else { // optional atoms
        if (auto opt_atom = dynamic_cast<const
          fhicl::OptionalAtom<double>* >(par))
        {
          return fhicl_optional_atom_to_json<>(opt_atom);
        }
        else if (auto opt_atom = dynamic_cast<const fhicl::OptionalAtom<
          std::string>* >(par))
        {
          return fhicl_optional_atom_to_json<>(opt_atom);
        }
        else throw cet::exception("MarleyGen") << "Failed to deduce"
          << " the type of the FHiCL optional atom " << par->key();
      }
    }
    else if (fhicl::is_sequence(type)) {

      // This is another ugly hack to allow us to fill the JSON array. We
      // have to do it type-by-type using the current (12/2016)
      // implementation of the fhicl::Sequence template class.
      if (auto seq = dynamic_cast<const fhicl::Sequence<double>* >(par))
        return fhicl_sequence_to_json<>(seq);
      else if (auto seq = dynamic_cast<const
        fhicl::Sequence<double, 3>* >(par))
      {
        return fhicl_sequence_to_json<>(seq);
      }
      else if (auto seq = dynamic_cast<const
        fhicl::Sequence<std::string>* >(par))
      {
        return fhicl_sequence_to_json<>(seq);
      }
      else throw cet::exception("MarleyGen") << "Failed to deduce"
        << " the type of the FHiCL sequence " << par->key();

      return marley::JSON();
    }
    else if (fhicl::is_table(type)) {

      // Downcast the parameter pointer into a TableBase pointer so that we
      // can iterate over its members
      auto table = dynamic_cast<const fhicl::detail::TableBase*>(par);
      if (!table) throw cet::exception("MarleyGen") << "Failed dynamic_cast"
        << " of fhicl::detail::ParameterBase* to fhicl::detail::TableBase*"
        << " when fhicl_is_table() was true";

      // Create an empty JSON object
      marley::JSON object = marley::JSON::make(
        marley::JSON::DataType::Object);

      // Load it with the members of the table recursively
      for (const auto& m : table->members()) {
        if (m && m->should_use()) {
          marley::JSON temp = fhicl_parameter_to_json(m.get());
          // Skip members that are null (typically unused optional FHiCL
          // parameters)
          if (!temp.is_null()) object[ m->name() ] = temp;
        }
      }
      return object;
    }
  }

  // Return a null JSON object if something strange happened.
  return marley::JSON();
}

template<typename S >

marley::JSON evgen::MARLEYHelper::fhicl_sequence_to_json ( const S *  sequence )

inlineprotected

Definition at line 284 of file MARLEYHelper.h.

References add_marley_particles(), lar::dump::array(), fhicl_parameter_to_json(), load_full_paths_into_json(), and track.

      {
        marley::JSON array = marley::JSON::make(
          marley::JSON::DataType::Array);

        for (size_t i = 0; i < sequence->size(); ++i) {
          array.append(sequence->operator()(i));
        }
        return array;
      }

std::string evgen::MARLEYHelper::find_file	(	const std::string &	fileName,
		const std::string &	fileType
	)

Definition at line 176 of file MARLEYHelper.cxx.

Referenced by get_generator(), and load_full_paths_into_json().

{
  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;
}

marley::Generator& evgen::MARLEYHelper::get_generator ( )

inline

Definition at line 260 of file MARLEYHelper.h.

References fMarleyGenerator.

{ return *fMarleyGenerator; }

const marley::Generator& evgen::MARLEYHelper::get_generator ( ) const

inline

Definition at line 261 of file MARLEYHelper.h.

References find_file(), and fMarleyGenerator.

        { return *fMarleyGenerator; }

void evgen::MARLEYHelper::load_full_paths_into_json ( marley::JSON &  json, const std::string &  array_name  )

protected

Definition at line 193 of file MARLEYHelper.cxx.

References find_file(), and fhicl::detail::atom::value().

Referenced by fhicl_sequence_to_json(), and reconfigure().

{
  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 throw cet::exception("MARLEYHelper") << "Missing \"" << key
    << "\" key in the MARLEY parameters.";
}

void evgen::MARLEYHelper::load_marley_dictionaries ( )

protected

Definition at line 247 of file MARLEYHelper.cxx.

References fHelperName, and LOG_INFO.

Referenced by MARLEYHelper().

{

  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) {
    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;
}

void evgen::MARLEYHelper::reconfigure

(

const fhicl::Table< Config > &

conf

)

Definition at line 215 of file MARLEYHelper.cxx.

References fHelperName, fhicl_parameter_to_json(), fMarleyGenerator, load_full_paths_into_json(), and LOG_INFO.

Referenced by MARLEYHelper().

{
  // Convert the FHiCL parameters into a JSON object that MARLEY can understand
  marley::JSON json = fhicl_parameter_to_json(&conf);

  // 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");
  load_full_paths_into_json(json, "structure");

  // 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
  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());
}

Member Data Documentation

std::string evgen::MARLEYHelper::fHelperName

protected

Definition at line 309 of file MARLEYHelper.h.

Referenced by create_MCTruth(), load_marley_dictionaries(), MARLEYHelper(), and reconfigure().

std::unique_ptr<marley::Generator> evgen::MARLEYHelper::fMarleyGenerator

protected

Definition at line 306 of file MARLEYHelper.h.

Referenced by create_MCTruth(), get_generator(), MARLEYHelper(), and reconfigure().

std::stringstream evgen::MARLEYHelper::fMarleyLogStream

protected

Definition at line 313 of file MARLEYHelper.h.

Referenced by create_MCTruth(), and MARLEYHelper().

---

The documentation for this class was generated from the following files:

• larsim/v06_53_00/source/larsim/EventGenerator/MARLEY/MARLEYHelper.h
• larsim/v06_53_00/source/larsim/EventGenerator/MARLEY/MARLEYHelper.cxx