MARLEYTimeGen_module.cc

Go to the documentation of this file.

// standard library includes
#include <fstream>
#include <limits>
#include <memory>
#include <regex>
#include <string>

// framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Persistency/Common/PtrMaker.h"
#include "art_root_io/TFileDirectory.h"
#include "art_root_io/TFileService.h"
#include "canvas/Persistency/Common/Assns.h"
#include "cetlib_except/exception.h"
#include "fhiclcpp/ParameterSet.h"
#include "fhiclcpp/types/OptionalAtom.h"
#include "fhiclcpp/types/Table.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// art extensions
#include "nurandom/RandomUtils/NuRandomService.h"

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "larcoreobj/SummaryData/RunData.h"
#include "lardataobj/Simulation/SupernovaTruth.h"
#include "larsim/EventGenerator/MARLEY/ActiveVolumeVertexSampler.h"
#include "larsim/EventGenerator/MARLEY/MARLEYHelper.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// ROOT includes
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TTree.h"

// MARLEY includes
#include "marley/Integrator.hh"
#include "marley/marley_root.hh"
#include "marley/marley_utils.hh"

// Anonymous namespace for definitions local to this source file
namespace {

  // The number of times to attempt to sample an energy uniformly
  // before giving up
  constexpr int MAX_UNIFORM_ENERGY_ITERATIONS = 1000;

  // Neutrino vertices generated using unbiased sampling are assigned
  // unit weight
  constexpr double ONE = 1.;

  // Number of sampling points to use for numerical integration
  // (via the Clenshaw-Curtis method)
  constexpr int NUM_INTEGRATION_SAMPLING_POINTS = 100;

  // Multiply an energy in MeV by this factor to convert to erg
  // (based on 2017 PDG "Physical Constants" article)
  constexpr double MeV_to_erg = 1.6021766208e-6; // erg / MeV

  // The PDG codes that represent one of the varieties of "nu_x"
  constexpr std::array<int, 4> NU_X_PDG_CODES{
    marley_utils::MUON_NEUTRINO,
    marley_utils::MUON_ANTINEUTRINO,
    marley_utils::TAU_NEUTRINO,
    marley_utils::TAU_ANTINEUTRINO,
  };

  // Helper function that tests whether a given PDG code represents
  // a "nu_x"
  bool is_nux(int pdg_code)
  {
    return std::find(std::begin(NU_X_PDG_CODES), std::end(NU_X_PDG_CODES), pdg_code) !=
           std::end(NU_X_PDG_CODES);
  }

  // Matches comment lines and empty lines in a "fit"-format spectrum file
  const std::regex rx_comment_or_empty = std::regex("\\s*(#.*)?");

  // Compute the flux-averaged total cross section (MeV^[-2]) for all
  // defined reactions and for a particular neutrino source
  double flux_averaged_total_xs(marley::NeutrinoSource& nu_source, marley::Generator& gen);

  // Overloaded version that allows access to the flux integral
  // (which will be loaded into source_integ) and the xs * flux integral
  // (which will be loaded into tot_xs_integ)
  double flux_averaged_total_xs(marley::NeutrinoSource& nu_source,
                                marley::Generator& gen,
                                double& source_integ,
                                double& tot_xs_integ);

  // Returns a numerical integral of the function f on the interval
  // [x_min, x_max]
  double integrate(const std::function<double(double)>& f, double x_min, double x_max);
}

namespace evgen {
  class MarleyTimeGen;
}

class evgen::MarleyTimeGen : public art::EDProducer {

public:
  using Name = fhicl::Name;
  using Comment = fhicl::Comment;

  struct KeysToIgnore {
    std::set<std::string> operator()() { return {"marley_parameters"}; }
  };

  struct Config {

    fhicl::Table<evgen::ActiveVolumeVertexSampler::Config> vertex_{
      Name("vertex"),
      Comment("Configuration for selecting the vertex location(s)")};

    fhicl::Atom<std::string> module_type_{
      Name("module_type"),
      Comment(""),
      "MARLEYTimeGen" // default value
    };

    fhicl::Atom<std::string> sampling_mode_{
      Name("sampling_mode"),
      Comment("Technique to use when sampling times and energies. Valid"
              " options are \"histogram\", \"uniform time\","
              " and \"uniform energy\""),
      std::string("histogram") // default value
    };

    fhicl::Atom<unsigned int> nu_per_event_{
      Name("nu_per_event"),
      Comment("The number of neutrino vertices to generate in"
              " each art::Event"),
      1u // default value
    };

    fhicl::Atom<std::string> spectrum_file_format_{
      Name("spectrum_file_format"),
      Comment("Format to assume for the neutrino spectrum file."
              " Valid options are \"th2d\" (a ROOT file containing a "
              " TH2D object) and \"fit\" (an ASCII text file containing"
              " fit parameters for each time bin). This parameter is not"
              " case-sensitive."),
      "th2d" // default value
    };

    fhicl::Atom<std::string> spectrum_file_{
      Name("spectrum_file"),
      Comment("Name of a file that contains a representation"
              " of the incident (not cross section weighted) neutrino spectrum"
              " as a function of time and energy.")};

