evgen::MarleyTimeGen Class Reference

Inheritance diagram for evgen::MarleyTimeGen:

Classes
struct	Config
	Collection of configuration parameters for the module. More...
class	FitParameters
	Stores parsed fit parameters from a single time bin and neutrino type in a "fit"-format spectrum file. More...
class	TimeFit
	Stores fitting parameters for all neutrino types from a single time bin in a "fit"-format spectrum file. More...

Public Types
using	Name = fhicl::Name
using	Comment = fhicl::Comment
using	Parameters = art::EDProducer::Table< Config >
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = ProducerBase::Table< UserConfig, KeysToIgnore >

Public Member Functions
	MarleyTimeGen (const Parameters &p)
virtual	~MarleyTimeGen ()
virtual void	produce (art::Event &e) override
virtual void	beginRun (art::Run &run) override
virtual void	reconfigure (const Parameters &p)
template<typename PROD , BranchType B = InEvent>
ProductID	getProductID (std::string const &instanceName={}) const
template<typename PROD , BranchType B>
ProductID	getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
bool	modifiesEvent () const
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename T , art::BranchType BT>
art::ProductToken< T >	consumes (InputTag const &it)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , art::BranchType BT>
art::ViewToken< T >	consumesView (InputTag const &it)
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename T , art::BranchType BT>
art::ProductToken< T >	mayConsume (InputTag const &it)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , art::BranchType BT>
art::ViewToken< T >	mayConsumeView (InputTag const &it)
base_engine_t &	createEngine (seed_t seed)
base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make)
base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
seed_t	get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)

Static Public Member Functions
static cet::exempt_ptr< Consumer >	non_module_context ()

Protected Types
enum	TimeGenSamplingMode { TimeGenSamplingMode::HISTOGRAM, TimeGenSamplingMode::UNIFORM_TIME, TimeGenSamplingMode::UNIFORM_ENERGY }
	Enumerated type that defines the allowed ways that a neutrino's energy and arrival time may be sampled. More...
enum	PinchParamType { PinchParamType::ALPHA, PinchParamType::BETA }
	Enumerated type that defines the pinching parameter conventions that are understood by this module. More...
enum	SpectrumFileFormat { SpectrumFileFormat::RootTH2D, SpectrumFileFormat::FIT }
	Enumerated type that defines the allowed neutrino spectrum input file formats. More...

Protected Member Functions
simb::MCTruth	make_uniform_energy_mctruth (double E_min, double E_max, double &E_nu, const TLorentzVector &vertex_pos)
	Creates a simb::MCTruth object using a uniformly-sampled neutrino energy. More...
std::unique_ptr< marley::NeutrinoSource >	source_from_time_fit (const TimeFit &fit)
	Create a MARLEY neutrino source object using a set of fit parameters for a particular time bin. More...
void	create_truths_th2d (simb::MCTruth &mc_truth, sim::SupernovaTruth &sn_truth, const TLorentzVector &vertex_pos)
	Create simb::MCTruth and sim::SupernovaTruth objects using spectrum information from a ROOT TH2D. More...
void	create_truths_time_fit (simb::MCTruth &mc_truth, sim::SupernovaTruth &sn_truth, const TLorentzVector &vertex_pos)
	Create simb::MCTruth and sim::SupernovaTruth objects using a neutrino spectrum described by a previously-parsed "fit"-format file. More...
void	make_final_timefit (double time)
	Helper function that makes a final dummy TimeFit object so that the final real time bin can have a right edge. More...
void	make_nu_emission_histograms () const
	Makes ROOT histograms showing the emitted neutrinos in each time bin when using a "fit"-format spectrum file. More...
CurrentProcessingContext const *	currentContext () const
void	validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
void	prepareForJob (fhicl::ParameterSet const &pset)
void	showMissingConsumes () const

Protected Attributes
std::unique_ptr< evgen::MARLEYHelper >	fMarleyHelper
	Object that provides an interface to the MARLEY event generator. More...
std::unique_ptr< evgen::ActiveVolumeVertexSampler >	fVertexSampler
	Algorithm that allows us to sample vertex locations within the active volume(s) of the detector. More...
std::unique_ptr< marley::Event >	fEvent
	unique_ptr to the current event created by MARLEY More...
std::unique_ptr< TH2D >	fSpectrumHist
	ROOT TH2D that contains the time-dependent spectrum to use when sampling neutrino times and energies. More...
std::vector< TimeFit >	fTimeFits
	Vector that contains the fit parameter information for each time bin when using a "fit"-format spectrum file. More...
TimeGenSamplingMode	fSamplingMode
	Represents the sampling mode to use when selecting neutrino times and energies. More...
PinchParamType	fPinchType
	The pinching parameter convention to use when interpreting the time-dependent fits. More...
SpectrumFileFormat	fSpectrumFileFormat
	Format to assume for the neutrino spectrum input file. More...
TTree *	fEventTree
	The event TTree created by MARLEY. More...
uint_fast32_t	fRunNumber
	Run number from the art::Event being processed. More...
uint_fast32_t	fSubRunNumber
	Subrun number from the art::Event being processed. More...
uint_fast32_t	fEventNumber
	Event number from the art::Event being processed. More...
double	fTNu
	Time since the supernova core bounce for the current MARLEY neutrino vertex. More...
double	fWeight
	Statistical weight for the current MARLEY neutrino vertex. More...
double	fFluxAveragedCrossSection
	Flux-averaged total cross section (fm2, average is taken over all energies and times for all defined reactions) used by MARLEY to generate neutrino vertices. More...
unsigned int	fNeutrinosPerEvent
	The number of MARLEY neutrino vertices to generate in each art::Event. More...
double	fFitEmin
	Minimum neutrino energy to consider when using a "fit"-format spectrum file. More...
double	fFitEmax
	Maximum neutrino energy to consider when using a "fit"-format spectrum file. More...

Detailed Description

Definition at line 115 of file MARLEYTimeGen_module.cc.

Member Typedef Documentation

using evgen::MarleyTimeGen::Comment = fhicl::Comment

Definition at line 120 of file MARLEYTimeGen_module.cc.

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 34 of file EDProducer.h.

using evgen::MarleyTimeGen::Name = fhicl::Name

Definition at line 119 of file MARLEYTimeGen_module.cc.

using evgen::MarleyTimeGen::Parameters = art::EDProducer::Table<Config>

Definition at line 227 of file MARLEYTimeGen_module.cc.

template<typename UserConfig , typename KeysToIgnore = void>

using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>

inherited

Definition at line 35 of file EDProducer.h.

Member Enumeration Documentation

enum evgen::MarleyTimeGen::PinchParamType

strongprotected

Enumerated type that defines the pinching parameter conventions that are understood by this module.

Enumerator
ALPHA
BETA

Definition at line 466 of file MARLEYTimeGen_module.cc.

{ ALPHA, BETA };

enum evgen::MarleyTimeGen::SpectrumFileFormat

strongprotected

Enumerated type that defines the allowed neutrino spectrum input file formats.

Enumerator
RootTH2D
FIT

Definition at line 487 of file MARLEYTimeGen_module.cc.

{ RootTH2D, FIT };

