Classes | Public Types | Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
evgen::MarleyTimeGen Class Reference
Inheritance diagram for evgen::MarleyTimeGen:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

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 Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 

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)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
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)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
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)
 
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...
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
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::MARLEYHelperfMarleyHelper
 Object that provides an interface to the MARLEY event generator. More...
 
std::unique_ptr< evgen::ActiveVolumeVertexSamplerfVertexSampler
 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< TimeFitfTimeFits
 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...
 

Additional Inherited Members

- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 

Detailed Description

Definition at line 111 of file MARLEYTimeGen_module.cc.

Member Typedef Documentation

Definition at line 116 of file MARLEYTimeGen_module.cc.

Definition at line 115 of file MARLEYTimeGen_module.cc.

Definition at line 226 of file MARLEYTimeGen_module.cc.

Member Enumeration Documentation

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

Enumerator
ALPHA 
BETA 

Definition at line 463 of file MARLEYTimeGen_module.cc.

463 { ALPHA, BETA };

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

Enumerator
RootTH2D 
FIT 

Definition at line 484 of file MARLEYTimeGen_module.cc.

484 { RootTH2D, FIT };

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

432 { HISTOGRAM, UNIFORM_TIME, UNIFORM_ENERGY };

Constructor & Destructor Documentation

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

Definition at line 528 of file MARLEYTimeGen_module.cc.

