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 >

using	ModuleType = EDProducer

template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, 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)

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

void	fillProductDescriptions ()

void	registerProducts (ProductDescriptions &productsToRegister)

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

std::unique_ptr< Worker >	makeWorker (WorkerParams const &wp)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

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...

ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

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 (fm², 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 114 of file MARLEYTimeGen_module.cc.

Member Typedef Documentation

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

Definition at line 118 of file MARLEYTimeGen_module.cc.

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 17 of file EDProducer.h.

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

Definition at line 117 of file MARLEYTimeGen_module.cc.

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

Definition at line 215 of file MARLEYTimeGen_module.cc.

template<typename UserConfig , typename KeysToIgnore = void>

using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 26 of file Producer.h.

Member Enumeration Documentation

enum evgen::MarleyTimeGen::PinchParamType

strongprotected

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

Enumerator
ALPHA
BETA

Definition at line 459 of file MARLEYTimeGen_module.cc.

459 { ALPHA, BETA };

enum evgen::MarleyTimeGen::SpectrumFileFormat

strongprotected

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

Enumerator
RootTH2D
FIT

Definition at line 480 of file MARLEYTimeGen_module.cc.

480 { 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 428 of file MARLEYTimeGen_module.cc.

428 { HISTOGRAM, UNIFORM_TIME, UNIFORM_ENERGY };

Constructor & Destructor Documentation

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

explicit

Definition at line 524 of file MARLEYTimeGen_module.cc.

References fEvent, fEventNumber, fEventTree, fFluxAveragedCrossSection, fRunNumber, fSubRunNumber, fTNu, fWeight, art::ProducerTable< UserConfig, ImplicitConfig, UserKeysToIgnore >::get_PSet(), and reconfigure().

   : 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

Reimplemented from art::EDProducer.

Definition at line 558 of file MARLEYTimeGen_module.cc.

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

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

template<typename T , BranchType BT>

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

protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

   {
     return collector_.consumes<T, BT>(tag);
   }

ConsumesCollector & art::ModuleBase::consumesCollector ( )

protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

   {
     return collector_;
   }

template<typename T , BranchType BT>

void art::ModuleBase::consumesMany ( )

protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

   {
     collector_.consumesMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

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

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::consumesView ( InputTag const & tag )

inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

   {
     return collector_.consumesView<T, BT>(tag);
   }

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.

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

Referenced by produce().

 {
   // Get a reference to the generator object created by MARLEY (we'll need
   // to do a few fancy things with it other than just creating events)
   marley::Generator& gen = fMarleyHelper->get_generator();
 
   if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM) {
     // Generate a MARLEY event using the time-integrated spectrum
     // (the generator was already configured to use it by reconfigure())
     mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());
 
     // Find the time distribution corresponding to the selected energy bin
     double E_nu = fEvent->projectile().total_energy();
     int E_bin_index = fSpectrumHist->GetYaxis()->FindBin(E_nu);
     TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist", E_bin_index, E_bin_index);
     double* time_bin_weights = t_hist->GetArray();
 
     // Sample a time bin from the distribution
     std::discrete_distribution<int> time_dist;
     std::discrete_distribution<int>::param_type time_params(
       time_bin_weights, &(time_bin_weights[t_hist->GetNbinsX() + 1]));
     int time_bin_index = gen.sample_from_distribution(time_dist, time_params);
 
     // Sample a time uniformly from within the selected time bin
     double t_min = t_hist->GetBinLowEdge(time_bin_index);
     double t_max = t_min + t_hist->GetBinWidth(time_bin_index);
     // sample a time on [ t_min, t_max )
     fTNu = gen.uniform_random_double(t_min, t_max, false);
     // Unbiased sampling was used, so assign this neutrino vertex a
     // unit statistical weight
     fWeight = ONE;
 
     sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection, sim::kUnbiased);
   }
 
   else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME) {
     // Generate a MARLEY event using the time-integrated spectrum
     // (the generator was already configured to use it by reconfigure())
     mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());
 
     // Sample a time uniformly
     TAxis* time_axis = fSpectrumHist->GetXaxis();
     // underflow bin has index zero
     double t_min = time_axis->GetBinLowEdge(1);
     double t_max = time_axis->GetBinLowEdge(time_axis->GetNbins() + 1);
     // sample a time on [ t_min, t_max )
     fTNu = gen.uniform_random_double(t_min, t_max, false);
 
     // Get the value of the true dependent probability density (probability
     // of the sampled time given the sampled energy) to use as a biasing
     // correction in the neutrino vertex weight.
     double E_nu = fEvent->projectile().total_energy();
     int E_bin_index = fSpectrumHist->GetYaxis()->FindBin(E_nu);
     TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist", E_bin_index, E_bin_index);
     int t_bin_index = t_hist->FindBin(fTNu);
     double weight_bias = t_hist->GetBinContent(t_bin_index) * (t_max - t_min) /
                          (t_hist->Integral() * t_hist->GetBinWidth(t_bin_index));
 
     fWeight = weight_bias;
 
     sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection, sim::kUniformTime);
   }
 
   else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_ENERGY) {
     // Select a time bin using the energy-integrated spectrum
     TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist");
     double* time_bin_weights = t_hist->GetArray();
 
     // Sample a time bin from the distribution
     std::discrete_distribution<int> time_dist;
     std::discrete_distribution<int>::param_type time_params(
       time_bin_weights, &(time_bin_weights[t_hist->GetNbinsX() + 1]));
     int time_bin_index = gen.sample_from_distribution(time_dist, time_params);
 
     // Sample a time uniformly from within the selected time bin
     double t_min = t_hist->GetBinLowEdge(time_bin_index);
     double t_max = t_min + t_hist->GetBinWidth(time_bin_index);
     // sample a time on [ t_min, t_max )
     fTNu = gen.uniform_random_double(t_min, t_max, false);
 
     // Sample an energy uniformly over the entire allowed range
     // underflow bin has index zero
     TAxis* energy_axis = fSpectrumHist->GetYaxis();
     double E_min = energy_axis->GetBinLowEdge(1);
     double E_max = energy_axis->GetBinLowEdge(energy_axis->GetNbins() + 1);
 
     double E_nu = std::numeric_limits<double>::lowest();
 
     // Generate a MARLEY event using a uniformly sampled energy
     mc_truth = make_uniform_energy_mctruth(E_min, E_max, E_nu, vertex_pos);
 
     // Get the value of the true dependent probability density (probability
     // of the sampled energy given the sampled time) to use as a biasing
     // correction in the neutrino vertex weight.
     //
     // Get a 1D projection of the energy spectrum for the sampled time bin
     TH1D* energy_spect = fSpectrumHist->ProjectionY("energy_spect", time_bin_index, time_bin_index);
 
     // Create a new MARLEY neutrino source object using this projection (this
     // will create a normalized probability density that we can use) and load
     // it into the generator.
     auto nu_source =
       marley_root::make_root_neutrino_source(marley_utils::ELECTRON_NEUTRINO, energy_spect);
     double new_source_E_min = nu_source->get_Emin();
     double new_source_E_max = nu_source->get_Emax();
     gen.set_source(std::move(nu_source));
     // NOTE: The marley::Generator object normalizes the E_pdf to unity
     // automatically, but just in case, we redo it here.
     double E_pdf_integ = integrate([&gen](double E_nu) -> double { return gen.E_pdf(E_nu); },
                                    new_source_E_min,
                                    new_source_E_max);
 
     // Compute the likelihood ratio that we need to bias the neutrino vertex
     // weight
     double weight_bias = (gen.E_pdf(E_nu) / E_pdf_integ) * (E_max - E_min);
 
     fWeight = weight_bias;
 
     sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection, sim::kUniformEnergy);
   }
 
   else {
     throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
                                           << " encountered in evgen::MarleyTimeGen::produce()";
   }
 }

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

protected

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

Definition at line 1042 of file MARLEYTimeGen_module.cc.

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

Referenced by produce().

 {
   // Get a reference to the generator object created by MARLEY (we'll need
   // to do a few fancy things with it other than just creating events)
   marley::Generator& gen = fMarleyHelper->get_generator();
 
   // Initialize the time bin index to something absurdly large. This will help
   // us detect strange bugs that arise when it is sampled incorrectly.
   size_t time_bin_index = std::numeric_limits<size_t>::max();
 
   // Create an object to represent the discrete time distribution given in
   // the spectrum file. Use the luminosity for each bin as its sampling weight.
   // This distribution will actually be used to sample a neutrino arrival time
   // unless we're using the uniform time sampling mode, in which case it will
   // be used to help calculate the neutrino vertex weight.
   int source_pdg_code = gen.get_source().get_pid();
 
   const auto fit_params_begin =
     TimeFit::make_FitParameters_iterator(source_pdg_code, fTimeFits.begin());
   const auto fit_params_end =
     TimeFit::make_FitParameters_iterator(source_pdg_code, fTimeFits.end());
 
   const auto lum_begin = FitParameters::make_luminosity_iterator(fit_params_begin);
   const auto lum_end = FitParameters::make_luminosity_iterator(fit_params_end);
 
   std::discrete_distribution<size_t> time_dist(lum_begin, lum_end);
 
   if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM ||
       fSamplingMode == TimeGenSamplingMode::UNIFORM_ENERGY) {
     time_bin_index = gen.sample_from_distribution(time_dist);
   }
 
   else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME) {
     int last_time_index = 0;
     if (fTimeFits.size() > 0) last_time_index = fTimeFits.size() - 1;
     std::uniform_int_distribution<size_t> uid(0, last_time_index);
     // for c2: time_bin_index has already been declared
     time_bin_index = gen.sample_from_distribution(uid);
   }
 
   else {
     throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
                                           << " encountered in evgen::MarleyTimeGen::produce()";
   }
 
   // Sample a time uniformly from within the selected time bin. Note that
   // the entries in fTimeFits use the time_ member to store the bin left
   // edges. The module creates fTimeFits in such a way that its last element
   // will always have luminosity_ == 0. (zero sampling weight), so we may
   // always add one to the sampled bin index without worrying about going
   // off the edge.
   double t_min = fTimeFits.at(time_bin_index).Time();
   double t_max = fTimeFits.at(time_bin_index + 1).Time();
   // sample a time on [ t_min, t_max )
   fTNu = gen.uniform_random_double(t_min, t_max, false);
 
   // Create a "beta-fit" neutrino source using the correct parameters for the
   // sampled time bin. This will be used to sample a neutrino energy unless
   // we're using the uniform time sampling mode. For uniform time sampling,
   // it will be used to determine the neutrino event weight.
   const auto& fit = fTimeFits.at(time_bin_index);
   std::unique_ptr<marley::NeutrinoSource> nu_source = source_from_time_fit(fit);
 
   if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM ||
       fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME) {
     // Replace the generator's old source with the new one for the current
     // time bin
     gen.set_source(std::move(nu_source));
 
     // Generate a MARLEY event using the updated source
     mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());
 
     if (fSamplingMode == TimeGenSamplingMode::HISTOGRAM) {
       // Unbiased sampling creates neutrino vertices with unit weight
       fWeight = ONE;
       sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection, sim::kUnbiased);
     }
     else {
       // fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME
 
       // Multiply by the likelihood ratio in order to correct for uniform
       // time sampling if we're using that biasing method
       double weight_bias = time_dist.probabilities().at(time_bin_index) / (t_max - t_min) *
                            (fTimeFits.back().Time() - fTimeFits.front().Time());
 
       fWeight = weight_bias;
       sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection, sim::kUniformTime);
     }
   }
 
   else if (fSamplingMode == TimeGenSamplingMode::UNIFORM_ENERGY) {
     double E_nu = std::numeric_limits<double>::lowest();
     mc_truth = make_uniform_energy_mctruth(fFitEmin, fFitEmax, E_nu, vertex_pos);
 
     // Get the value of the true dependent probability density (probability
     // of the sampled energy given the sampled time) to use as a biasing
     // correction in the neutrino vertex weight.
 
     // Load the generator with the neutrino source that represents the
     // true (i.e., unbiased) energy probability distribution. This will
     // create a normalized probability density that we can use to determine
     // the neutrino vertex weight.
     double nu_source_E_min = nu_source->get_Emin();
     double nu_source_E_max = nu_source->get_Emax();
     gen.set_source(std::move(nu_source));
 
     // NOTE: The marley::Generator object normalizes the E_pdf to unity
     // automatically, but just in case, we redo it here.
     double E_pdf_integ = integrate(
       [&gen](double E_nu) -> double { return gen.E_pdf(E_nu); }, nu_source_E_min, nu_source_E_max);
 
     // Compute the likelihood ratio that we need to bias the neutrino vertex
     // weight
     double weight_bias = (gen.E_pdf(E_nu) / E_pdf_integ) * (fFitEmax - fFitEmin);
 
     fWeight = weight_bias;
     sn_truth = sim::SupernovaTruth(fTNu, fWeight, fFluxAveragedCrossSection, sim::kUniformEnergy);
   }
 
   else {
     throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
                                           << " encountered in evgen::MarleyTimeGen::produce()";
   }
 }

void art::detail::Producer::doBeginJob ( SharedResources const & resources )

inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

   {
     setupQueues(resources);
     ProcessingFrame const frame{ScheduleID{}};
     beginJobWithFrame(frame);
   }

bool art::detail::Producer::doBeginRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

   {
     auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

bool art::detail::Producer::doBeginSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

   {
     auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

void art::detail::Producer::doEndJob ( )

inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     endJobWithFrame(frame);
   }

bool art::detail::Producer::doEndRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

   {
     auto r = rp.makeRun(mc, rp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

bool art::detail::Producer::doEndSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

   {
     auto sr = srp.makeSubRun(mc, srp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

bool art::detail::Producer::doEvent	(	EventPrincipal &	ep,
		ModuleContext const &	mc,
		std::atomic< std::size_t > &	counts_run,
		std::atomic< std::size_t > &	counts_passed,
		std::atomic< std::size_t > &	counts_failed
	)

inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

   {
     auto e = ep.makeEvent(mc);
     ++counts_run;
     ProcessingFrame const frame{mc.scheduleID()};
     produceWithFrame(e, frame);
     e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
     ++counts_passed;
     return true;
   }

void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const & fb )

inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const & fb )

inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseOutputFilesWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const & fb )

inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const & fb )

inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenOutputFilesWithFrame(fb, frame);
   }

void art::Modifier::fillProductDescriptions ( )

inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

   {
     ProductRegistryHelper::fillDescriptions(moduleDescription());
   }

std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const

inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

   {
     return collector_.getConsumables();
   }

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

References fTimeFits.

Referenced by reconfigure().

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

void evgen::MarleyTimeGen::make_nu_emission_histograms ( ) const

protected

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

Definition at line 1256 of file MARLEYTimeGen_module.cc.

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

Referenced by reconfigure().

 {
   // To make a histogram with variable size bins, ROOT needs the bin
   // low edges to be contiguous in memory. This is not true of the
   // stored times in fTimeFits, so we'll need to make a temporary copy
   // of them.
   std::vector<double> temp_times;
   for (const auto& fit : fTimeFits)
     temp_times.push_back(fit.Time());
 
   // If, for some reason, there are fewer than two time bins, just return
   // without making the histograms.
   // TODO: consider throwing an exception here
   if (temp_times.size() < 2) return;
 
   // Get the number of time bins
   int num_bins = temp_times.size() - 1;
 
   // Create some ROOT TH1D objects using the TFileService. These will store
   // the number of emitted neutrinos of each type in each time bin
   art::ServiceHandle<art::TFileService 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 1198 of file MARLEYTimeGen_module.cc.

References fEvent, and fMarleyHelper.

Referenced by create_truths_th2d(), and create_truths_time_fit().

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

std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const & wp )

inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

   {
     return doMakeWorker(wp);
   }

template<typename T , BranchType BT>

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

protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

   {
     return collector_.mayConsume<T, BT>(tag);
   }

template<typename T , BranchType BT>

void art::ModuleBase::mayConsumeMany ( )

protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

   {
     collector_.mayConsumeMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

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

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const & tag )

inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

   {
     return collector_.mayConsumeView<T, BT>(tag);
   }

ModuleDescription const & art::ModuleBase::moduleDescription ( ) const

inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

   {
     if (md_.has_value()) {
       return *md_;
     }
 
     throw Exception{errors::LogicError,
                     "There was an error while calling moduleDescription().\n"}
       << "The moduleDescription() base-class member function cannot be called\n"
          "during module construction.  To determine which module is "
          "responsible\n"
          "for calling it, find the '<module type>:<module "
          "label>@Construction'\n"
          "tag in the message prefix above.  Please contact artists@fnal.gov\n"
          "for guidance.\n";
   }

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

overridevirtual

Implements art::EDProducer.

Definition at line 695 of file MARLEYTimeGen_module.cc.

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

 {
   art::ServiceHandle<geo::Geometry 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 762 of file MARLEYTimeGen_module.cc.

References ALPHA, BETA, fFitEmax, fFitEmin, fFluxAveragedCrossSection, FIT, fMarleyHelper, fNeutrinosPerEvent, fPinchType, fSamplingMode, fSpectrumFileFormat, fSpectrumHist, fTimeFits, fVertexSampler, art::ProducerTable< UserConfig, ImplicitConfig, UserKeysToIgnore >::get_PSet(), evgen::MarleyTimeGen::TimeFit::GetFitParameters(), HISTOGRAM, evgen::MarleyTimeGen::FitParameters::Luminosity(), make_final_timefit(), make_nu_emission_histograms(), MF_LOG_INFO, MF_LOG_WARNING, RootTH2D, source_from_time_fit(), UNIFORM_ENERGY, and UNIFORM_TIME.

Referenced by MarleyTimeGen().

 {
   const auto& seed_service = art::ServiceHandle<rndm::NuRandomService>();
   const auto& geom_service = art::ServiceHandle<geo::Geometry 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 art::Modifier::registerProducts ( ProductDescriptions & productsToRegister )

inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

   {
     ProductRegistryHelper::registerProducts(productsToRegister,
                                             moduleDescription());
   }

void art::ModuleBase::setModuleDescription ( ModuleDescription const & md )

inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

   {
     md_ = md;
   }

void art::ModuleBase::sortConsumables ( std::string const & current_process_name )

inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

   {
     // Now that we know we have seen all the consumes declarations,
     // sort the results for fast lookup later.
     collector_.sortConsumables(current_process_name);
   }

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

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

Referenced by create_truths_time_fit(), and reconfigure().

 {
   // Create a "beta-fit" neutrino source using the given fit parameters.
 
   // The two common fitting schemes (alpha and beta) differ in their
   // definitions by \beta = \alpha + 1.
   // TODO: add handling for non-nu_e neutrino types
   const FitParameters& nue_parameters = fit.GetFitParameters(marley_utils::ELECTRON_NEUTRINO);
 
   double beta = nue_parameters.Alpha();
   if (fPinchType == PinchParamType::ALPHA)
     beta += 1.;
   else if (fPinchType != PinchParamType::BETA) {
     throw cet::exception("MARLEYTimeGen")
       << "Unreognized pinching parameter"
       << " type encountered in evgen::MarleyTimeGen::source_from_time_fit()";
   }
 
   // Create the new source
   std::unique_ptr<marley::NeutrinoSource> nu_source =
     std::make_unique<marley::BetaFitNeutrinoSource>(
       marley_utils::ELECTRON_NEUTRINO, fFitEmin, fFitEmax, nue_parameters.Emean(), beta);
 
   return nu_source;
 }

Member Data Documentation

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

protected

unique_ptr to the current event created by MARLEY

Definition at line 399 of file MARLEYTimeGen_module.cc.

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

uint_fast32_t evgen::MarleyTimeGen::fEventNumber

protected

Event number from the art::Event being processed.

Definition at line 496 of file MARLEYTimeGen_module.cc.

Referenced by MarleyTimeGen(), and produce().

TTree* evgen::MarleyTimeGen::fEventTree

protected

The event TTree created by MARLEY.

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

Definition at line 489 of file MARLEYTimeGen_module.cc.

Referenced by MarleyTimeGen(), and produce().

double evgen::MarleyTimeGen::fFitEmax

protected

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

Definition at line 520 of file MARLEYTimeGen_module.cc.

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

double evgen::MarleyTimeGen::fFitEmin

protected

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

Definition at line 516 of file MARLEYTimeGen_module.cc.

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

double evgen::MarleyTimeGen::fFluxAveragedCrossSection

protected

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

Definition at line 508 of file MARLEYTimeGen_module.cc.

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