enum evgen::MarleyTimeGen::TimeGenSamplingMode

strongprotected

Enumerated type that defines the allowed ways that a neutrino's energy and arrival time may be sampled.

Enumerator
HISTOGRAM
UNIFORM_TIME
UNIFORM_ENERGY

Definition at line 435 of file MARLEYTimeGen_module.cc.

{ HISTOGRAM, UNIFORM_TIME, UNIFORM_ENERGY };

Constructor & Destructor Documentation

evgen::MarleyTimeGen::MarleyTimeGen ( const Parameters &  p )

explicit

Definition at line 531 of file MARLEYTimeGen_module.cc.

References fEvent, fEventNumber, fEventTree, fRunNumber, fSubRunNumber, fTNu, fWeight, art::TFileDirectory::make(), and reconfigure().

  : fEvent(new marley::Event), fRunNumber(0), fSubRunNumber(0),
  fEventNumber(0), fTNu(0.), fFluxAveragedCrossSection(0.)
{
  // Configure the module (including MARLEY itself) using the FHiCL parameters
  this->reconfigure(p);

  // Create a ROOT TTree using the TFileService that will store the MARLEY
  // event objects (useful for debugging purposes)
  art::ServiceHandle<art::TFileService> tfs;
  fEventTree = tfs->make<TTree>("MARLEY_event_tree",
    "Neutrino events generated by MARLEY");
  fEventTree->Branch("event", "marley::Event", fEvent.get());

  // Add branches that give the art::Event run, subrun, and event numbers for
  // easy match-ups between the MARLEY and art TTrees. All three are recorded
  // as 32-bit unsigned integers.
  fEventTree->Branch("run_number", &fRunNumber, "run_number/i");
  fEventTree->Branch("subrun_number", &fSubRunNumber, "subrun_number/i");
  fEventTree->Branch("event_number", &fEventNumber, "event_number/i");
  fEventTree->Branch("tSN", &fTNu, "tSN/D");
  fEventTree->Branch("weight", &fWeight, "weight/D");

  produces< std::vector<simb::MCTruth> >();
  produces< std::vector<sim::SupernovaTruth> >();
  produces< art::Assns<simb::MCTruth, sim::SupernovaTruth> >();
  produces< sumdata::RunData, art::InRun >();
}

evgen::MarleyTimeGen::~MarleyTimeGen ( )

virtual

Definition at line 561 of file MARLEYTimeGen_module.cc.

{
}

Member Function Documentation

void evgen::MarleyTimeGen::beginRun ( art::Run &  run )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 566 of file MARLEYTimeGen_module.cc.

References geo::GeometryCore::DetectorName(), and art::Run::put().

                                             {
  // grab the geometry object to see what geometry we are using
  art::ServiceHandle<geo::Geometry> geo;
  std::unique_ptr<sumdata::RunData>
    runcol(new sumdata::RunData(geo->DetectorName()));

  run.put(std::move(runcol));
}

template<typename T , BranchType = InEvent>

ProductToken<T> art::Consumer::consumes ( InputTag const &  )

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it )

inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

{
  if (!moduleContext_)
    return ProductToken<T>::invalid();

  consumables_[BT].emplace_back(ConsumableType::Product,
                                TypeID{typeid(T)},
                                it.label(),
                                it.instance(),
                                it.process());
  return ProductToken<T>{it};
}

template<typename T , art::BranchType BT>

void art::Consumer::consumesMany ( )

inherited

Definition at line 162 of file Consumer.h.

{
  if (!moduleContext_)
    return;

  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
}

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it )

inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

{
  if (!moduleContext_)
    return ViewToken<T>::invalid();

  consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                TypeID{typeid(T)},
                                it.label(),
                                it.instance(),
                                it.process());
  return ViewToken<T>{it};
}

void evgen::MarleyTimeGen::create_truths_th2d ( simb::MCTruth &  mc_truth, sim::SupernovaTruth &  sn_truth, const TLorentzVector &  vertex_pos  )

protected

Create simb::MCTruth and sim::SupernovaTruth objects using spectrum information from a ROOT TH2D.

Definition at line 576 of file MARLEYTimeGen_module.cc.

References fEvent, fFluxAveragedCrossSection, fMarleyHelper, fSamplingMode, fSpectrumHist, fTNu, fWeight, HISTOGRAM, sim::kUnbiased, sim::kUniformEnergy, sim::kUniformTime, make_uniform_energy_mctruth(), UNIFORM_ENERGY, and UNIFORM_TIME.

Referenced by produce().