529  : EDProducer{ p.get_PSet() }, fEvent(new marley::Event), fRunNumber(0),
531 {
532  // Configure the module (including MARLEY itself) using the FHiCL parameters
533  this->reconfigure(p);
534 
535  // Create a ROOT TTree using the TFileService that will store the MARLEY
536  // event objects (useful for debugging purposes)
538  fEventTree = tfs->make<TTree>("MARLEY_event_tree",
539  "Neutrino events generated by MARLEY");
540  fEventTree->Branch("event", "marley::Event", fEvent.get());
541 
542  // Add branches that give the art::Event run, subrun, and event numbers for
543  // easy match-ups between the MARLEY and art TTrees. All three are recorded
544  // as 32-bit unsigned integers.
545  fEventTree->Branch("run_number", &fRunNumber, "run_number/i");
546  fEventTree->Branch("subrun_number", &fSubRunNumber, "subrun_number/i");
547  fEventTree->Branch("event_number", &fEventNumber, "event_number/i");
548  fEventTree->Branch("tSN", &fTNu, "tSN/D");
549  fEventTree->Branch("weight", &fWeight, "weight/D");
550 
551  produces< std::vector<simb::MCTruth> >();
552  produces< std::vector<sim::SupernovaTruth> >();
553  produces< art::Assns<simb::MCTruth, sim::SupernovaTruth> >();
554  produces< sumdata::RunData, art::InRun >();
555 }
uint_fast32_t fEventNumber
Event number from the art::Event being processed.
virtual void reconfigure(const Parameters &p)
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
double fFluxAveragedCrossSection
Flux-averaged total cross section (fm2, average is taken over all energies and times for all defined ...
double fWeight
Statistical weight for the current MARLEY neutrino vertex.
TTree * fEventTree
The event TTree created by MARLEY.
uint_fast32_t fSubRunNumber
Subrun number from the art::Event being processed.
uint_fast32_t fRunNumber
Run number from the art::Event being processed.
double fTNu
Time since the supernova core bounce for the current MARLEY neutrino vertex.
p
Definition: test.py:223
std::unique_ptr< marley::Event > fEvent
unique_ptr to the current event created by MARLEY

Member Function Documentation

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

Reimplemented from art::EDProducer.

Definition at line 558 of file MARLEYTimeGen_module.cc.

558  {
560  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
561 }
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
LArSoft geometry interface.
Definition: ChannelGeo.h:16
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 564 of file MARLEYTimeGen_module.cc.

566 {
567  // Get a reference to the generator object created by MARLEY (we'll need
568  // to do a few fancy things with it other than just creating events)
569  marley::Generator& gen = fMarleyHelper->get_generator();
570 
572  {
573  // Generate a MARLEY event using the time-integrated spectrum
574  // (the generator was already configured to use it by reconfigure())
575  mc_truth = fMarleyHelper->create_MCTruth(vertex_pos,
576  fEvent.get());
577 
578  // Find the time distribution corresponding to the selected energy bin
579  double E_nu = fEvent->projectile().total_energy();
580  int E_bin_index = fSpectrumHist->GetYaxis()->FindBin(E_nu);
581  TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist", E_bin_index,
582  E_bin_index);
583  double* time_bin_weights = t_hist->GetArray();
584 
585  // Sample a time bin from the distribution
586  std::discrete_distribution<int> time_dist;
587  std::discrete_distribution<int>::param_type time_params(time_bin_weights,
588  &(time_bin_weights[t_hist->GetNbinsX() + 1]));
589  int time_bin_index = gen.sample_from_distribution(time_dist, time_params);
590 
591  // Sample a time uniformly from within the selected time bin
592  double t_min = t_hist->GetBinLowEdge(time_bin_index);
593  double t_max = t_min + t_hist->GetBinWidth(time_bin_index);
594  // sample a time on [ t_min, t_max )
595  fTNu = gen.uniform_random_double(t_min, t_max, false);
596  // Unbiased sampling was used, so assign this neutrino vertex a
597  // unit statistical weight
598  fWeight = ONE;
599 
602  }
603 
605  {
606  // Generate a MARLEY event using the time-integrated spectrum
607  // (the generator was already configured to use it by reconfigure())
608  mc_truth = fMarleyHelper->create_MCTruth(vertex_pos,
609  fEvent.get());
610 
611  // Sample a time uniformly
612  TAxis* time_axis = fSpectrumHist->GetXaxis();
613  // underflow bin has index zero
614  double t_min = time_axis->GetBinLowEdge(1);
615  double t_max = time_axis->GetBinLowEdge(time_axis->GetNbins() + 1);
616  // sample a time on [ t_min, t_max )
617  fTNu = gen.uniform_random_double(t_min, t_max, false);
618 
619  // Get the value of the true dependent probability density (probability
620  // of the sampled time given the sampled energy) to use as a biasing
621  // correction in the neutrino vertex weight.
622  double E_nu = fEvent->projectile().total_energy();
623  int E_bin_index = fSpectrumHist->GetYaxis()->FindBin(E_nu);
624  TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist", E_bin_index,
625  E_bin_index);
626  int t_bin_index = t_hist->FindBin(fTNu);
627  double weight_bias = t_hist->GetBinContent(t_bin_index) * (t_max - t_min)
628  / ( t_hist->Integral() * t_hist->GetBinWidth(t_bin_index) );
629 
630  fWeight = weight_bias;
631 
634  }
635 
637  {
638  // Select a time bin using the energy-integrated spectrum
639  TH1D* t_hist = fSpectrumHist->ProjectionX("dummy_time_hist");
640  double* time_bin_weights = t_hist->GetArray();
641 
642  // Sample a time bin from the distribution
643  std::discrete_distribution<int> time_dist;
644  std::discrete_distribution<int>::param_type time_params(time_bin_weights,
645  &(time_bin_weights[t_hist->GetNbinsX() + 1]));
646  int time_bin_index = gen.sample_from_distribution(time_dist, time_params);
647 
648  // Sample a time uniformly from within the selected time bin
649  double t_min = t_hist->GetBinLowEdge(time_bin_index);
650  double t_max = t_min + t_hist->GetBinWidth(time_bin_index);
651  // sample a time on [ t_min, t_max )
652  fTNu = gen.uniform_random_double(t_min, t_max, false);
653 
654  // Sample an energy uniformly over the entire allowed range
655  // underflow bin has index zero
656  TAxis* energy_axis = fSpectrumHist->GetYaxis();
657  double E_min = energy_axis->GetBinLowEdge(1);
658  double E_max = energy_axis->GetBinLowEdge(energy_axis->GetNbins() + 1);
659 
660  double E_nu = std::numeric_limits<double>::lowest();
661 
662  // Generate a MARLEY event using a uniformly sampled energy
663  mc_truth = make_uniform_energy_mctruth(E_min, E_max, E_nu, vertex_pos);
664 
665  // Get the value of the true dependent probability density (probability
666  // of the sampled energy given the sampled time) to use as a biasing
667  // correction in the neutrino vertex weight.
668  //
669  // Get a 1D projection of the energy spectrum for the sampled time bin
670  TH1D* energy_spect = fSpectrumHist->ProjectionY("energy_spect",
671  time_bin_index, time_bin_index);
672 
673  // Create a new MARLEY neutrino source object using this projection (this
674  // will create a normalized probability density that we can use) and load
675  // it into the generator.
676  auto nu_source = marley_root::make_root_neutrino_source(
677  marley_utils::ELECTRON_NEUTRINO, energy_spect);
678  double new_source_E_min = nu_source->get_Emin();
679  double new_source_E_max = nu_source->get_Emax();
680  gen.set_source(std::move(nu_source));
681  // NOTE: The marley::Generator object normalizes the E_pdf to unity
682  // automatically, but just in case, we redo it here.
683  double E_pdf_integ = integrate([&gen](double E_nu)
684  -> double { return gen.E_pdf(E_nu); }, new_source_E_min,
685  new_source_E_max);
686 
687  // Compute the likelihood ratio that we need to bias the neutrino vertex
688  // weight
689  double weight_bias = (gen.E_pdf(E_nu) / E_pdf_integ) * (E_max - E_min);
690 
691  fWeight = weight_bias;
692 
695  }
696 
697  else {
698  throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
699  << " encountered in evgen::MarleyTimeGen::produce()";
700  }
701 
702 }
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 ...
double fWeight
Statistical weight for the current MARLEY neutrino vertex.
Arrival times were sampled uniformly.
def move(depos, offset)
Definition: depos.py:107
double fTNu
Time since the supernova core bounce for the current MARLEY neutrino vertex.
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.
gen
Definition: demo.py:24
std::unique_ptr< TH2D > fSpectrumHist
ROOT TH2D that contains the time-dependent spectrum to use when sampling neutrino times and energies...
std::unique_ptr< marley::Event > fEvent
unique_ptr to the current event created by MARLEY
Sample directly from cross-section weighted spectrum.
TimeGenSamplingMode fSamplingMode
Represents the sampling mode to use when selecting neutrino times and energies.
Neutrino energies were sampled uniformly.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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 1044 of file MARLEYTimeGen_module.cc.