    fhicl::OptionalAtom<std::string> pinching_parameter_type_{
      Name("pinching_parameter_type"),
      Comment("Type of pinching parameter to assume when parsing"
              " the time-dependent fit parameters for the incident neutrino"
              " spectrum. Valid options are \"alpha\" and \"beta\". This"
              " parameter is not case-sensitive."),
      [this]() -> bool {
        auto spectrum_file_format = marley_utils::to_lowercase(spectrum_file_format_());
        return (spectrum_file_format == "fit");
      }};

    fhicl::OptionalAtom<std::string> namecycle_{
      Name("namecycle"),
      Comment("Name of the TH2D object to use to represent the"
              " incident neutrino spectrum. This value should match the"
              " name of the TH2D as given in the ROOT file specified"
              " in the \"spectrum_file\" parameter. The TH2D should use "
              " time bins on the X axis (seconds) and energy bins on the "
              " Y axis (MeV)."),
      [this]() -> bool {
        auto spectrum_file_format = marley_utils::to_lowercase(spectrum_file_format_());
        return (spectrum_file_format == "th2d");
      }};

    fhicl::OptionalAtom<double> fit_Emin_{
      Name("fit_Emin"),
      Comment("Minimum allowed neutrino energy (MeV) for a \"fit\" format"
              " spectrum file"),
      [this]() -> bool {
        auto spectrum_file_format = marley_utils::to_lowercase(spectrum_file_format_());
        return (spectrum_file_format == "fit");
      }};

    fhicl::OptionalAtom<double> fit_Emax_{
      Name("fit_Emax"),
      Comment("Maximum allowed neutrino energy (MeV) for a \"fit\" format"
              " spectrum file"),
      [this]() -> bool {
        auto spectrum_file_format = marley_utils::to_lowercase(spectrum_file_format_());
        return (spectrum_file_format == "fit");
      }};

  }; // struct Config

  // Type to enable FHiCL parameter validation by art
  using Parameters = art::EDProducer::Table<Config, KeysToIgnore>;

  // @brief Configuration-checking constructor
  explicit MarleyTimeGen(const Parameters& p);

  virtual void produce(art::Event& e) override;
  virtual void beginRun(art::Run& run) override;

  virtual void reconfigure(const Parameters& p);

protected:
  class FitParameters {
  public:
    FitParameters(double Emean, double alpha, double luminosity)
      : fEmean(Emean), fAlpha(alpha), fLuminosity(luminosity)
    {}

    double Emean() const { return fEmean; }
    double Alpha() const { return fAlpha; }
    double Luminosity() const { return fLuminosity; }

    void set_Emean(double Emean) { fEmean = Emean; }
    void set_Alpha(double alpha) { fAlpha = alpha; }
    void set_Luminosity(double lum) { fLuminosity = lum; }

    template <typename It>
    static marley::IteratorToMember<It, double> make_luminosity_iterator(It it)
    {
      return marley::IteratorToMember<It, double>(it, &FitParameters::fLuminosity);
    }

  protected:
    double fEmean;      
    double fAlpha;      
    double fLuminosity; 
  };

  class TimeFit {
  public:
    TimeFit(double time,
            double Emean_nue,
            double alpha_nue,
            double lum_nue,
            double Emean_nuebar,
            double alpha_nuebar,
            double lum_nuebar,
            double Emean_nux,
            double alpha_nux,
            double lum_nux)
      : fTime(time)
      , fNueFitParams(Emean_nue, alpha_nue, lum_nue)
      , fNuebarFitParams(Emean_nuebar, alpha_nuebar, lum_nuebar)
      , fNuxFitParams(Emean_nux, alpha_nux, lum_nux)
    {}

    double Time() const { return fTime; }

  protected:
    static FitParameters TimeFit::*GetFitParametersMemberPointer(int pdg_code)
    {
      if (pdg_code == marley_utils::ELECTRON_NEUTRINO) { return &TimeFit::fNueFitParams; }
      else if (pdg_code == marley_utils::ELECTRON_ANTINEUTRINO) {
        return &TimeFit::fNuebarFitParams;
      }
      else if (is_nux(pdg_code)) {
        // The PDG code represents one of the varieties of nu_x
        return &TimeFit::fNuxFitParams;
      }
      else
        throw cet::exception("MARLEYTimeGen")
          << "Invalid neutrino"
          << " PDG code " << pdg_code << " encountered in MARLEYTimeGen"
          << "::TimeFit::GetFitParametersMemberPointer()";
      return nullptr;
    }

  public:
    const FitParameters& GetFitParameters(int pdg_code) const
    {
      return this->*GetFitParametersMemberPointer(pdg_code);
    }

    template <typename It>
    static marley::IteratorToMember<It, FitParameters> make_FitParameters_iterator(int pdg_code,
                                                                                   It iterator)
    {
      return marley::IteratorToMember<It, FitParameters>(iterator,
                                                         GetFitParametersMemberPointer(pdg_code));
    }

  protected:
    double fTime;