{
  // Get a reference to the generator object created by MARLEY (we'll need
  // to do a few fancy things with it other than just creating events)
  marley::Generator& gen = fMarleyHelper->get_generator();

  if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM)
  {
    // Generate a MARLEY event using the time-integrated spectrum
    // (the generator was already configured to use it by reconfigure())
    mc_truth = fMarleyHelper->create_MCTruth(vertex_pos,
      fEvent.get());

    // Find the time distribution corresponding to the selected energy bin
    double E_nu = fEvent->projectile().total_energy();
    int E_bin_index = fSpectrumHist->GetYaxis()->FindBin(E_nu);
    TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist", E_bin_index,
      E_bin_index);
    double* time_bin_weights = t_hist->GetArray();

    // Sample a time bin from the distribution
    std::discrete_distribution<int> time_dist;
    std::discrete_distribution<int>::param_type time_params(time_bin_weights,
      &(time_bin_weights[t_hist->GetNbinsX() + 1]));
    int time_bin_index = gen.sample_from_distribution(time_dist, time_params);

    // Sample a time uniformly from within the selected time bin
    double t_min = t_hist->GetBinLowEdge(time_bin_index);
    double t_max = t_min + t_hist->GetBinWidth(time_bin_index);
    // sample a time on [ t_min, t_max )
    fTNu = gen.uniform_random_double(t_min, t_max, false);
    // Unbiased sampling was used, so assign this neutrino vertex a
    // unit statistical weight
    fWeight = ONE;

    sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection,
      sim::kUnbiased);
  }

  else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME)
  {
    // Generate a MARLEY event using the time-integrated spectrum
    // (the generator was already configured to use it by reconfigure())
    mc_truth = fMarleyHelper->create_MCTruth(vertex_pos,
      fEvent.get());

    // Sample a time uniformly
    TAxis* time_axis = fSpectrumHist->GetXaxis();
    // underflow bin has index zero
    double t_min = time_axis->GetBinLowEdge(1);
    double t_max = time_axis->GetBinLowEdge(time_axis->GetNbins() + 1);
    // sample a time on [ t_min, t_max )
    fTNu = gen.uniform_random_double(t_min, t_max, false);

    // Get the value of the true dependent probability density (probability
    // of the sampled time given the sampled energy) to use as a biasing
    // correction in the neutrino vertex weight.
    double E_nu = fEvent->projectile().total_energy();
    int E_bin_index = fSpectrumHist->GetYaxis()->FindBin(E_nu);
    TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist", E_bin_index,
      E_bin_index);
    int t_bin_index = t_hist->FindBin(fTNu);
    double weight_bias = t_hist->GetBinContent(t_bin_index) * (t_max - t_min)
      / ( t_hist->Integral() * t_hist->GetBinWidth(t_bin_index) );

    fWeight = weight_bias;

    sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection,
      sim::kUniformTime);
  }

  else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_ENERGY)
  {
    // Select a time bin using the energy-integrated spectrum
    TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist");
    double* time_bin_weights = t_hist->GetArray();

    // Sample a time bin from the distribution
    std::discrete_distribution<int> time_dist;
    std::discrete_distribution<int>::param_type time_params(time_bin_weights,
      &(time_bin_weights[t_hist->GetNbinsX() + 1]));
    int time_bin_index = gen.sample_from_distribution(time_dist, time_params);

    // Sample a time uniformly from within the selected time bin
    double t_min = t_hist->GetBinLowEdge(time_bin_index);
    double t_max = t_min + t_hist->GetBinWidth(time_bin_index);
    // sample a time on [ t_min, t_max )
    fTNu = gen.uniform_random_double(t_min, t_max, false);

    // Sample an energy uniformly over the entire allowed range
    // underflow bin has index zero
    TAxis* energy_axis = fSpectrumHist->GetYaxis();
    double E_min = energy_axis->GetBinLowEdge(1);
    double E_max = energy_axis->GetBinLowEdge(energy_axis->GetNbins() + 1);

    double E_nu = std::numeric_limits<double>::lowest();

    // Generate a MARLEY event using a uniformly sampled energy
    mc_truth = make_uniform_energy_mctruth(E_min, E_max, E_nu, vertex_pos);

    // Get the value of the true dependent probability density (probability
    // of the sampled energy given the sampled time) to use as a biasing
    // correction in the neutrino vertex weight.
    //
    // Get a 1D projection of the energy spectrum for the sampled time bin
    TH1D* energy_spect = fSpectrumHist->ProjectionY("energy_spect",
      time_bin_index, time_bin_index);

    // Create a new MARLEY neutrino source object using this projection (this
    // will create a normalized probability density that we can use) and load
    // it into the generator.
    auto nu_source = marley_root::make_root_neutrino_source(
      marley_utils::ELECTRON_NEUTRINO, energy_spect);
    double new_source_E_min = nu_source->get_Emin();
    double new_source_E_max = nu_source->get_Emax();
    gen.set_source(std::move(nu_source));
    // NOTE: The marley::Generator object normalizes the E_pdf to unity
    // automatically, but just in case, we redo it here.
    double E_pdf_integ = integrate([&gen](double E_nu)
      -> double { return gen.E_pdf(E_nu); }, new_source_E_min,
      new_source_E_max);

    // Compute the likelihood ratio that we need to bias the neutrino vertex
    // weight
    double weight_bias = (gen.E_pdf(E_nu) / E_pdf_integ) * (E_max - E_min);

    fWeight = weight_bias;

    sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection,
      sim::kUniformEnergy);
  }

  else {
    throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
      << " encountered in evgen::MarleyTimeGen::produce()";
  }

}

void evgen::MarleyTimeGen::create_truths_time_fit ( simb::MCTruth &  mc_truth, sim::SupernovaTruth &  sn_truth, const TLorentzVector &  vertex_pos  )

protected

Create simb::MCTruth and sim::SupernovaTruth objects using a neutrino spectrum described by a previously-parsed "fit"-format file.

Definition at line 1054 of file MARLEYTimeGen_module.cc.

References fEvent, fFitEmax, fFitEmin, fFluxAveragedCrossSection, fMarleyHelper, fSamplingMode, fTimeFits, fTNu, fWeight, HISTOGRAM, sim::kUnbiased, sim::kUniformEnergy, sim::kUniformTime, evgen::MarleyTimeGen::TimeFit::make_FitParameters_iterator(), evgen::MarleyTimeGen::FitParameters::make_luminosity_iterator(), make_uniform_energy_mctruth(), max, source_from_time_fit(), UNIFORM_ENERGY, and UNIFORM_TIME.

Referenced by produce().

{
  // Get a reference to the generator object created by MARLEY (we'll need
  // to do a few fancy things with it other than just creating events)
  marley::Generator& gen = fMarleyHelper->get_generator();

  // Initialize the time bin index to something absurdly large. This will help
  // us detect strange bugs that arise when it is sampled incorrectly.
  size_t time_bin_index = std::numeric_limits<size_t>::max();

  // Create an object to represent the discrete time distribution given in
  // the spectrum file. Use the luminosity for each bin as its sampling weight.
  // This distribution will actually be used to sample a neutrino arrival time
  // unless we're using the uniform time sampling mode, in which case it will
  // be used to help calculate the neutrino vertex weight.
  int source_pdg_code = gen.get_source().get_pid();

  const auto fit_params_begin = TimeFit::make_FitParameters_iterator(
    source_pdg_code, fTimeFits.begin() );
  const auto fit_params_end = TimeFit::make_FitParameters_iterator(
    source_pdg_code, fTimeFits.end() );

  const auto lum_begin = FitParameters::make_luminosity_iterator(
    fit_params_begin);
  const auto lum_end = FitParameters::make_luminosity_iterator(
    fit_params_end);

  std::discrete_distribution<size_t> time_dist(lum_begin, lum_end);

  if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM
    || fSamplingMode == TimeGenSamplingMode::UNIFORM_ENERGY)
  {
    time_bin_index = gen.sample_from_distribution(time_dist);
  }

  else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME) {
    int last_time_index = 0;
    if (fTimeFits.size() > 0) last_time_index = fTimeFits.size() - 1;
    std::uniform_int_distribution<size_t> uid(0, last_time_index);
    // for c2: time_bin_index has already been declared
    time_bin_index = gen.sample_from_distribution(uid);
  }

  else {
    throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
      << " encountered in evgen::MarleyTimeGen::produce()";
  }

  // Sample a time uniformly from within the selected time bin. Note that
  // the entries in fTimeFits use the time_ member to store the bin left
  // edges. The module creates fTimeFits in such a way that its last element
  // will always have luminosity_ == 0. (zero sampling weight), so we may
  // always add one to the sampled bin index without worrying about going
  // off the edge.
  double t_min = fTimeFits.at(time_bin_index).Time();
  double t_max = fTimeFits.at(time_bin_index + 1).Time();
  // sample a time on [ t_min, t_max )
  fTNu = gen.uniform_random_double(t_min, t_max, false);

  // Create a "beta-fit" neutrino source using the correct parameters for the
  // sampled time bin. This will be used to sample a neutrino energy unless
  // we're using the uniform time sampling mode. For uniform time sampling,
  // it will be used to determine the neutrino event weight.
  const auto& fit = fTimeFits.at(time_bin_index);
  std::unique_ptr<marley::NeutrinoSource> nu_source = source_from_time_fit(fit);

  if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM
    || fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME)
  {
    // Replace the generator's old source with the new one for the current
    // time bin
    gen.set_source(std::move(nu_source));

    // Generate a MARLEY event using the updated source
    mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());

    if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM) {
      // Unbiased sampling creates neutrino vertices with unit weight
      fWeight = ONE;
      sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection,
        sim::kUnbiased);
    }
    else {
      // fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME

      // Multiply by the likelihood ratio in order to correct for uniform
      // time sampling if we're using that biasing method
      double weight_bias = time_dist.probabilities().at(time_bin_index)
        / (t_max - t_min) * ( fTimeFits.back().Time()
        - fTimeFits.front().Time() );

      fWeight = weight_bias;
      sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection,
        sim::kUniformTime);
    }
  }

  else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_ENERGY)
  {
    double E_nu = std::numeric_limits<double>::lowest();
    mc_truth = make_uniform_energy_mctruth(fFitEmin, fFitEmax, E_nu,
      vertex_pos);

    // Get the value of the true dependent probability density (probability
    // of the sampled energy given the sampled time) to use as a biasing
    // correction in the neutrino vertex weight.

    // Load the generator with the neutrino source that represents the
    // true (i.e., unbiased) energy probability distribution. This will
    // create a normalized probability density that we can use to determine
    // the neutrino vertex weight.
    double nu_source_E_min = nu_source->get_Emin();
    double nu_source_E_max = nu_source->get_Emax();
    gen.set_source(std::move(nu_source));

    // NOTE: The marley::Generator object normalizes the E_pdf to unity
    // automatically, but just in case, we redo it here.
    double E_pdf_integ = integrate([&gen](double E_nu)
      -> double { return gen.E_pdf(E_nu); }, nu_source_E_min,
      nu_source_E_max);

    // Compute the likelihood ratio that we need to bias the neutrino vertex
    // weight
    double weight_bias = (gen.E_pdf(E_nu) / E_pdf_integ)
      * (fFitEmax - fFitEmin);

    fWeight = weight_bias;
    sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection,
      sim::kUniformEnergy);
  }

  else {
    throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
      << " encountered in evgen::MarleyTimeGen::produce()";
  }

}

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed )

inherited

Definition at line 26 of file EngineCreator.cc.

References art::EngineCreator::rng().

Referenced by evgen::CosmicsGen::CosmicsGen(), rndm::NuRandomService::createEngine(), cluster::fuzzyCluster::fuzzyCluster(), cluster::HoughLineFinder::HoughLineFinder(), art::MixFilter< T >::initEngine_(), larg4::LArG4::LArG4(), evgen::LightSource::LightSource(), evgen::NeutronOsc::NeutronOsc(), evgen::NucleonDecay::NucleonDecay(), opdet::OpMCDigi::OpMCDigi(), opdet::OptDetDigitizer::OptDetDigitizer(), phot::PhotonLibraryPropagation::PhotonLibraryPropagation(), detsim::SimDriftElectrons::SimDriftElectrons(), evgen::SingleGen::SingleGen(), evgen::SNNueAr40CCGen::SNNueAr40CCGen(), ToyOneShowerGen::ToyOneShowerGen(), and trkf::Track3DKalman::Track3DKalman().

{
  return rng()->createEngine(placeholder_schedule_id(), seed);
}

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed, std::string const &  kind_of_engine_to_make  )

inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

{
  return rng()->createEngine(
    placeholder_schedule_id(), seed, kind_of_engine_to_make);
}

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed, std::string const &  kind_of_engine_to_make, label_t const &  engine_label  )

inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

{
  return rng()->createEngine(
    placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
}

CurrentProcessingContext const * art::EDProducer::currentContext ( ) const

protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

  {
    return current_context_.get();
  }

EngineCreator::seed_t EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset, char const  key[] = "seed", seed_t const  implicit_seed = -1  )

inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

{
  auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
  return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
}

template<typename PROD , BranchType B>

ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {} ) const

inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

  {
    return ProducerBase::getProductID<PROD, B>(moduleDescription_,
                                               instanceName);
  }

template<typename PROD , BranchType B>

ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription, std::string const &  instanceName  ) const

inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

  {
    auto const& pd =
      get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
    return pd.productID();
  }

void evgen::MarleyTimeGen::make_final_timefit ( double  time )

protected

Helper function that makes a final dummy TimeFit object so that the final real time bin can have a right edge.

Definition at line 1269 of file MARLEYTimeGen_module.cc.

References fTimeFits.

Referenced by reconfigure().

{
  // The final time bin has zero luminosity, and therefore zero sampling
  // weight. We need it to be present so that the last nonzero weight bin
  // has a right edge.
  fTimeFits.emplace_back(time, 0., 0., 0., 0., 0., 0., 0., 0., 0.);
}

void evgen::MarleyTimeGen::make_nu_emission_histograms ( ) const

protected

Makes ROOT histograms showing the emitted neutrinos in each time bin when using a "fit"-format spectrum file.

Definition at line 1278 of file MARLEYTimeGen_module.cc.

References DEFINE_ART_MODULE, f, fTimeFits, evgen::MarleyTimeGen::TimeFit::GetFitParameters(), evgen::MarleyTimeGen::FitParameters::Luminosity(), art::TFileDirectory::make(), evgen::MarleyTimeGen::TimeFit::Time(), x_max, and x_min.

Referenced by reconfigure().

{
  // To make a histogram with variable size bins, ROOT needs the bin
  // low edges to be contiguous in memory. This is not true of the
  // stored times in fTimeFits, so we'll need to make a temporary copy
  // of them.
  std::vector<double> temp_times;
  for (const auto& fit : fTimeFits) temp_times.push_back(fit.Time());

  // If, for some reason, there are fewer than two time bins, just return
  // without making the histograms.
  // TODO: consider throwing an exception here
  if (temp_times.size() < 2) return;

  // Get the number of time bins
  int num_bins = temp_times.size() - 1;

  // Create some ROOT TH1D objects using the TFileService. These will store
  // the number of emitted neutrinos of each type in each time bin
  art::ServiceHandle<art::TFileService> tfs;
  TH1D* nue_hist = tfs->make<TH1D>("NueEmission",
    "Number of emitted #nu_{e}; time (s); number of neutrinos in time bin",
    num_bins, temp_times.data());
  TH1D* nuebar_hist = tfs->make<TH1D>("NuebarEmission",
    "Number of emitted #bar{#nu}_{e}; time (s); number of neutrinos in"
    " time bin", num_bins, temp_times.data());
  TH1D* nux_hist = tfs->make<TH1D>("NuxEmission",
    "Number of emitted #nu_{x}; time (s); number of neutrinos in time bin",
    num_bins, temp_times.data());

  // Load the histograms with the emitted neutrino counts
  for (size_t b = 1; b < temp_times.size(); ++b) {
    const TimeFit& current_fit = fTimeFits.at(b - 1);
    const TimeFit& next_fit = fTimeFits.at(b);
    double bin_deltaT = next_fit.Time() - current_fit.Time();

    const auto& nue_params = current_fit.GetFitParameters(
      marley_utils::ELECTRON_NEUTRINO);
    const auto& nuebar_params = current_fit.GetFitParameters(
      marley_utils::ELECTRON_ANTINEUTRINO);
    const auto& nux_params = current_fit.GetFitParameters(
      marley_utils::MUON_NEUTRINO);

    // Convert from bin luminosity to number of neutrinos by
    // multiplying by the bin time interval and dividing by the
    // mean neutrino energy
    double num_nue = 0.;
    double num_nuebar = 0.;
    double num_nux = 0.;
    if (nue_params.Emean() != 0.) num_nue = nue_params.Luminosity()
      * bin_deltaT / (nue_params.Emean() * MeV_to_erg);
    if (nuebar_params.Emean() != 0.) num_nuebar = nuebar_params.Luminosity()
      * bin_deltaT / (nuebar_params.Emean() * MeV_to_erg);
    if (nux_params.Emean() != 0.) num_nux = nux_params.Luminosity()
      * bin_deltaT / (nux_params.Emean() * MeV_to_erg);

    nue_hist->SetBinContent(b, num_nue);
    nuebar_hist->SetBinContent(b, num_nuebar);
    nux_hist->SetBinContent(b, num_nux);
  }

}