1046 {
1047  // Get a reference to the generator object created by MARLEY (we'll need
1048  // to do a few fancy things with it other than just creating events)
1049  marley::Generator& gen = fMarleyHelper->get_generator();
1050 
1051  // Initialize the time bin index to something absurdly large. This will help
1052  // us detect strange bugs that arise when it is sampled incorrectly.
1053  size_t time_bin_index = std::numeric_limits<size_t>::max();
1054 
1055  // Create an object to represent the discrete time distribution given in
1056  // the spectrum file. Use the luminosity for each bin as its sampling weight.
1057  // This distribution will actually be used to sample a neutrino arrival time
1058  // unless we're using the uniform time sampling mode, in which case it will
1059  // be used to help calculate the neutrino vertex weight.
1060  int source_pdg_code = gen.get_source().get_pid();
1061 
1062  const auto fit_params_begin = TimeFit::make_FitParameters_iterator(
1063  source_pdg_code, fTimeFits.begin() );
1064  const auto fit_params_end = TimeFit::make_FitParameters_iterator(
1065  source_pdg_code, fTimeFits.end() );
1066 
1067  const auto lum_begin = FitParameters::make_luminosity_iterator(
1068  fit_params_begin);
1069  const auto lum_end = FitParameters::make_luminosity_iterator(
1070  fit_params_end);
1071 
1072  std::discrete_distribution<size_t> time_dist(lum_begin, lum_end);
1073 
1076  {
1077  time_bin_index = gen.sample_from_distribution(time_dist);
1078  }
1079 
1081  int last_time_index = 0;
1082  if (fTimeFits.size() > 0) last_time_index = fTimeFits.size() - 1;
1083  std::uniform_int_distribution<size_t> uid(0, last_time_index);
1084  // for c2: time_bin_index has already been declared
1085  time_bin_index = gen.sample_from_distribution(uid);
1086  }
1087 
1088  else {
1089  throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
1090  << " encountered in evgen::MarleyTimeGen::produce()";
1091  }
1092 
1093  // Sample a time uniformly from within the selected time bin. Note that
1094  // the entries in fTimeFits use the time_ member to store the bin left
1095  // edges. The module creates fTimeFits in such a way that its last element
1096  // will always have luminosity_ == 0. (zero sampling weight), so we may
1097  // always add one to the sampled bin index without worrying about going
1098  // off the edge.
1099  double t_min = fTimeFits.at(time_bin_index).Time();
1100  double t_max = fTimeFits.at(time_bin_index + 1).Time();
1101  // sample a time on [ t_min, t_max )
1102  fTNu = gen.uniform_random_double(t_min, t_max, false);
1103 
1104  // Create a "beta-fit" neutrino source using the correct parameters for the
1105  // sampled time bin. This will be used to sample a neutrino energy unless
1106  // we're using the uniform time sampling mode. For uniform time sampling,
1107  // it will be used to determine the neutrino event weight.
1108  const auto& fit = fTimeFits.at(time_bin_index);
1109  std::unique_ptr<marley::NeutrinoSource> nu_source = source_from_time_fit(fit);
1110 
1113  {
1114  // Replace the generator's old source with the new one for the current
1115  // time bin
1116  gen.set_source(std::move(nu_source));
1117 
1118  // Generate a MARLEY event using the updated source
1119  mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());
1120 
1122  // Unbiased sampling creates neutrino vertices with unit weight
1123  fWeight = ONE;
1125  sim::kUnbiased);
1126  }
1127  else {
1128  // fSamplingMode == TimeGenSamplingMode::UNIFORM_TIME
1129 
1130  // Multiply by the likelihood ratio in order to correct for uniform
1131  // time sampling if we're using that biasing method
1132  double weight_bias = time_dist.probabilities().at(time_bin_index)
1133  / (t_max - t_min) * ( fTimeFits.back().Time()
1134  - fTimeFits.front().Time() );
1135 
1136  fWeight = weight_bias;
1139  }
1140  }
1141 
1143  {
1144  double E_nu = std::numeric_limits<double>::lowest();
1145  mc_truth = make_uniform_energy_mctruth(fFitEmin, fFitEmax, E_nu,
1146  vertex_pos);
1147 
1148  // Get the value of the true dependent probability density (probability
1149  // of the sampled energy given the sampled time) to use as a biasing
1150  // correction in the neutrino vertex weight.
1151 
1152  // Load the generator with the neutrino source that represents the
1153  // true (i.e., unbiased) energy probability distribution. This will
1154  // create a normalized probability density that we can use to determine
1155  // the neutrino vertex weight.
1156  double nu_source_E_min = nu_source->get_Emin();
1157  double nu_source_E_max = nu_source->get_Emax();
1158  gen.set_source(std::move(nu_source));
1159 
1160  // NOTE: The marley::Generator object normalizes the E_pdf to unity
1161  // automatically, but just in case, we redo it here.
1162  double E_pdf_integ = integrate([&gen](double E_nu)
1163  -> double { return gen.E_pdf(E_nu); }, nu_source_E_min,
1164  nu_source_E_max);
1165 
1166  // Compute the likelihood ratio that we need to bias the neutrino vertex
1167  // weight
1168  double weight_bias = (gen.E_pdf(E_nu) / E_pdf_integ)
1169  * (fFitEmax - fFitEmin);
1170 
1171  fWeight = weight_bias;
1174  }
1175 
1176  else {
1177  throw cet::exception("MARLEYTimeGen") << "Unrecognized sampling mode"
1178  << " encountered in evgen::MarleyTimeGen::produce()";
1179  }
1180 
1181 }
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 ...
double fWeight
Statistical weight for the current MARLEY neutrino vertex.
Arrival times were sampled uniformly.
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...
static marley::IteratorToMember< It, double > make_luminosity_iterator(It it)
Converts an iterator that points to a FitParameters object into an iterator that points to that objec...
def move(depos, offset)
Definition: depos.py:107
double fTNu
Time since the supernova core bounce for the current MARLEY neutrino vertex.
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.
gen
Definition: demo.py:24
double fFitEmax
Maximum neutrino energy to consider when using a "fit"-format spectrum file.
static marley::IteratorToMember< It, FitParameters > make_FitParameters_iterator(int pdg_code, It iterator)
Converts an iterator that points to a TimeFit object into an iterator that points to the appropriate ...
static int max(int a, int b)
std::unique_ptr< marley::Event > fEvent
unique_ptr to the current event created by MARLEY
Sample directly from cross-section weighted spectrum.
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.
std::vector< TimeFit > fTimeFits
Vector that contains the fit parameter information for each time bin when using a "fit"-format spectr...
Neutrino energies were sampled uniformly.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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 1259 of file MARLEYTimeGen_module.cc.