    FitParameters fNueFitParams;

    FitParameters fNuebarFitParams;

    FitParameters fNuxFitParams;
  };

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

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

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

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

  void make_final_timefit(double time);

  void make_nu_emission_histograms() const;

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

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

  std::unique_ptr<marley::Event> fEvent;

  std::unique_ptr<TH2D> fSpectrumHist;

  std::vector<TimeFit> fTimeFits;

  enum class TimeGenSamplingMode { HISTOGRAM, UNIFORM_TIME, UNIFORM_ENERGY };

  TimeGenSamplingMode fSamplingMode;

  enum class PinchParamType { ALPHA, BETA };

  PinchParamType fPinchType;

  // TODO: add information about the bin value units
  // TODO: add a "fit"-format description here
  enum class SpectrumFileFormat { RootTH2D, FIT };

  SpectrumFileFormat fSpectrumFileFormat;

  TTree* fEventTree;

  uint_fast32_t fRunNumber;
  uint_fast32_t fSubRunNumber;
  uint_fast32_t fEventNumber;

  double fTNu;

  double fWeight;

  double fFluxAveragedCrossSection; // fm^2

  unsigned int fNeutrinosPerEvent;

  double fFitEmin;

  double fFitEmax;
};

//------------------------------------------------------------------------------
evgen::MarleyTimeGen::MarleyTimeGen(const Parameters& p)
  : EDProducer{p.get_PSet()}
  , 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 const> 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>();
}

//------------------------------------------------------------------------------
void evgen::MarleyTimeGen::beginRun(art::Run& run)
{
  art::ServiceHandle<geo::Geometry const> geo;
  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
}

//------------------------------------------------------------------------------
void evgen::MarleyTimeGen::create_truths_th2d(simb::MCTruth& mc_truth,
                                              sim::SupernovaTruth& sn_truth,
                                              const TLorentzVector& vertex_pos)
{
  // 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::produce(art::Event& e)
{
  art::ServiceHandle<geo::Geometry const> 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
    truth_assns->addSingle(art::PtrMaker<simb::MCTruth>{e}(truthcol->size() - 1),
                           art::PtrMaker<sim::SupernovaTruth>{e}(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)
{
  const auto& seed_service = art::ServiceHandle<rndm::NuRandomService>();
  const auto& geom_service = art::ServiceHandle<geo::Geometry const>();

  // 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
  const fhicl::ParameterSet marley_pset =
    p.get_PSet().get<fhicl::ParameterSet>("marley_parameters");
  fMarleyHelper = std::make_unique<MARLEYHelper>(marley_pset, *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.";

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

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

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

//------------------------------------------------------------------------------
void evgen::MarleyTimeGen::create_truths_time_fit(simb::MCTruth& mc_truth,
                                                  sim::SupernovaTruth& sn_truth,
                                                  const TLorentzVector& vertex_pos)
{
  // 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()";
  }
}

//------------------------------------------------------------------------------
std::unique_ptr<marley::NeutrinoSource> evgen::MarleyTimeGen::source_from_time_fit(
  const TimeFit& fit)
{
  // 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;
}

//------------------------------------------------------------------------------
simb::MCTruth evgen::MarleyTimeGen::make_uniform_energy_mctruth(double E_min,
                                                                double E_max,
                                                                double& E_nu,
                                                                const TLorentzVector& vertex_pos)
{
  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;
}

//------------------------------------------------------------------------------
void evgen::MarleyTimeGen::make_final_timefit(double time)
{
  // 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
{
  // 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 const> 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);
  }
}

//------------------------------------------------------------------------------
// Anonymous namespace function definitions
namespace {

  double flux_averaged_total_xs(marley::NeutrinoSource& nu_source,
                                marley::Generator& gen,
                                double& source_integ,
                                double& tot_xs_integ)
  {
    // Get an integral of the source PDF (in case it isn't normalized to 1)
    source_integ = integrate([&nu_source](double E_nu) -> double { return nu_source.pdf(E_nu); },
                             nu_source.get_Emin(),
                             nu_source.get_Emax());

    tot_xs_integ = integrate(
      [&nu_source, &gen](double E_nu) -> double {
        double xs = 0.;
        for (const auto& react : gen.get_reactions()) {
          xs += react->total_xs(marley_utils::ELECTRON_NEUTRINO, E_nu);
        }
        return xs * nu_source.pdf(E_nu);
      },
      nu_source.get_Emin(),
      nu_source.get_Emax());

    return tot_xs_integ / source_integ;
  }

  double flux_averaged_total_xs(marley::NeutrinoSource& nu_source, marley::Generator& gen)
  {
    double dummy1, dummy2;
    return flux_averaged_total_xs(nu_source, gen, dummy1, dummy2);
  }

  double integrate(const std::function<double(double)>& f, double x_min, double x_max)
  {
    static marley::Integrator integrator(NUM_INTEGRATION_SAMPLING_POINTS);
    return integrator.num_integrate(f, x_min, x_max);
  }

}

DEFINE_ART_MODULE(evgen::MarleyTimeGen)