simb::MCTruth evgen::MarleyTimeGen::make_uniform_energy_mctruth ( double  E_min, double  E_max, double &  E_nu, const TLorentzVector &  vertex_pos  )

protected

Creates a simb::MCTruth object using a uniformly-sampled neutrino energy.

This function samples a reacting neutrino energy uniformly from the full range of energies allowed by the incident spectrum and the currently defined reactions.

Definition at line 1222 of file MARLEYTimeGen_module.cc.

References fEvent, and fMarleyHelper.

Referenced by create_truths_th2d(), and create_truths_time_fit().

{
  marley::Generator& gen = fMarleyHelper->get_generator();

  // Sample an energy uniformly over the entire allowed range
  double total_xs;
  int j = 0;
  do {
    // We have to check that the cross section is nonzero for the sampled
    // energy (otherwise we'll generate an unphysical event). However, if the
    // user has given us a histogram that is below threshold, the
    // program could get stuck here endlessly, sampling rejected energy
    // after rejected energy. Just in case, we cap the total number of tries
    // and quit if things don't work out.
    if (j >= MAX_UNIFORM_ENERGY_ITERATIONS) {
      throw cet::exception("MARLEYTimeGen") << "Exceeded the maximum of "
        << MAX_UNIFORM_ENERGY_ITERATIONS << " while attempting to sample"
        << " a neutrino energy uniformly.";
    }
    // Sample an energy uniformly on [ E_min, E_max )
    E_nu = gen.uniform_random_double(E_min, E_max, false);
    total_xs = 0.;
    // Check that at least one defined reaction has a non-zero total
    // cross section at the sampled energy. If this is not the case, try
    // again.
    for (const auto& react : gen.get_reactions()) {
      total_xs += react->total_xs(marley_utils::ELECTRON_NEUTRINO, E_nu);
    }

    ++j;
  } while (total_xs <= 0.);

  // Replace the existing neutrino source with a monoenergetic one at the
  // neutrino energy that was just sampled above
  std::unique_ptr<marley::NeutrinoSource> nu_source
    = std::make_unique<marley::MonoNeutrinoSource>(
    marley_utils::ELECTRON_NEUTRINO, E_nu);
  gen.set_source(std::move(nu_source));

  // Generate a MARLEY event using the new monoenergetic source
  auto mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());

  return mc_truth;
}

template<typename T , BranchType = InEvent>

ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it )

inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

{
  if (!moduleContext_)
    return ProductToken<T>::invalid();

  consumables_[BT].emplace_back(ConsumableType::Product,
                                TypeID{typeid(T)},
                                it.label(),
                                it.instance(),
                                it.process());
  return ProductToken<T>{it};
}

template<typename T , art::BranchType BT>

void art::Consumer::mayConsumeMany ( )

inherited

Definition at line 205 of file Consumer.h.

{
  if (!moduleContext_)
    return;

  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
}

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it )

inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

{
  if (!moduleContext_)
    return ViewToken<T>::invalid();

  consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                TypeID{typeid(T)},
                                it.label(),
                                it.instance(),
                                it.process());
  return ViewToken<T>{it};
}

bool art::ProducerBase::modifiesEvent ( ) const

inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

    {
      return true;
    }

cet::exempt_ptr< art::Consumer > art::Consumer::non_module_context ( )

staticinherited

Definition at line 76 of file Consumer.cc.

Referenced by art::RootOutput::beginSubRun(), art::OutputModule::doBeginRun(), art::OutputModule::doBeginSubRun(), art::OutputModule::doEndRun(), art::OutputModule::doEndSubRun(), art::ProducingService::doPostReadEvent(), art::ProducingService::doPostReadRun(), art::ProducingService::doPostReadSubRun(), art::OutputModule::doWriteEvent(), art::ProcessPackage< L >::postScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::Run >::End::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::End::postScheduleSignal(), art::ProcessPackage< L >::preScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::preScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::preScheduleSignal(), art::EventProcessor::readEvent(), art::EventProcessor::readRun(), art::EmptyEvent::readRun_(), art::EventProcessor::readSubRun(), and art::EmptyEvent::readSubRun_().

{
  static Consumer invalid{InvalidTag{}};
  return &invalid;
}

void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset )

protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

{
  if (!moduleContext_)
    return;

  pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
  for (auto& consumablesPerBranch : consumables_) {
    cet::sort_all(consumablesPerBranch);
  }
}

void evgen::MarleyTimeGen::produce ( art::Event &  e )

overridevirtual

Implements art::EDProducer.

Definition at line 717 of file MARLEYTimeGen_module.cc.

References create_truths_th2d(), create_truths_time_fit(), util::CreateAssn(), art::Event::event(), fEventNumber, fEventTree, FIT, fNeutrinosPerEvent, fRunNumber, fSpectrumFileFormat, fSubRunNumber, fTNu, fVertexSampler, n, art::Event::put(), RootTH2D, art::Event::run(), and art::Event::subRun().

{
  art::ServiceHandle<geo::Geometry> geo;

  // Get the run, subrun, and event numbers from the current art::Event
  fRunNumber = e.run();
  fSubRunNumber = e.subRun();
  fEventNumber = e.event();

  // Prepare associations and vectors of truth objects that will be produced
  // and loaded into the current art::Event
  std::unique_ptr< std::vector<simb::MCTruth> >
    truthcol(new std::vector<simb::MCTruth>);

  std::unique_ptr< std::vector<sim::SupernovaTruth> >
    sn_truthcol(new std::vector<sim::SupernovaTruth>);

  std::unique_ptr< art::Assns<simb::MCTruth, sim::SupernovaTruth> >
    truth_assns(new art::Assns<simb::MCTruth, sim::SupernovaTruth>);

  // Create temporary truth objects that we will use to load the event
  simb::MCTruth truth;
  sim::SupernovaTruth sn_truth;

  for (unsigned int n = 0; n < fNeutrinosPerEvent; ++n) {

    // Sample a primary vertex location for this event
    TLorentzVector vertex_pos = fVertexSampler->sample_vertex_pos(*geo);

    // Reset the neutrino's time-since-supernova to a bogus value (for now)
    fTNu = std::numeric_limits<double>::lowest();

    if (fSpectrumFileFormat == SpectrumFileFormat::RootTH2D) {
      create_truths_th2d(truth, sn_truth, vertex_pos);
    }
    else if (fSpectrumFileFormat == SpectrumFileFormat::FIT) {
      create_truths_time_fit(truth, sn_truth, vertex_pos);
    }
    else {
      throw cet::exception("MARLEYTimeGen") << "Invalid spectrum file"
        << " format encountered in evgen::MarleyTimeGen::produce()";
    }

    // Write the marley::Event object to the event tree
    fEventTree->Fill();

    // Add the truth objects to the appropriate vectors
    truthcol->push_back(truth);

    sn_truthcol->push_back(sn_truth);

    // Associate the last entries in each of the truth object vectors (the
    // truth objects that were just created for the current neutrino vertex)
    // with each other
    util::CreateAssn(*this, e, *truthcol, *sn_truthcol, *truth_assns,
      truthcol->size() - 1, truthcol->size()/*, sn_truthcol->size() - 1*/);
  }

  // Load the completed truth object vectors and associations into the event
  e.put(std::move(truthcol));

  e.put(std::move(sn_truthcol));

  e.put(std::move(truth_assns));
}