1260 {
1261  // The final time bin has zero luminosity, and therefore zero sampling
1262  // weight. We need it to be present so that the last nonzero weight bin
1263  // has a right edge.
1264  fTimeFits.emplace_back(time, 0., 0., 0., 0., 0., 0., 0., 0., 0.);
1265 }
std::vector< TimeFit > fTimeFits
Vector that contains the fit parameter information for each time bin when using a "fit"-format spectr...
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 1268 of file MARLEYTimeGen_module.cc.

1269 {
1270  // To make a histogram with variable size bins, ROOT needs the bin
1271  // low edges to be contiguous in memory. This is not true of the
1272  // stored times in fTimeFits, so we'll need to make a temporary copy
1273  // of them.
1274  std::vector<double> temp_times;
1275  for (const auto& fit : fTimeFits) temp_times.push_back(fit.Time());
1276 
1277  // If, for some reason, there are fewer than two time bins, just return
1278  // without making the histograms.
1279  // TODO: consider throwing an exception here
1280  if (temp_times.size() < 2) return;
1281 
1282  // Get the number of time bins
1283  int num_bins = temp_times.size() - 1;
1284 
1285  // Create some ROOT TH1D objects using the TFileService. These will store
1286  // the number of emitted neutrinos of each type in each time bin
1288  TH1D* nue_hist = tfs->make<TH1D>("NueEmission",
1289  "Number of emitted #nu_{e}; time (s); number of neutrinos in time bin",
1290  num_bins, temp_times.data());
1291  TH1D* nuebar_hist = tfs->make<TH1D>("NuebarEmission",
1292  "Number of emitted #bar{#nu}_{e}; time (s); number of neutrinos in"
1293  " time bin", num_bins, temp_times.data());
1294  TH1D* nux_hist = tfs->make<TH1D>("NuxEmission",
1295  "Number of emitted #nu_{x}; time (s); number of neutrinos in time bin",
1296  num_bins, temp_times.data());
1297 
1298  // Load the histograms with the emitted neutrino counts
1299  for (size_t b = 1; b < temp_times.size(); ++b) {
1300  const TimeFit& current_fit = fTimeFits.at(b - 1);
1301  const TimeFit& next_fit = fTimeFits.at(b);
1302  double bin_deltaT = next_fit.Time() - current_fit.Time();
1303 
1304  const auto& nue_params = current_fit.GetFitParameters(
1305  marley_utils::ELECTRON_NEUTRINO);
1306  const auto& nuebar_params = current_fit.GetFitParameters(
1307  marley_utils::ELECTRON_ANTINEUTRINO);
1308  const auto& nux_params = current_fit.GetFitParameters(
1309  marley_utils::MUON_NEUTRINO);
1310 
1311  // Convert from bin luminosity to number of neutrinos by
1312  // multiplying by the bin time interval and dividing by the
1313  // mean neutrino energy
1314  double num_nue = 0.;
1315  double num_nuebar = 0.;
1316  double num_nux = 0.;
1317  if (nue_params.Emean() != 0.) num_nue = nue_params.Luminosity()
1318  * bin_deltaT / (nue_params.Emean() * MeV_to_erg);
1319  if (nuebar_params.Emean() != 0.) num_nuebar = nuebar_params.Luminosity()
1320  * bin_deltaT / (nuebar_params.Emean() * MeV_to_erg);
1321  if (nux_params.Emean() != 0.) num_nux = nux_params.Luminosity()
1322  * bin_deltaT / (nux_params.Emean() * MeV_to_erg);
1323 
1324  nue_hist->SetBinContent(b, num_nue);
1325  nuebar_hist->SetBinContent(b, num_nuebar);
1326  nux_hist->SetBinContent(b, num_nux);
1327  }
1328 
1329 }
static bool * b
Definition: config.cpp:1043
std::vector< TimeFit > fTimeFits
Vector that contains the fit parameter information for each time bin when using a "fit"-format spectr...
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 1212 of file MARLEYTimeGen_module.cc.

