evgen::MARLEYHelper Class Reference

#include "MARLEYHelper.h"

Public Member Functions
	MARLEYHelper (const fhicl::ParameterSet &pset, rndm::NuRandomService &rand_service, const std::string &generator_name)
void	reconfigure (const fhicl::ParameterSet &pset)
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
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, bool missing_ok=false)
void	load_marley_dictionaries ()

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

Detailed Description

Definition at line 45 of file MARLEYHelper.h.

Constructor & Destructor Documentation

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

Definition at line 35 of file MARLEYHelper.cxx.

References fHelperName, fMarleyGenerator, fMarleyLogStream, load_marley_dictionaries(), MF_LOG_INFO, reconfigure(), rndm::NuRandomService::registerEngine(), and seed.

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

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 84 of file MARLEYHelper.cxx.

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

Referenced by create_MCTruth().

{
  // 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 120 of file MARLEYHelper.cxx.

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

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

Definition at line 183 of file MARLEYHelper.cxx.

Referenced by 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 60 of file MARLEYHelper.h.

{ return *fMarleyGenerator; }

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

inline

Definition at line 61 of file MARLEYHelper.h.

References track.

{ return *fMarleyGenerator; }

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

protected

Definition at line 198 of file MARLEYHelper.cxx.

References find_file(), and value.

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

void evgen::MARLEYHelper::load_marley_dictionaries ( )

protected

Definition at line 255 of file MARLEYHelper.cxx.

References fHelperName, and MF_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) {
    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;
}

void evgen::MARLEYHelper::reconfigure

(

const fhicl::ParameterSet &

pset

)

Definition at line 223 of file MARLEYHelper.cxx.

References fHelperName, fMarleyGenerator, evgen::MarleyParameterSetWalker::get_json(), load_full_paths_into_json(), MF_LOG_INFO, and fhicl::ParameterSet::walk().

Referenced by MARLEYHelper().

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

Member Data Documentation

std::string evgen::MARLEYHelper::fHelperName

protected

Definition at line 78 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 75 of file MARLEYHelper.h.

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

std::stringstream evgen::MARLEYHelper::fMarleyLogStream

protected

Definition at line 82 of file MARLEYHelper.h.

Referenced by create_MCTruth(), and MARLEYHelper().

---

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

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