void evgen::MarleyTimeGen::reconfigure ( const Parameters &  p )

virtual

Definition at line 784 of file MARLEYTimeGen_module.cc.

References ALPHA, BETA, fFitEmax, fFitEmin, fFluxAveragedCrossSection, FIT, fMarleyHelper, fNeutrinosPerEvent, fPinchType, fSamplingMode, fSpectrumFileFormat, fSpectrumHist, fTimeFits, fVertexSampler, evgen::MarleyTimeGen::TimeFit::GetFitParameters(), HISTOGRAM, LOG_INFO, LOG_WARNING, evgen::MarleyTimeGen::FitParameters::Luminosity(), make_final_timefit(), make_nu_emission_histograms(), RootTH2D, s, source_from_time_fit(), UNIFORM_ENERGY, and UNIFORM_TIME.

Referenced by MarleyTimeGen().

{
  const auto& seed_service = art::ServiceHandle<rndm::NuRandomService>();
  const auto& geom_service = art::ServiceHandle<geo::Geometry>();

  // Create a new evgen::ActiveVolumeVertexSampler object based on the current
  // configuration
  fVertexSampler = std::make_unique<evgen::ActiveVolumeVertexSampler>(
    p().vertex_, *seed_service, *geom_service, "MARLEY_Vertex_Sampler");

  // Create a new marley::Generator object based on the current configuration
  fMarleyHelper = std::make_unique<MARLEYHelper>(p().marley_parameters_,
    *seed_service, "MARLEY");

  // Get the number of neutrino vertices per event from the FHiCL parameters
  fNeutrinosPerEvent = p().nu_per_event_();

  // Determine the current sampling mode from the FHiCL parameters
  const std::string& samp_mode_str = p().sampling_mode_();
  if (samp_mode_str == "histogram")
    fSamplingMode = TimeGenSamplingMode::HISTOGRAM;
  else if (samp_mode_str == "uniform time")
    fSamplingMode = TimeGenSamplingMode::UNIFORM_TIME;
  else if (samp_mode_str == "uniform energy")
    fSamplingMode = TimeGenSamplingMode::UNIFORM_ENERGY;
  else throw cet::exception("MARLEYTimeGen")
      << "Invalid sampling mode \"" << samp_mode_str << "\""
      << " specified for the MARLEYTimeGen module.";

  LOG_INFO("MARLEYTimeGen") << fNeutrinosPerEvent << " neutrino vertices"
    << " will be generated for each art::Event using the \"" << samp_mode_str
    << "\" sampling mode.";

  // Retrieve the time-dependent neutrino spectrum from the spectrum file.
  // Use different methods depending on the file's format.
  std::string spectrum_file_format = marley_utils::to_lowercase(
    p().spectrum_file_format_());

  if (spectrum_file_format == "th2d")
    fSpectrumFileFormat = SpectrumFileFormat::RootTH2D;
  else if (spectrum_file_format == "fit") {
    fSpectrumFileFormat = SpectrumFileFormat::FIT;

    std::string pinch_type;
    if ( !p().pinching_parameter_type_(pinch_type) ) {
      throw cet::exception("MARLEYTimeGen") << "Missing pinching parameter"
        << " type for a \"fit\" format spectrum file";
    }

    marley_utils::to_lowercase_inplace(pinch_type);
    if (pinch_type == "alpha") fPinchType = PinchParamType::ALPHA;
    else if (pinch_type == "beta") fPinchType = PinchParamType::BETA;
    else throw cet::exception("MARLEYTimeGen")
      << "Invalid pinching parameter type \"" << pinch_type
      << "\" specified for the MARLEYTimeGen module.";

    if ( !p().fit_Emin_(fFitEmin) ) throw cet::exception("MARLEYTimeGen")
      << "Missing minimum energy for a \"fit\" format spectrum"
      << " used by the MARLEYTimeGen module.";

    if ( !p().fit_Emax_(fFitEmax) ) throw cet::exception("MARLEYTimeGen")
      << "Missing maximum energy for a \"fit\" format spectrum"
      << " used by the MARLEYTimeGen module.";

    if (fFitEmax < fFitEmin) throw cet::exception("MARLEYTimeGen")
      << "Maximum energy is less than the minimum energy for"
      << " a \"fit\" format spectrum used by the MARLEYTimeGen module.";
  }
  else throw cet::exception("MARLEYTimeGen")
      << "Invalid spectrum file format \"" << p().spectrum_file_format_()
      << "\" specified for the MARLEYTimeGen module.";

  // Determine the full file name (including path) of the spectrum file
  std::string full_spectrum_file_name
    = fMarleyHelper->find_file(p().spectrum_file_(), "spectrum");

  marley::Generator& gen = fMarleyHelper->get_generator();

  if (fSpectrumFileFormat == SpectrumFileFormat::RootTH2D) {

    // Retrieve the time-dependent neutrino flux from a ROOT file
    std::unique_ptr<TFile> spectrum_file
      = std::make_unique<TFile>(full_spectrum_file_name.c_str(), "read");
    TH2D* temp_h2 = nullptr;
    std::string namecycle;
    if ( !p().namecycle_(namecycle) ) {
      throw cet::exception("MARLEYTimeGen") << "Missing namecycle for"
        << " a TH2D spectrum file";
    }

    spectrum_file->GetObject(namecycle.c_str(), temp_h2);
    fSpectrumHist.reset(temp_h2);

    // Disassociate the TH2D from its parent TFile. If we fail to do this,
    // then ROOT will auto-delete the TH2D when the TFile goes out of scope.
    fSpectrumHist->SetDirectory(nullptr);

    // Compute the flux-averaged total cross section using MARLEY. This will be
    // used to compute neutrino vertex weights for the sim::SupernovaTruth
    // objects.

    // Get a 1D projection of the energy spectrum (integrated over time)
    TH1D* energy_spect = fSpectrumHist->ProjectionY("energy_spect");

    // Create a new MARLEY neutrino source object using this projection
    // TODO: replace the hard-coded electron neutrino PDG code here (and in
    // several other places in this source file) when you're ready to use
    // MARLEY with multiple neutrino flavors
    std::unique_ptr<marley::NeutrinoSource> nu_source
      = marley_root::make_root_neutrino_source(marley_utils::ELECTRON_NEUTRINO,
      energy_spect);

    // Factor of hbar_c^2 converts from MeV^(-2) to fm^2
    fFluxAveragedCrossSection = marley_utils::hbar_c2
      * flux_averaged_total_xs(*nu_source, gen);

    // For speed, sample energies first whenever possible (and then sample from
    // an energy-dependent timing distribution). This avoids unnecessary calls
    // to MARLEY to change the energy spectrum.
    if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM ||
      fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME)
    {
      gen.set_source(std::move(nu_source));
    }

  } // spectrum_file_format == "th2d"

  else if (fSpectrumFileFormat == SpectrumFileFormat::FIT) {

    // Clear out the old parameterized spectrum, if one exists
    fTimeFits.clear();

    std::ifstream fit_file(full_spectrum_file_name);
    std::string line;

    bool found_end = false;

    // current line number
    int line_num = 0;
    // number of lines checked in last call to marley_utils::get_next_line()
    int lines_checked = 0;

    double old_time = std::numeric_limits<double>::lowest();

    while ( line = marley_utils::get_next_line(fit_file, rx_comment_or_empty,
      false, lines_checked), line_num += lines_checked, !line.empty() )
    {
      if (found_end) {
        LOG_WARNING("MARLEYTimeGen") << "Trailing content after last time"
          << " bin found on line " << line_num << " of the spectrum file "
          << full_spectrum_file_name;
      }

      double time;
      double Emean_nue, alpha_nue, luminosity_nue;
      double Emean_nuebar, alpha_nuebar, luminosity_nuebar;
      double Emean_nux, alpha_nux, luminosity_nux;
      std::istringstream iss(line);
      bool ok_first = static_cast<bool>( iss >> time );

      if (time <= old_time) throw cet::exception("MARLEYTimeGen")
        << "Time bin left edges given in the spectrum file must be"
        << " strictly increasing. Invalid time bin value found on line "
        << line_num << " of the spectrum file " << full_spectrum_file_name;
      else old_time = time;

      bool ok_rest = static_cast<bool>( iss >> Emean_nue >> alpha_nue
        >> luminosity_nue >> Emean_nuebar >> alpha_nuebar >> luminosity_nuebar
        >> Emean_nux >> alpha_nux >> luminosity_nux
      );

      if (ok_first) {
        // We haven't reached the final bin, so add another time bin
        // in the typical way.
        if (ok_rest) {
          fTimeFits.emplace_back(time, Emean_nue, alpha_nue, luminosity_nue,
            Emean_nuebar, alpha_nuebar, luminosity_nuebar, Emean_nux,
            alpha_nux, luminosity_nux);
        }
        else {
          make_final_timefit(time);
          found_end = true;
        }
      }
      else throw cet::exception("MARLEYTimeGen") << "Parse error on line "
          << line_num << " of the spectrum file " << full_spectrum_file_name;
    }

    if (!found_end) {

      size_t num_time_bins = fTimeFits.size();
      if (num_time_bins < 2) throw cet::exception("MARLEYTimeGen")
        << "Missing right edge for the final time bin in the spectrum file "
        << full_spectrum_file_name << ". Unable to guess a bin width for the "
        << " final bin.";

      // Guess that the width of the penultimate time bin and the width of
      // the final time bin are the same
      double delta_t_bin = fTimeFits.back().Time()
        - fTimeFits.at(num_time_bins - 2).Time();

      double last_bin_right_edge = fTimeFits.back().Time() + delta_t_bin;

      make_final_timefit(last_bin_right_edge);

      LOG_WARNING("MARLEYTimeGen") << "Missing right"
      << " edge for the final time bin in the spectrum file "
      << full_spectrum_file_name << ". Assuming a width of "
      << delta_t_bin << " s for the final bin.";
    }


    // Compute the flux-averaged total cross section for the fitted spectrum.
    // We will need this to compute neutrino vertex weights.
    std::vector< std::unique_ptr< marley::NeutrinoSource > > fit_sources;
    for (const auto& fit : fTimeFits) {
      fit_sources.emplace_back( source_from_time_fit(fit) );
    }

    // TODO: add handling for non-nu_e neutrino types when suitable data become
    // available in MARLEY
    auto temp_source = std::make_unique<marley::FunctionNeutrinoSource>(
      marley_utils::ELECTRON_NEUTRINO, fFitEmin, fFitEmax,
      [&fit_sources, this](double E_nu) -> double {
        double flux = 0.;
        for (size_t s = 0; s < fit_sources.size(); ++s) {

          const TimeFit& current_fit = this->fTimeFits.at(s);
          const auto& current_source = fit_sources.at(s);
          const FitParameters& params = current_fit.GetFitParameters(
            current_source->get_pid() );

          double lum = params.Luminosity();

          // Skip entries with zero luminosity, since they won't contribute
          // anything to the overall integral. Skip negative luminosity ones as
          // well, just in case.
          if (lum <= 0.) continue;

          flux += lum * current_source->pdf(E_nu);
        }
        return flux;
      }
    );

    double flux_integ = 0.;
    double tot_xs_integ = 0.;
    flux_averaged_total_xs(*temp_source, gen, flux_integ, tot_xs_integ);

    // Factor of hbar_c^2 converts from MeV^(-2) to fm^2
    fFluxAveragedCrossSection =  marley_utils::hbar_c2
      * tot_xs_integ / flux_integ;

    make_nu_emission_histograms();

  } // spectrum_file_format == "fit"

  else {
    throw cet::exception("MARLEYTimeGen") << "Unrecognized neutrino spectrum"
      << " file format \"" << p().spectrum_file_format_() << "\" encountered"
      << " in evgen::MarleyTimeGen::reconfigure()";
  }

  LOG_INFO("MARLEYTimeGen") << "The flux-averaged total cross section"
    << " predicted by MARLEY for the current supernova spectrum is "
    << fFluxAveragedCrossSection << " fm^2";

}