1214 {
1215  marley::Generator& gen = fMarleyHelper->get_generator();
1216 
1217  // Sample an energy uniformly over the entire allowed range
1218  double total_xs;
1219  int j = 0;
1220  do {
1221  // We have to check that the cross section is nonzero for the sampled
1222  // energy (otherwise we'll generate an unphysical event). However, if the
1223  // user has given us a histogram that is below threshold, the
1224  // program could get stuck here endlessly, sampling rejected energy
1225  // after rejected energy. Just in case, we cap the total number of tries
1226  // and quit if things don't work out.
1227  if (j >= MAX_UNIFORM_ENERGY_ITERATIONS) {
1228  throw cet::exception("MARLEYTimeGen") << "Exceeded the maximum of "
1229  << MAX_UNIFORM_ENERGY_ITERATIONS << " while attempting to sample"
1230  << " a neutrino energy uniformly.";
1231  }
1232  // Sample an energy uniformly on [ E_min, E_max )
1233  E_nu = gen.uniform_random_double(E_min, E_max, false);
1234  total_xs = 0.;
1235  // Check that at least one defined reaction has a non-zero total
1236  // cross section at the sampled energy. If this is not the case, try
1237  // again.
1238  for (const auto& react : gen.get_reactions()) {
1239  total_xs += react->total_xs(marley_utils::ELECTRON_NEUTRINO, E_nu);
1240  }
1241 
1242  ++j;
1243  } while (total_xs <= 0.);
1244 
1245  // Replace the existing neutrino source with a monoenergetic one at the
1246  // neutrino energy that was just sampled above
1247  std::unique_ptr<marley::NeutrinoSource> nu_source
1248  = std::make_unique<marley::MonoNeutrinoSource>(
1249  marley_utils::ELECTRON_NEUTRINO, E_nu);
1250  gen.set_source(std::move(nu_source));
1251 
1252  // Generate a MARLEY event using the new monoenergetic source
1253  auto mc_truth = fMarleyHelper->create_MCTruth(vertex_pos, fEvent.get());
1254 
1255  return mc_truth;
1256 }
std::unique_ptr< evgen::MARLEYHelper > fMarleyHelper
Object that provides an interface to the MARLEY event generator.
def move(depos, offset)
Definition: depos.py:107
gen
Definition: demo.py:24
std::unique_ptr< marley::Event > fEvent
unique_ptr to the current event created by MARLEY
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgen::MarleyTimeGen::produce ( art::Event e)
overridevirtual

Implements art::EDProducer.

Definition at line 705 of file MARLEYTimeGen_module.cc.

