LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
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

Definition at line 53 of file MARLEYHelper.h.

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.

34  : fHelperName(helper_name)
35 {
36  // Configure MARLEY using the FHiCL parameters
37  this->reconfigure(conf);
38 
39  // Register this MARLEY generator with the NuRandomService. For simplicity,
40  // we use a lambda as the seeder function (see NuRandomService.h for
41  // details). This allows the SeedService to automatically re-seed MARLEY
42  // whenever necessary. The user can set an explicit seed for MARLEY in the
43  // FHiCL configuration using the "seed" parameter. If you need to get the
44  // seed for MARLEY from the SeedService, note that we're using use the value
45  // of the input variable helper_name as its generator instance name.
46  rndm::NuRandomService::seed_t marley_seed = rand_service.registerEngine(
47  [this](rndm::NuRandomService::EngineId const& /* unused */,
48  rndm::NuRandomService::seed_t lar_seed) -> void
49  {
50  if (fMarleyGenerator && fMarleyGenerator.get()) {
51  auto seed = static_cast<uint_fast64_t>(lar_seed);
52  fMarleyGenerator.get()->reseed(seed);
53  }
54  },
55  fHelperName, conf.get_PSet(), { "seed" }
56  );
57 
58  // Unless I'm mistaken, the call to registerEngine should seed the generator
59  // with the seed from the FHiCL configuration file if one is included, but it
60  // doesn't appear to do so (as of 16 Aug 2016, larsoft v06_03_00). As a
61  // workaround, I manually reseed the generator (if needed) here using the
62  // result of the call to registerEngine, which will be the seed from the
63  // FHiCL file if one was given.
64  // TODO: figure out what's going on here, and remove this workaround as
65  // needed
66  uint_fast64_t marley_cast_seed = static_cast<uint_fast64_t>(marley_seed);
67  if (marley_cast_seed != fMarleyGenerator->get_seed()) {
68  fMarleyGenerator->reseed(marley_cast_seed);
69  }
70 
71  // Log initialization information from the MARLEY generator
73  fMarleyLogStream = std::stringstream();
74 
75  // Do any needed setup of the MARLEY class dictionaries
77 }
#define LOG_INFO(category)
ParameterSet const & get_PSet() const
Definition: Table.h:69
art::RandomNumberGenerator::seed_t seed_t
std::stringstream fMarleyLogStream
Definition: MARLEYHelper.h:313
std::unique_ptr< marley::Generator > fMarleyGenerator
Definition: MARLEYHelper.h:306
std::string fHelperName
Definition: MARLEYHelper.h:309
long seed
Definition: chem4.cc:68
Identifier for a engine, made of module name and optional instance name.
Definition: EngineId.h:22
void reconfigure(const fhicl::Table< Config > &conf)
seed_t registerEngine(SeedMaster_t::Seeder_t seeder, std::string instance="")
Registers an existing engine with NuRandomService.
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().