void art::Consumer::showMissingConsumes ( ) const

protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

{
  if (!moduleContext_)
    return;

  // If none of the branches have missing consumes statements, exit early.
  if (std::all_of(cbegin(missingConsumes_),
                  cend(missingConsumes_),
                  [](auto const& perBranch) { return perBranch.empty(); }))
    return;

  constexpr cet::HorizontalRule rule{60};
  mf::LogPrint log{"MTdiagnostics"};
  log << '\n'
      << rule('=') << '\n'
      << "The following consumes (or mayConsume) statements are missing from\n"
      << module_context(moduleDescription_) << '\n'
      << rule('-') << '\n';

  cet::for_all_with_index(
    missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
      for (auto const& pi : perBranch) {
        log << "  "
            << assemble_consumes_statement(static_cast<BranchType>(i), pi)
            << '\n';
      }
    });
  log << rule('=');
}

std::unique_ptr< marley::NeutrinoSource > evgen::MarleyTimeGen::source_from_time_fit ( const TimeFit &  fit )

protected

Create a MARLEY neutrino source object using a set of fit parameters for a particular time bin.

Definition at line 1195 of file MARLEYTimeGen_module.cc.

References evgen::MarleyTimeGen::FitParameters::Alpha(), ALPHA, beta, BETA, evgen::MarleyTimeGen::FitParameters::Emean(), fFitEmax, fFitEmin, fPinchType, and evgen::MarleyTimeGen::TimeFit::GetFitParameters().

