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...
struct	KeysToIgnore
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, KeysToIgnore >

Public Member Functions
	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 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...

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 112 of file MARLEYTimeGen_module.cc.

Member Typedef Documentation

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

Definition at line 117 of file MARLEYTimeGen_module.cc.

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

Definition at line 116 of file MARLEYTimeGen_module.cc.

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

Definition at line 227 of file MARLEYTimeGen_module.cc.

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 464 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 485 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 433 of file MARLEYTimeGen_module.cc.

{ HISTOGRAM, UNIFORM_TIME, UNIFORM_ENERGY };

Constructor & Destructor Documentation

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

explicit

Definition at line 529 of file MARLEYTimeGen_module.cc.

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

Member Function Documentation

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

overridevirtual

Definition at line 559 of file MARLEYTimeGen_module.cc.

                                             {
  art::ServiceHandle<geo::Geometry const> geo;
  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
}

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 565 of file MARLEYTimeGen_module.cc.

{
  // 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 1045 of file MARLEYTimeGen_module.cc.

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

}

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 1260 of file MARLEYTimeGen_module.cc.

{
  // 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 1269 of file MARLEYTimeGen_module.cc.

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

}

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 1213 of file MARLEYTimeGen_module.cc.

{
  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::produce ( art::Event &  e )

overridevirtual

Definition at line 706 of file MARLEYTimeGen_module.cc.

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

virtual

Definition at line 773 of file MARLEYTimeGen_module.cc.

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

}

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 1186 of file MARLEYTimeGen_module.cc.

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

Member Data Documentation

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

protected

unique_ptr to the current event created by MARLEY

Definition at line 404 of file MARLEYTimeGen_module.cc.

uint_fast32_t evgen::MarleyTimeGen::fEventNumber

protected

Event number from the art::Event being processed.

Definition at line 501 of file MARLEYTimeGen_module.cc.

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 494 of file MARLEYTimeGen_module.cc.

double evgen::MarleyTimeGen::fFitEmax

protected

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

Definition at line 525 of file MARLEYTimeGen_module.cc.

double evgen::MarleyTimeGen::fFitEmin

protected

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

Definition at line 521 of file MARLEYTimeGen_module.cc.

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 513 of file MARLEYTimeGen_module.cc.

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

protected

Object that provides an interface to the MARLEY event generator.

Definition at line 397 of file MARLEYTimeGen_module.cc.

unsigned int evgen::MarleyTimeGen::fNeutrinosPerEvent

protected

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

Definition at line 517 of file MARLEYTimeGen_module.cc.

PinchParamType evgen::MarleyTimeGen::fPinchType

protected

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

Definition at line 468 of file MARLEYTimeGen_module.cc.

uint_fast32_t evgen::MarleyTimeGen::fRunNumber

protected

Run number from the art::Event being processed.

Definition at line 497 of file MARLEYTimeGen_module.cc.

TimeGenSamplingMode evgen::MarleyTimeGen::fSamplingMode

protected

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

Definition at line 437 of file MARLEYTimeGen_module.cc.

SpectrumFileFormat evgen::MarleyTimeGen::fSpectrumFileFormat

protected

Format to assume for the neutrino spectrum input file.

Definition at line 488 of file MARLEYTimeGen_module.cc.

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 410 of file MARLEYTimeGen_module.cc.

uint_fast32_t evgen::MarleyTimeGen::fSubRunNumber

protected

Subrun number from the art::Event being processed.

Definition at line 499 of file MARLEYTimeGen_module.cc.

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 416 of file MARLEYTimeGen_module.cc.

double evgen::MarleyTimeGen::fTNu

protected

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

Definition at line 505 of file MARLEYTimeGen_module.cc.

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 401 of file MARLEYTimeGen_module.cc.

double evgen::MarleyTimeGen::fWeight

protected

Statistical weight for the current MARLEY neutrino vertex.

Definition at line 508 of file MARLEYTimeGen_module.cc.

---

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

• MARLEYTimeGen_module.cc