250  rand_service, generator_name) {}
MARLEYHelper(const fhicl::Table< Config > &conf, rndm::NuRandomService &rand_service, const std::string &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().

83 {
84  // Loop over the vector of MARLEY particles and add simb::MCParticle
85  // versions of each of them to the MCTruth object.
86  for (const marley::Particle* p : particles) {
87  // Treat all of these particles as primaries, which have negative
88  // track IDs by convention
89  int trackID = -1*(truth.NParticles() + 1);
90 
91  int pdg = p->pdg_code();
92  double mass = p->mass() * MeV_to_GeV;
93  double px = p->px() * MeV_to_GeV;
94  double py = p->py() * MeV_to_GeV;
95  double pz = p->pz() * MeV_to_GeV;
96  double E = p->total_energy() * MeV_to_GeV;
97  TLorentzVector mom(px, py, pz, E);
98 
99  int status = 0; // don't track the particles in LArG4 by default
100  if (track) status = 1;
101 
102  simb::MCParticle part(trackID /* trackID to use in Geant4 */, pdg,
103  "MARLEY", -1 /* primary particle */, mass, status);
104 
105  part.AddTrajectoryPoint(vtx_pos, mom);
106  truth.Add(part);
107  }
108 }
Float_t E
Definition: plot.C:23
int NParticles() const
Definition: MCTruth.h:72
void Add(simb::MCParticle &part)
Definition: MCTruth.h:77
TString part[npart]
Definition: Style.C:32
Float_t track
Definition: plot.C:34
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().

113 {
114  simb::MCTruth truth;
115 
117 
118  marley::Event event = fMarleyGenerator->create_event();
119 
120  // Add the initial and final state particles to the MCTruth object.
121  add_marley_particles(truth, event.get_initial_particles(), vtx_pos, false);
122  add_marley_particles(truth, event.get_final_particles(), vtx_pos, true);
123 
124  // calculate a few parameters for the call to SetNeutrino
125  const marley::Particle& nu = event.projectile();
126  const marley::Particle& lep = event.ejectile();
127  double qx = nu.px() - lep.px();
128  double qy = nu.py() - lep.py();
129  double qz = nu.pz() - lep.pz();
130  double Enu = nu.total_energy();
131  double Elep = lep.total_energy();
132  double Q2 = qx*qx + qy*qy + qz*qz - std::pow(Enu - Elep, 2);
133 
134  // For definitions of Bjorken x, etc., a good reference is Mark Thomson's
135  // set of slides on deep inelastic scattering (http://tinyurl.com/hcn5n6l)
136  double bjorken_x = Q2 / (2 * event.target().mass() * (Enu - Elep));
137  double inelasticity_y = 1. - Elep / Enu;
138 
139  // Include the initial excitation energy of the final-state nucleus when
140  // calculating W (the final-state invariant mass of the hadronic system)
141  // since the other parameters (X, Y) also take into account the 2-2
142  // scattering reaction only.
143  const marley::Particle& res = event.residue();
144  double hadronic_mass_W = res.mass() + event.Ex();
145 
146  // TODO: do a more careful job of setting the parameters here
147  truth.SetNeutrino(
148  simb::kCC, // change when MARLEY can handle NC
149  simb::kUnknownInteraction, // not sure what the mode should be
150  simb::kUnknownInteraction, // not sure what the interaction type should be
151  marley_utils::get_nucleus_pid(18, 40), // Ar-40 PDG code
152  marley_utils::NEUTRON, // nucleon PDG
153  0, // MARLEY handles low enough energies that we shouldn't need HitQuark
154  hadronic_mass_W * MeV_to_GeV,
155  bjorken_x, // dimensionless
156  inelasticity_y, // dimensionless
157  Q2 * std::pow(MeV_to_GeV, 2)
158  );
159 
160  if (marley_event) *marley_event = event;
161 
162  // Process the MARLEY logging messages (if any) captured by our
163  // stringstream and forward them to the messagefacility logger
164  std::string line;
165  while(std::getline(fMarleyLogStream, line)) {
166  LOG_INFO(fHelperName) << line;
167  }
168 
169  // Reset the MARLEY log stream
170  fMarleyLogStream = std::stringstream();
171 
172  return truth;
173 }
#define LOG_INFO(category)
void add_marley_particles(simb::MCTruth &truth, const std::vector< marley::Particle * > &particles, const TLorentzVector &vtx_pos, bool track)
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:78
void SetNeutrino(int CCNC, int mode, int interactionType, int target, int nucleon, int quark, double w, double x, double y, double qsqr)
Definition: MCTruth.cxx:30
std::stringstream fMarleyLogStream
Definition: MARLEYHelper.h:313
std::unique_ptr< marley::Generator > fMarleyGenerator
Definition: MARLEYHelper.h:306
std::string fHelperName
Definition: MARLEYHelper.h:309
n.b.: this group is similar but not quite, entirely unlike GENIE ScatteringType convention ...
Definition: MCNeutrino.h:84
Supernova neutrinos.
Definition: MCTruth.h:23
Event generator information.
Definition: MCTruth.h:30
Event finding and building.
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.

271  {
272  return marley::JSON(atom->operator()());
273  }
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().

277  {
278  T value;
279  if (opt_atom->operator()(value)) return marley::JSON(value);
280  // Return a null JSON object if the fhicl::OptionalAtom wasn't used
281  else return marley::JSON();
282  }
std::string value(boost::any const &)
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().

291 {
292  // Don't bother with the rest of the analysis if we've been handed a
293  // nullptr or if the parameter is disabled in the current configuration
294  if (par && par->should_use()) {
295 
296  auto type = par->parameter_type();
297  if (fhicl::is_atom(type)) {
298 
299  // This is an ugly hack to load FHiCL atoms into JSON objects. We have
300  // to do it type-by-type using the current (12/2016) implementation of
301  // the fhicl::Atom and fhicl::OptionalAtom template classes.
302  if (!par->is_optional()) {
303  if (auto atm = dynamic_cast<const fhicl::Atom<double>* >(par))
304  return fhicl_atom_to_json<>(atm);
305  else if (auto atm = dynamic_cast<const
306  fhicl::Atom<std::string>* >(par))
307  {
308  return fhicl_atom_to_json<>(atm);
309  }
310  else throw cet::exception("MarleyGen") << "Failed to deduce"
311  << " the type of the FHiCL atom " << par->key();
312  }
313  else { // optional atoms
314  if (auto opt_atom = dynamic_cast<const
316  {
317  return fhicl_optional_atom_to_json<>(opt_atom);
318  }
319  else if (auto opt_atom = dynamic_cast<const fhicl::OptionalAtom<
320  std::string>* >(par))
321  {
322  return fhicl_optional_atom_to_json<>(opt_atom);
323  }
324  else throw cet::exception("MarleyGen") << "Failed to deduce"
325  << " the type of the FHiCL optional atom " << par->key();
326  }
327  }
328  else if (fhicl::is_sequence(type)) {
329 
330  // This is another ugly hack to allow us to fill the JSON array. We
331  // have to do it type-by-type using the current (12/2016)
332  // implementation of the fhicl::Sequence template class.
333  if (auto seq = dynamic_cast<const fhicl::Sequence<double>* >(par))
334  return fhicl_sequence_to_json<>(seq);
335  else if (auto seq = dynamic_cast<const
337  {
338  return fhicl_sequence_to_json<>(seq);
339  }
340  else if (auto seq = dynamic_cast<const
342  {
343  return fhicl_sequence_to_json<>(seq);
344  }
345  else throw cet::exception("MarleyGen") << "Failed to deduce"
346  << " the type of the FHiCL sequence " << par->key();
347 
348  return marley::JSON();
349  }
350  else if (fhicl::is_table(type)) {
351 
352  // Downcast the parameter pointer into a TableBase pointer so that we
353  // can iterate over its members
354  auto table = dynamic_cast<const fhicl::detail::TableBase*>(par);
355  if (!table) throw cet::exception("MarleyGen") << "Failed dynamic_cast"
356  << " of fhicl::detail::ParameterBase* to fhicl::detail::TableBase*"
357  << " when fhicl_is_table() was true";
358 
359  // Create an empty JSON object
360  marley::JSON object = marley::JSON::make(
361  marley::JSON::DataType::Object);
362 
363  // Load it with the members of the table recursively
364  for (const auto& m : table->members()) {
365  if (m && m->should_use()) {
366  marley::JSON temp = fhicl_parameter_to_json(m.get());
367  // Skip members that are null (typically unused optional FHiCL
368  // parameters)
369  if (!temp.is_null()) object[ m->name() ] = temp;
370  }
371  }
372  return object;
373  }
374  }
375 
376  // Return a null JSON object if something strange happened.
377  return marley::JSON();
378 }
bool is_atom(par_type const pt)
marley::JSON fhicl_parameter_to_json(const fhicl::detail::ParameterBase *par)
bool is_sequence(par_type const pt)
par_type parameter_type() const
Definition: ParameterBase.h:74
bool is_table(par_type const pt)
std::string key() const
Definition: ParameterBase.h:44
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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.

286  {
287  marley::JSON array = marley::JSON::make(
288  marley::JSON::DataType::Array);
289 
290  for (size_t i = 0; i < sequence->size(); ++i) {
291  array.append(sequence->operator()(i));
292  }
293  return array;
294  }
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:228
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().

178 {
179  cet::search_path searchPath("FW_SEARCH_PATH");
180 
181  std::string fullName;
182  searchPath.find_file(fileName, fullName);
183 
184  if (fullName.empty())
185  throw cet::exception("MARLEYHelper")
186  << "Cannot find MARLEY " << fileType << " data file '"
187  << fileName << '\'';
188 
189  return fullName;
190 }
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
marley::Generator& evgen::MARLEYHelper::get_generator ( )
inline

Definition at line 260 of file MARLEYHelper.h.

References fMarleyGenerator.

260 { return *fMarleyGenerator; }
std::unique_ptr< marley::Generator > fMarleyGenerator
Definition: MARLEYHelper.h:306
const marley::Generator& evgen::MARLEYHelper::get_generator ( ) const
inline

Definition at line 261 of file MARLEYHelper.h.

References find_file(), and fMarleyGenerator.

262  { return *fMarleyGenerator; }
std::unique_ptr< marley::Generator > fMarleyGenerator
Definition: MARLEYHelper.h:306
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().

195 {
196  if ( json.has_key(key) ) {
197 
198  marley::JSON& value = json.at(key);
199 
200  if ( value.is_array() ) {
201  // Replace each file name (which may appear in the FHiCL configuration
202  // without a full path) with the full path found using cetlib
203  for (auto& element : value.array_range()) {
204  element = find_file(element.to_string(), key);
205  }
206  }
207 
208  else value = find_file(value.to_string(), key);
209  }
210  else throw cet::exception("MARLEYHelper") << "Missing \"" << key
211  << "\" key in the MARLEY parameters.";
212 }
std::string find_file(const std::string &fileName, const std::string &fileType)
std::string value(boost::any const &)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgen::MARLEYHelper::load_marley_dictionaries ( )
protected

Definition at line 247 of file MARLEYHelper.cxx.

References fHelperName, and LOG_INFO.

Referenced by MARLEYHelper().

248 {
249 
250  static bool already_loaded_marley_dict = false;
251 
252  if (already_loaded_marley_dict) return;
253 
254  // Current (24 July 2016) versions of ROOT 6 require runtime
255  // loading of headers for custom classes in order to use
256  // dictionaries correctly. If we're running ROOT 6+, do the
257  // loading here, and give the user guidance if there are any
258  // problems.
259  //
260  // This is the same technique used in the MARLEY source code
261  // for the executable (src/marley.cc). If you change how this
262  // code works, please sync changes with the executable as well.
263  if (gROOT->GetVersionInt() >= 60000) {
264  LOG_INFO("MARLEYHelper " + fHelperName) << "ROOT 6 or greater"
265  << " detected. Loading class information\nfrom headers"
266  << " \"marley/Particle.hh\" and \"marley/Event.hh\"";
267  TInterpreter::EErrorCode* ec = new TInterpreter::EErrorCode();
268  gInterpreter->ProcessLine("#include \"marley/Particle.hh\"", ec);
269  if (*ec != 0) throw cet::exception("MARLEYHelper " + fHelperName)
270  << "Error loading MARLEY header Particle.hh. For MARLEY headers stored"
271  << " in /path/to/include/marley/, please add /path/to/include"
272  << " to your ROOT_INCLUDE_PATH environment variable and"
273  << " try again.";
274  gInterpreter->ProcessLine("#include \"marley/Event.hh\"");
275  if (*ec != 0) throw cet::exception("MARLEYHelper") << "Error loading"
276  << " MARLEY header Event.hh. For MARLEY headers stored in"
277  << " /path/to/include/marley/, please add /path/to/include"
278  << " to your ROOT_INCLUDE_PATH environment variable and"
279  << " try again.";
280  }
281 
282  // No further action is required for ROOT 5 because the compiled
283  // dictionaries (which are linked to this algorithm) contain all of
284  // the needed information
285  already_loaded_marley_dict = true;
286 }
#define LOG_INFO(category)
std::string fHelperName
Definition: MARLEYHelper.h:309
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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().

217 {
218  // Convert the FHiCL parameters into a JSON object that MARLEY can understand
219  marley::JSON json = fhicl_parameter_to_json(&conf);
220 
221  // Update the reaction and structure data file names to the full paths
222  // using cetlib to search for them
223  load_full_paths_into_json(json, "reactions");
224  load_full_paths_into_json(json, "structure");
225 
226  // Also update the path for a neutrino source spectrum given in a ROOT
227  // TFile
228  if ( json.has_key("source") ) {
229  marley::JSON& source_object = json.at("source");
230 
231  if ( source_object.has_key("tfile") ) {
232  load_full_paths_into_json(source_object, "tfile");
233  }
234  }
235 
236  // Create a new MARLEY configuration based on the JSON parameters
237  LOG_INFO("MARLEYHelper " + fHelperName) << "MARLEY will now use"
238  " the JSON configuration\n" << json.dump_string() << '\n';
239  marley::RootJSONConfig config(json);
240 
241  // Create a new marley::Generator object based on the current configuration
242  fMarleyGenerator = std::make_unique<marley::Generator>(
243  config.create_generator());
244 }
#define LOG_INFO(category)
marley::JSON fhicl_parameter_to_json(const fhicl::detail::ParameterBase *par)
std::unique_ptr< marley::Generator > fMarleyGenerator
Definition: MARLEYHelper.h:306
std::string fHelperName
Definition: MARLEYHelper.h:309
void load_full_paths_into_json(marley::JSON &json, const std::string &array_name)

Member Data Documentation

std::string evgen::MARLEYHelper::fHelperName
protected
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: