779 fVertexSampler = std::make_unique<evgen::ActiveVolumeVertexSampler>(
780 p().vertex_, *seed_service, *geom_service,
"MARLEY_Vertex_Sampler");
786 *seed_service,
"MARLEY" );
792 const std::string& samp_mode_str =
p().sampling_mode_();
793 if (samp_mode_str ==
"histogram")
795 else if (samp_mode_str ==
"uniform time")
797 else if (samp_mode_str ==
"uniform energy")
800 <<
"Invalid sampling mode \"" << samp_mode_str <<
"\"" 801 <<
" specified for the MARLEYTimeGen module.";
803 MF_LOG_INFO(
"MARLEYTimeGen") << fNeutrinosPerEvent <<
" neutrino vertices" 804 <<
" will be generated for each art::Event using the \"" << samp_mode_str
805 <<
"\" sampling mode.";
809 std::string spectrum_file_format = marley_utils::to_lowercase(
810 p().spectrum_file_format_());
812 if (spectrum_file_format ==
"th2d")
814 else if (spectrum_file_format ==
"fit") {
818 if ( !
p().pinching_parameter_type_(pinch_type) ) {
819 throw cet::exception(
"MARLEYTimeGen") <<
"Missing pinching parameter" 820 <<
" type for a \"fit\" format spectrum file";
823 marley_utils::to_lowercase_inplace(pinch_type);
827 <<
"Invalid pinching parameter type \"" << pinch_type
828 <<
"\" specified for the MARLEYTimeGen module.";
831 <<
"Missing minimum energy for a \"fit\" format spectrum" 832 <<
" used by the MARLEYTimeGen module.";
835 <<
"Missing maximum energy for a \"fit\" format spectrum" 836 <<
" used by the MARLEYTimeGen module.";
839 <<
"Maximum energy is less than the minimum energy for" 840 <<
" a \"fit\" format spectrum used by the MARLEYTimeGen module.";
843 <<
"Invalid spectrum file format \"" <<
p().spectrum_file_format_()
844 <<
"\" specified for the MARLEYTimeGen module.";
848 = fMarleyHelper->find_file(
p().spectrum_file_(),
"spectrum");
850 marley::Generator&
gen = fMarleyHelper->get_generator();
855 std::unique_ptr<TFile> spectrum_file
856 = std::make_unique<TFile>(full_spectrum_file_name.c_str(),
"read");
857 TH2D* temp_h2 =
nullptr;
859 if ( !
p().namecycle_(namecycle) ) {
861 <<
" a TH2D spectrum file";
864 spectrum_file->GetObject(namecycle.c_str(), temp_h2);
876 TH1D* energy_spect =
fSpectrumHist->ProjectionY(
"energy_spect");
882 std::unique_ptr<marley::NeutrinoSource> nu_source
883 = marley_root::make_root_neutrino_source(marley_utils::ELECTRON_NEUTRINO,
888 * flux_averaged_total_xs(*nu_source, gen);
906 std::ifstream fit_file(full_spectrum_file_name);
909 bool found_end =
false;
914 int lines_checked = 0;
916 double old_time = std::numeric_limits<double>::lowest();
918 while ( line = marley_utils::get_next_line(fit_file, rx_comment_or_empty,
919 false, lines_checked), line_num += lines_checked, !line.empty() )
922 MF_LOG_WARNING(
"MARLEYTimeGen") <<
"Trailing content after last time" 923 <<
" bin found on line " << line_num <<
" of the spectrum file " 924 << full_spectrum_file_name;
928 double Emean_nue, alpha_nue, luminosity_nue;
929 double Emean_nuebar, alpha_nuebar, luminosity_nuebar;
930 double Emean_nux, alpha_nux, luminosity_nux;
931 std::istringstream iss(line);
932 bool ok_first =
static_cast<bool>( iss >> time );
935 <<
"Time bin left edges given in the spectrum file must be" 936 <<
" strictly increasing. Invalid time bin value found on line " 937 << line_num <<
" of the spectrum file " << full_spectrum_file_name;
938 else old_time = time;
940 bool ok_rest =
static_cast<bool>( iss >> Emean_nue >> alpha_nue
941 >> luminosity_nue >> Emean_nuebar >> alpha_nuebar >> luminosity_nuebar
942 >> Emean_nux >> alpha_nux >> luminosity_nux
949 fTimeFits.emplace_back(time, Emean_nue, alpha_nue, luminosity_nue,
950 Emean_nuebar, alpha_nuebar, luminosity_nuebar, Emean_nux,
951 alpha_nux, luminosity_nux);
958 else throw cet::exception(
"MARLEYTimeGen") <<
"Parse error on line " 959 << line_num <<
" of the spectrum file " << full_spectrum_file_name;
966 <<
"Missing right edge for the final time bin in the spectrum file " 967 << full_spectrum_file_name <<
". Unable to guess a bin width for the " 972 double delta_t_bin =
fTimeFits.back().Time()
973 -
fTimeFits.at(num_time_bins - 2).Time();
975 double last_bin_right_edge =
fTimeFits.back().Time() + delta_t_bin;
980 <<
" edge for the final time bin in the spectrum file " 981 << full_spectrum_file_name <<
". Assuming a width of " 982 << delta_t_bin <<
" s for the final bin.";
988 std::vector< std::unique_ptr< marley::NeutrinoSource > > fit_sources;
995 auto temp_source = std::make_unique<marley::FunctionNeutrinoSource>(
997 [&fit_sources,
this](
double E_nu) ->
double {
999 for (
size_t s = 0;
s < fit_sources.
size(); ++
s) {
1001 const TimeFit& current_fit = this->fTimeFits.at(
s);
1002 const auto& current_source = fit_sources.at(
s);
1003 const FitParameters&
params = current_fit.GetFitParameters(
1004 current_source->get_pid() );
1006 double lum = params.Luminosity();
1011 if (lum <= 0.)
continue;
1013 flux += lum * current_source->pdf(E_nu);
1019 double flux_integ = 0.;
1020 double tot_xs_integ = 0.;
1021 flux_averaged_total_xs(*temp_source, gen, flux_integ, tot_xs_integ);
1025 * tot_xs_integ / flux_integ;
1032 throw cet::exception(
"MARLEYTimeGen") <<
"Unrecognized neutrino spectrum" 1033 <<
" file format \"" <<
p().spectrum_file_format_() <<
"\" encountered" 1034 <<
" in evgen::MarleyTimeGen::reconfigure()";
1037 MF_LOG_INFO(
"MARLEYTimeGen") <<
"The flux-averaged total cross section" 1038 <<
" predicted by MARLEY for the current supernova spectrum is "
SpectrumFileFormat fSpectrumFileFormat
Format to assume for the neutrino spectrum input file.
PinchParamType fPinchType
The pinching parameter convention to use when interpreting the time-dependent fits.
std::unique_ptr< evgen::ActiveVolumeVertexSampler > fVertexSampler
Algorithm that allows us to sample vertex locations within the active volume(s) of the detector...
std::unique_ptr< evgen::MARLEYHelper > fMarleyHelper
Object that provides an interface to the MARLEY event generator.
double fFluxAveragedCrossSection
Flux-averaged total cross section (fm2, average is taken over all energies and times for all defined ...
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...
double fFitEmax
Maximum neutrino energy to consider when using a "fit"-format spectrum file.
#define MF_LOG_INFO(category)
std::unique_ptr< TH2D > fSpectrumHist
ROOT TH2D that contains the time-dependent spectrum to use when sampling neutrino times and energies...
unsigned int fNeutrinosPerEvent
The number of MARLEY neutrino vertices to generate in each art::Event.
void line(double t, double *p, double &x, double &y, double &z)
TimeGenSamplingMode fSamplingMode
Represents the sampling mode to use when selecting neutrino times and energies.
double fFitEmin
Minimum neutrino energy to consider when using a "fit"-format spectrum file.
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 ri...
std::vector< TimeFit > fTimeFits
Vector that contains the fit parameter information for each time bin when using a "fit"-format spectr...
#define MF_LOG_WARNING(category)
void make_nu_emission_histograms() const
Makes ROOT histograms showing the emitted neutrinos in each time bin when using a "fit"-format spectr...
cet::coded_exception< error, detail::translate > exception