706 {
708 
709  // Get the run, subrun, and event numbers from the current art::Event
710  fRunNumber = e.run();
711  fSubRunNumber = e.subRun();
712  fEventNumber = e.event();
713 
714  // Prepare associations and vectors of truth objects that will be produced
715  // and loaded into the current art::Event
716  std::unique_ptr< std::vector<simb::MCTruth> >
717  truthcol(new std::vector<simb::MCTruth>);
718 
719  std::unique_ptr< std::vector<sim::SupernovaTruth> >
720  sn_truthcol(new std::vector<sim::SupernovaTruth>);
721 
722  std::unique_ptr< art::Assns<simb::MCTruth, sim::SupernovaTruth> >
724 
725  // Create temporary truth objects that we will use to load the event
726  simb::MCTruth truth;
727  sim::SupernovaTruth sn_truth;
728 
729  for (unsigned int n = 0; n < fNeutrinosPerEvent; ++n) {
730 
731  // Sample a primary vertex location for this event
732  TLorentzVector vertex_pos = fVertexSampler->sample_vertex_pos(*geo);
733 
734  // Reset the neutrino's time-since-supernova to a bogus value (for now)
735  fTNu = std::numeric_limits<double>::lowest();
736 
738  create_truths_th2d(truth, sn_truth, vertex_pos);
739  }
741  create_truths_time_fit(truth, sn_truth, vertex_pos);
742  }
743  else {
744  throw cet::exception("MARLEYTimeGen") << "Invalid spectrum file"
745  << " format encountered in evgen::MarleyTimeGen::produce()";
746  }
747 
748  // Write the marley::Event object to the event tree
749  fEventTree->Fill();
750 
751  // Add the truth objects to the appropriate vectors
752  truthcol->push_back(truth);
753 
754  sn_truthcol->push_back(sn_truth);
755 
756  // Associate the last entries in each of the truth object vectors (the
757  // truth objects that were just created for the current neutrino vertex)
758  // with each other
759  truth_assns->addSingle(art::PtrMaker<simb::MCTruth>{e}(truthcol->size() - 1),
760  art::PtrMaker<sim::SupernovaTruth>{e}(sn_truthcol->size() - 1));
761  }
762 
763  // Load the completed truth object vectors and associations into the event
764  e.put(std::move(truthcol));
765 
766  e.put(std::move(sn_truthcol));
767 
768  e.put(std::move(truth_assns));
769 }
EventNumber_t event() const
Definition: DataViewImpl.cc:85
SpectrumFileFormat fSpectrumFileFormat
Format to assume for the neutrino spectrum input file.
uint_fast32_t fEventNumber
Event number from the art::Event being processed.
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 previou...
std::unique_ptr< evgen::ActiveVolumeVertexSampler > fVertexSampler
Algorithm that allows us to sample vertex locations within the active volume(s) of the detector...
TTree * fEventTree
The event TTree created by MARLEY.
uint_fast32_t fSubRunNumber
Subrun number from the art::Event being processed.
uint_fast32_t fRunNumber
Run number from the art::Event being processed.
std::void_t< T > n
def move(depos, offset)
Definition: depos.py:107
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...
double fTNu
Time since the supernova core bounce for the current MARLEY neutrino vertex.
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
RunNumber_t run() const
Definition: DataViewImpl.cc:71
unsigned int fNeutrinosPerEvent
The number of MARLEY neutrino vertices to generate in each art::Event.
Event generator information.
Definition: MCTruth.h:32
LArSoft geometry interface.
Definition: ChannelGeo.h:16
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgen::MarleyTimeGen::reconfigure ( const Parameters p)
virtual

Definition at line 772 of file MARLEYTimeGen_module.cc.