Referenced by create_truths_time_fit(), and reconfigure().

{
  // Create a "beta-fit" neutrino source using the given fit parameters.

  // The two common fitting schemes (alpha and beta) differ in their
  // definitions by \beta = \alpha + 1.
  // TODO: add handling for non-nu_e neutrino types
  const FitParameters& nue_parameters
    = fit.GetFitParameters(marley_utils::ELECTRON_NEUTRINO);

  double beta = nue_parameters.Alpha();
  if (fPinchType == PinchParamType::ALPHA) beta += 1.;
  else if (fPinchType != PinchParamType::BETA) {
    throw cet::exception("MARLEYTimeGen") << "Unreognized pinching parameter"
      << " type encountered in evgen::MarleyTimeGen::source_from_time_fit()";
  }

  // Create the new source
  std::unique_ptr<marley::NeutrinoSource> nu_source
    = std::make_unique<marley::BetaFitNeutrinoSource>(
    marley_utils::ELECTRON_NEUTRINO, fFitEmin, fFitEmax,
    nue_parameters.Emean(), beta);

  return nu_source;
}

void art::Consumer::validateConsumedProduct ( BranchType const  bt, ProductInfo const &  pi  )

protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

{
  // Early exits if consumes tracking has been disabled or if the
  // consumed product is an allowed consumable.
  if (!moduleContext_)
    return;

  if (cet::binary_search_all(consumables_[bt], pi))
    return;

  if (requireConsumes_) {
    throw Exception(errors::ProductRegistrationFailure,
                    "Consumer: an error occurred during validation of a "
                    "retrieved product\n\n")
      << "The following consumes (or mayConsume) statement is missing from\n"
      << module_context(moduleDescription_) << ":\n\n"
      << "  " << assemble_consumes_statement(bt, pi) << "\n\n";
  }

  missingConsumes_[bt].insert(pi);
}

Member Data Documentation

std::unique_ptr<marley::Event> evgen::MarleyTimeGen::fEvent

protected

unique_ptr to the current event created by MARLEY

Definition at line 406 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), create_truths_time_fit(), make_uniform_energy_mctruth(), and MarleyTimeGen().

uint_fast32_t evgen::MarleyTimeGen::fEventNumber

protected

Event number from the art::Event being processed.

Definition at line 503 of file MARLEYTimeGen_module.cc.

Referenced by MarleyTimeGen(), and produce().

TTree* evgen::MarleyTimeGen::fEventTree

protected

The event TTree created by MARLEY.

This tree will be saved to the "hist" output file for validation purposes. The tree contains the same information as the generated simb::MCTruth objects, but in MARLEY's internal format

Definition at line 496 of file MARLEYTimeGen_module.cc.

Referenced by MarleyTimeGen(), and produce().

double evgen::MarleyTimeGen::fFitEmax

protected

Maximum neutrino energy to consider when using a "fit"-format spectrum file.

Definition at line 527 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_time_fit(), reconfigure(), and source_from_time_fit().

double evgen::MarleyTimeGen::fFitEmin

protected

Minimum neutrino energy to consider when using a "fit"-format spectrum file.

Definition at line 523 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_time_fit(), reconfigure(), and source_from_time_fit().

double evgen::MarleyTimeGen::fFluxAveragedCrossSection

protected

Flux-averaged total cross section (fm², average is taken over all energies and times for all defined reactions) used by MARLEY to generate neutrino vertices.

Definition at line 515 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), create_truths_time_fit(), and reconfigure().

std::unique_ptr<evgen::MARLEYHelper> evgen::MarleyTimeGen::fMarleyHelper

protected

Object that provides an interface to the MARLEY event generator.

Definition at line 399 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), create_truths_time_fit(), make_uniform_energy_mctruth(), and reconfigure().

unsigned int evgen::MarleyTimeGen::fNeutrinosPerEvent

protected

The number of MARLEY neutrino vertices to generate in each art::Event.

Definition at line 519 of file MARLEYTimeGen_module.cc.

Referenced by produce(), and reconfigure().

PinchParamType evgen::MarleyTimeGen::fPinchType

protected

The pinching parameter convention to use when interpreting the time-dependent fits.

Definition at line 470 of file MARLEYTimeGen_module.cc.

Referenced by reconfigure(), and source_from_time_fit().

uint_fast32_t evgen::MarleyTimeGen::fRunNumber

protected

Run number from the art::Event being processed.

Definition at line 499 of file MARLEYTimeGen_module.cc.

Referenced by MarleyTimeGen(), and produce().

TimeGenSamplingMode evgen::MarleyTimeGen::fSamplingMode

protected

Represents the sampling mode to use when selecting neutrino times and energies.

Definition at line 439 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), create_truths_time_fit(), and reconfigure().

SpectrumFileFormat evgen::MarleyTimeGen::fSpectrumFileFormat

protected

Format to assume for the neutrino spectrum input file.

Definition at line 490 of file MARLEYTimeGen_module.cc.

Referenced by produce(), and reconfigure().

std::unique_ptr<TH2D> evgen::MarleyTimeGen::fSpectrumHist

protected

ROOT TH2D that contains the time-dependent spectrum to use when sampling neutrino times and energies.

This member is only used when reading the spectrum from a ROOT file.

Definition at line 412 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), and reconfigure().

uint_fast32_t evgen::MarleyTimeGen::fSubRunNumber

protected

Subrun number from the art::Event being processed.

Definition at line 501 of file MARLEYTimeGen_module.cc.

Referenced by MarleyTimeGen(), and produce().

std::vector<TimeFit> evgen::MarleyTimeGen::fTimeFits

protected

Vector that contains the fit parameter information for each time bin when using a "fit"-format spectrum file.

This member is unused when the spectrum is read from a ROOT file.

Definition at line 418 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_time_fit(), make_final_timefit(), make_nu_emission_histograms(), and reconfigure().

double evgen::MarleyTimeGen::fTNu

protected

Time since the supernova core bounce for the current MARLEY neutrino vertex.

Definition at line 507 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), create_truths_time_fit(), MarleyTimeGen(), and produce().

std::unique_ptr<evgen::ActiveVolumeVertexSampler> evgen::MarleyTimeGen::fVertexSampler

protected

Algorithm that allows us to sample vertex locations within the active volume(s) of the detector.

Definition at line 403 of file MARLEYTimeGen_module.cc.

Referenced by produce(), and reconfigure().

double evgen::MarleyTimeGen::fWeight

protected

Statistical weight for the current MARLEY neutrino vertex.

Definition at line 510 of file MARLEYTimeGen_module.cc.

Referenced by create_truths_th2d(), create_truths_time_fit(), and MarleyTimeGen().

---

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

• /cvmfs/larsoft.opensciencegrid.org/products/larsim/v07_09_00/source/larsim/EventGenerator/MARLEY/MARLEYTimeGen_module.cc