773 {
774  const auto& seed_service = art::ServiceHandle<rndm::NuRandomService>();
775  const auto& geom_service = art::ServiceHandle<geo::Geometry const>();
776 
777  // Create a new evgen::ActiveVolumeVertexSampler object based on the current
778  // configuration
779  fVertexSampler = std::make_unique<evgen::ActiveVolumeVertexSampler>(
780  p().vertex_, *seed_service, *geom_service, "MARLEY_Vertex_Sampler");
781 
782  // Create a new marley::Generator object based on the current configuration
783  const fhicl::ParameterSet marley_pset = p.get_PSet()
784  .get< fhicl::ParameterSet >( "marley_parameters" );
785  fMarleyHelper = std::make_unique<MARLEYHelper>( marley_pset,
786  *seed_service, "MARLEY" );
787 
788  // Get the number of neutrino vertices per event from the FHiCL parameters
789  fNeutrinosPerEvent = p().nu_per_event_();
790 
791  // Determine the current sampling mode from the FHiCL parameters
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")
799  else throw cet::exception("MARLEYTimeGen")
800  << "Invalid sampling mode \"" << samp_mode_str << "\""
801  << " specified for the MARLEYTimeGen module.";
802 
803  MF_LOG_INFO("MARLEYTimeGen") << fNeutrinosPerEvent << " neutrino vertices"
804  << " will be generated for each art::Event using the \"" << samp_mode_str
805  << "\" sampling mode.";
806 
807  // Retrieve the time-dependent neutrino spectrum from the spectrum file.
808  // Use different methods depending on the file's format.
809  std::string spectrum_file_format = marley_utils::to_lowercase(
810  p().spectrum_file_format_());
811 
812  if (spectrum_file_format == "th2d")
814  else if (spectrum_file_format == "fit") {
816 
817  std::string pinch_type;
818  if ( !p().pinching_parameter_type_(pinch_type) ) {
819  throw cet::exception("MARLEYTimeGen") << "Missing pinching parameter"
820  << " type for a \"fit\" format spectrum file";
821  }
822 
823  marley_utils::to_lowercase_inplace(pinch_type);
824  if (pinch_type == "alpha") fPinchType = PinchParamType::ALPHA;
825  else if (pinch_type == "beta") fPinchType = PinchParamType::BETA;
826  else throw cet::exception("MARLEYTimeGen")
827  << "Invalid pinching parameter type \"" << pinch_type
828  << "\" specified for the MARLEYTimeGen module.";
829 
830  if ( !p().fit_Emin_(fFitEmin) ) throw cet::exception("MARLEYTimeGen")
831  << "Missing minimum energy for a \"fit\" format spectrum"
832  << " used by the MARLEYTimeGen module.";
833 
834  if ( !p().fit_Emax_(fFitEmax) ) throw cet::exception("MARLEYTimeGen")
835  << "Missing maximum energy for a \"fit\" format spectrum"
836  << " used by the MARLEYTimeGen module.";
837 
838  if (fFitEmax < fFitEmin) throw cet::exception("MARLEYTimeGen")
839  << "Maximum energy is less than the minimum energy for"
840  << " a \"fit\" format spectrum used by the MARLEYTimeGen module.";
841  }
842  else throw cet::exception("MARLEYTimeGen")
843  << "Invalid spectrum file format \"" << p().spectrum_file_format_()
844  << "\" specified for the MARLEYTimeGen module.";
845 
846  // Determine the full file name (including path) of the spectrum file
847  std::string full_spectrum_file_name
848  = fMarleyHelper->find_file(p().spectrum_file_(), "spectrum");
849 
850  marley::Generator& gen = fMarleyHelper->get_generator();
851 
853 
854  // Retrieve the time-dependent neutrino flux from a ROOT file
855  std::unique_ptr<TFile> spectrum_file
856  = std::make_unique<TFile>(full_spectrum_file_name.c_str(), "read");
857  TH2D* temp_h2 = nullptr;
858  std::string namecycle;
859  if ( !p().namecycle_(namecycle) ) {
860  throw cet::exception("MARLEYTimeGen") << "Missing namecycle for"
861  << " a TH2D spectrum file";
862  }
863 
864  spectrum_file->GetObject(namecycle.c_str(), temp_h2);
865  fSpectrumHist.reset(temp_h2);
866 
867  // Disassociate the TH2D from its parent TFile. If we fail to do this,
868  // then ROOT will auto-delete the TH2D when the TFile goes out of scope.
869  fSpectrumHist->SetDirectory(nullptr);
870 
871  // Compute the flux-averaged total cross section using MARLEY. This will be
872  // used to compute neutrino vertex weights for the sim::SupernovaTruth
873  // objects.
874 
875  // Get a 1D projection of the energy spectrum (integrated over time)
876  TH1D* energy_spect = fSpectrumHist->ProjectionY("energy_spect");
877 
878  // Create a new MARLEY neutrino source object using this projection
879  // TODO: replace the hard-coded electron neutrino PDG code here (and in
880  // several other places in this source file) when you're ready to use
881  // MARLEY with multiple neutrino flavors
882  std::unique_ptr<marley::NeutrinoSource> nu_source
883  = marley_root::make_root_neutrino_source(marley_utils::ELECTRON_NEUTRINO,
884  energy_spect);
885 
886  // Factor of hbar_c^2 converts from MeV^(-2) to fm^2
887  fFluxAveragedCrossSection = marley_utils::hbar_c2
888  * flux_averaged_total_xs(*nu_source, gen);
889 
890  // For speed, sample energies first whenever possible (and then sample from
891  // an energy-dependent timing distribution). This avoids unnecessary calls
892  // to MARLEY to change the energy spectrum.
895  {
896  gen.set_source(std::move(nu_source));
897  }
898 
899  } // spectrum_file_format == "th2d"
900 
902 
903  // Clear out the old parameterized spectrum, if one exists
904  fTimeFits.clear();
905 
906  std::ifstream fit_file(full_spectrum_file_name);
908 
909  bool found_end = false;
910 
911  // current line number
912  int line_num = 0;
913  // number of lines checked in last call to marley_utils::get_next_line()
914  int lines_checked = 0;
915 
916  double old_time = std::numeric_limits<double>::lowest();
917 
918  while ( line = marley_utils::get_next_line(fit_file, rx_comment_or_empty,
919  false, lines_checked), line_num += lines_checked, !line.empty() )
920  {
921  if (found_end) {
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;
925  }
926 
927  double time;
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 );
933 
934  if (time <= old_time) throw cet::exception("MARLEYTimeGen")
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;
939 
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
943  );
944 
945  if (ok_first) {
946  // We haven't reached the final bin, so add another time bin
947  // in the typical way.
948  if (ok_rest) {
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);
952  }
953  else {
954  make_final_timefit(time);
955  found_end = true;
956  }
957  }
958  else throw cet::exception("MARLEYTimeGen") << "Parse error on line "
959  << line_num << " of the spectrum file " << full_spectrum_file_name;
960  }
961 
962  if (!found_end) {
963 
964  size_t num_time_bins = fTimeFits.size();
965  if (num_time_bins < 2) throw cet::exception("MARLEYTimeGen")
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 "
968  << " final bin.";
969 
970  // Guess that the width of the penultimate time bin and the width of
971  // the final time bin are the same
972  double delta_t_bin = fTimeFits.back().Time()
973  - fTimeFits.at(num_time_bins - 2).Time();
974 
975  double last_bin_right_edge = fTimeFits.back().Time() + delta_t_bin;
976 
977  make_final_timefit(last_bin_right_edge);
978 
979  MF_LOG_WARNING("MARLEYTimeGen") << "Missing right"
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.";
983  }
984 
985 
986  // Compute the flux-averaged total cross section for the fitted spectrum.
987  // We will need this to compute neutrino vertex weights.
988  std::vector< std::unique_ptr< marley::NeutrinoSource > > fit_sources;
989  for (const auto& fit : fTimeFits) {
990  fit_sources.emplace_back( source_from_time_fit(fit) );
991  }
992 
993  // TODO: add handling for non-nu_e neutrino types when suitable data become
994  // available in MARLEY
995  auto temp_source = std::make_unique<marley::FunctionNeutrinoSource>(
996  marley_utils::ELECTRON_NEUTRINO, fFitEmin, fFitEmax,
997  [&fit_sources, this](double E_nu) -> double {
998  double flux = 0.;
999  for (size_t s = 0; s < fit_sources.size(); ++s) {
1000 
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() );
1005 
1006  double lum = params.Luminosity();
1007 
1008  // Skip entries with zero luminosity, since they won't contribute
1009  // anything to the overall integral. Skip negative luminosity ones as
1010  // well, just in case.
1011  if (lum <= 0.) continue;
1012 
1013  flux += lum * current_source->pdf(E_nu);
1014  }
1015  return flux;
1016  }
1017  );
1018 
1019  double flux_integ = 0.;
1020  double tot_xs_integ = 0.;
1021  flux_averaged_total_xs(*temp_source, gen, flux_integ, tot_xs_integ);
1022 
1023  // Factor of hbar_c^2 converts from MeV^(-2) to fm^2
1024  fFluxAveragedCrossSection = marley_utils::hbar_c2
1025  * tot_xs_integ / flux_integ;
1026 
1028 
1029  } // spectrum_file_format == "fit"
1030 
1031  else {
1032  throw cet::exception("MARLEYTimeGen") << "Unrecognized neutrino spectrum"
1033  << " file format \"" << p().spectrum_file_format_() << "\" encountered"
1034  << " in evgen::MarleyTimeGen::reconfigure()";
1035  }
1036 
1037  MF_LOG_INFO("MARLEYTimeGen") << "The flux-averaged total cross section"
1038  << " predicted by MARLEY for the current supernova spectrum is "
1039  << fFluxAveragedCrossSection << " fm^2";
1040 
1041 }
SpectrumFileFormat fSpectrumFileFormat
Format to assume for the neutrino spectrum input file.
std::string string
Definition: nybbler.cc:12
PinchParamType fPinchType
The pinching parameter convention to use when interpreting the time-dependent fits.
uint size() const
Definition: qcstring.h:201
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...
def move(depos, offset)
Definition: depos.py:107
p
Definition: test.py:223
gen
Definition: demo.py:24
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...
static QCString * s
Definition: config.cpp:1042
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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 1185 of file MARLEYTimeGen_module.cc.

1186 {
1187  // Create a "beta-fit" neutrino source using the given fit parameters.
1188 
1189  // The two common fitting schemes (alpha and beta) differ in their
1190  // definitions by \beta = \alpha + 1.
1191  // TODO: add handling for non-nu_e neutrino types
1192  const FitParameters& nue_parameters
1193  = fit.GetFitParameters(marley_utils::ELECTRON_NEUTRINO);
1194 
1195  double beta = nue_parameters.Alpha();
1196  if (fPinchType == PinchParamType::ALPHA) beta += 1.;
1197  else if (fPinchType != PinchParamType::BETA) {
1198  throw cet::exception("MARLEYTimeGen") << "Unreognized pinching parameter"
1199  << " type encountered in evgen::MarleyTimeGen::source_from_time_fit()";
1200  }
1201 
1202  // Create the new source
1203  std::unique_ptr<marley::NeutrinoSource> nu_source
1204  = std::make_unique<marley::BetaFitNeutrinoSource>(
1205  marley_utils::ELECTRON_NEUTRINO, fFitEmin, fFitEmax,
1206  nue_parameters.Emean(), beta);
1207 
1208  return nu_source;
1209 }
double beta(double KE, const simb::MCParticle *part)
PinchParamType fPinchType
The pinching parameter convention to use when interpreting the time-dependent fits.
double fFitEmax
Maximum neutrino energy to consider when using a "fit"-format spectrum file.
double fFitEmin
Minimum neutrino energy to consider when using a "fit"-format spectrum file.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

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

unique_ptr to the current event created by MARLEY

Definition at line 403 of file MARLEYTimeGen_module.cc.

uint_fast32_t evgen::MarleyTimeGen::fEventNumber
protected

Event number from the art::Event being processed.

Definition at line 500 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 493 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 524 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 520 of file MARLEYTimeGen_module.cc.

double evgen::MarleyTimeGen::fFluxAveragedCrossSection
protected

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.

Definition at line 512 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 396 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 516 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 467 of file MARLEYTimeGen_module.cc.

uint_fast32_t evgen::MarleyTimeGen::fRunNumber
protected

Run number from the art::Event being processed.

Definition at line 496 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 436 of file MARLEYTimeGen_module.cc.

SpectrumFileFormat evgen::MarleyTimeGen::fSpectrumFileFormat
protected

Format to assume for the neutrino spectrum input file.

Definition at line 487 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 409 of file MARLEYTimeGen_module.cc.

uint_fast32_t evgen::MarleyTimeGen::fSubRunNumber
protected

Subrun number from the art::Event being processed.

Definition at line 498 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 415 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 504 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 400 of file MARLEYTimeGen_module.cc.

double evgen::MarleyTimeGen::fWeight
protected

Statistical weight for the current MARLEY neutrino vertex.

Definition at line 507 of file MARLEYTimeGen_module.cc.


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