From fd5c2a789e6e3ed6b82f65048d4ce0ec144658a9 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Fri, 10 Jul 2026 13:53:02 -0500 Subject: [PATCH 1/7] Move binned spectrum fraction calculation, allow neutron tail sim --- EventGenerator/inc/ParticleGeneratorTool.hh | 33 +++++++++++++++++++ EventGenerator/src/DIOGenerator_tool.cc | 33 +++++-------------- .../src/MuCapNeutronGenerator_tool.cc | 8 ++++- MCDataProducts/src/classes_def.xml | 2 ++ 4 files changed, 50 insertions(+), 26 deletions(-) diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh index e0995f47a4..6b8684459c 100644 --- a/EventGenerator/inc/ParticleGeneratorTool.hh +++ b/EventGenerator/inc/ParticleGeneratorTool.hh @@ -14,6 +14,7 @@ #include "Offline/MCDataProducts/inc/ProcessCode.hh" #include "Offline/MCDataProducts/inc/GenParticle.hh" #include "Offline/MCDataProducts/inc/SpectrumConfig.hh" +#include "Offline/Mu2eUtilities/inc/BinnedSpectrum.hh" #include "Offline/GeneralUtilities/inc/RSNTIO.hh" @@ -43,6 +44,38 @@ namespace mu2e { virtual ~ParticleGeneratorTool() noexcept = default; + double calculateBinnedSpectrumEnergyFraction(fhicl::ParameterSet pset, + std::string var_low = "elow", // default to energy spectrum variables + std::string var_high = "ehi", + double full_var_low = 0., // if elow == ehi, it does the full spectrum + double full_var_high = 0.) const { + + // Initialize the spectra with and without the (possible) energy restriction + BinnedSpectrum spectrum(pset); + pset.erase(var_low); + pset.erase(var_high); + pset.put(var_low, full_var_low); + pset.put(var_high, full_var_high); + BinnedSpectrum fullSpectrum(pset); + + // Calculate the integrals + double integral = 0.; + double fullIntegral = 0.; + for(size_t ibin=0;ibin < spectrum.getNbins();++ibin){ + integral += spectrum.getPDF(ibin); + } + + // The "full" integral may omit the 0 - first bin center contribution + for(size_t ibin=0;ibin < fullSpectrum.getNbins();++ibin){ + fullIntegral += fullSpectrum.getPDF(ibin); + } + + // Evaluate the fraction being sampled + const double fraction = (fullIntegral > 0.) ? integral / fullIntegral : 0.; + + return fraction; + } + bool _isPrimary = true; // flag to indicate if this is for primary generation or not }; } diff --git a/EventGenerator/src/DIOGenerator_tool.cc b/EventGenerator/src/DIOGenerator_tool.cc index 756d3a0555..9e6126820a 100644 --- a/EventGenerator/src/DIOGenerator_tool.cc +++ b/EventGenerator/src/DIOGenerator_tool.cc @@ -48,29 +48,12 @@ namespace mu2e { auto fullconfig = conf().spectrum.get(); _emin = fullconfig.get("elow", _spectrum.getXMin()); _emax = fullconfig.get("ehi", _spectrum.getXMax()); - fullconfig.erase(std::string("elow")); - fullconfig.erase(std::string("ehi")); - fullconfig.put(std::string("elow"),double(0.0)); - fullconfig.put(std::string("ehi"),double(0.0)); - BinnedSpectrum fullspect(fullconfig); - double fullintegral(0.0); - for(size_t ibin=0;ibin < fullspect.getNbins();++ibin){ - fullintegral += fullspect.getPDF(ibin); - } - // correct for the missing prediction near threshold, assuming the rate falls to 0 linearly. - double pmin = _spectrum.getAbscissa(0); - double pdfmin = _spectrum.getPDF(0); - double binsize = _spectrum.getBinWidth(); - fullintegral += 0.5*pdfmin*pmin/binsize; - _energy_fraction = (fullintegral > 0.) ? integral / fullintegral : 0.; - std::cout << "Cos(theta_z) min " << _czmin << " max " << _czmax << std::endl; - std::cout << "Restricted Spectrum min " << _spectrum.getAbscissa(0) << " max " << _spectrum.getAbscissa(_spectrum.getNbins()-1) << std::endl; - std::cout << "Full Spectrum min " << fullspect.getAbscissa(0) << " max " << fullspect.getAbscissa(fullspect.getNbins()-1) << std::endl; - std::cout << "Restricted Spectrum integral " << integral << std::endl; - std::cout << "Restricted Spectrum integral*cos(theta_z) restriction " << integral*((_czmax - _czmin)/2.) << std::endl; - std::cout << "Full Spectrum integral " << fullintegral << std::endl; - std::cout << "Sampled spectrum fraction " << integral/fullintegral << std::endl; - std::cout << "Sampled spectrum fraction (with cos(theta_z)) " << (integral/fullintegral)*((_czmax - _czmin)/2.) << std::endl; + _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig); + + std::cout << "[" << __func__ << "] Cos(theta_z) min " << _czmin << " max " << _czmax << std::endl; + std::cout << "[" << __func__ << "] Restricted Spectrum min " << _spectrum.getAbscissa(0) << " max " << _spectrum.getAbscissa(_spectrum.getNbins()-1) << std::endl; + std::cout << "[" << __func__ << "] Sampled spectrum fraction " << _energyFraction << std::endl; + std::cout << "[" << __func__ << "] Sampled spectrum fraction (with cos(theta_z)) " << (_energyFraction)*((_czmax - _czmin)/2.) << std::endl; } @@ -94,7 +77,7 @@ namespace mu2e { BinnedSpectrum _spectrum; double _emin; double _emax; - double _energy_fraction; + double _energyFraction; bool _flatSpectrum; std::unique_ptr _randomUnitSphere; @@ -135,7 +118,7 @@ namespace mu2e { std::unique_ptr DIOGenerator::spectrumConfig() { auto config = std::make_unique(); - config->add_var(SpectrumConfig::RestrictedVar("energy", _energy_fraction , _emin , _emax, + config->add_var(SpectrumConfig::RestrictedVar("energy", _energyFraction , _emin , _emax, _flatSpectrum ? SpectrumConfig::Type::kFlat : SpectrumConfig::Type::kPhysical)); config->add_var(SpectrumConfig::RestrictedVar("cosz" , (_czmax - _czmin)/2., _czmin, _czmax)); return config; diff --git a/EventGenerator/src/MuCapNeutronGenerator_tool.cc b/EventGenerator/src/MuCapNeutronGenerator_tool.cc index 042d93c977..a131b65555 100644 --- a/EventGenerator/src/MuCapNeutronGenerator_tool.cc +++ b/EventGenerator/src/MuCapNeutronGenerator_tool.cc @@ -40,9 +40,14 @@ namespace mu2e { _czMax(conf().czMax()), _spectrumXMin(_spectrum.getXMin()), _spectrumXMax(_spectrum.getXMax()), + _energyFraction(1.), _flatSpectrum(conf().spectrum.get().get("spectrumShape", "") == "flat") { if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax; + auto fullconfig = conf().spectrum.get(); + _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig); + std::cout << "[" << __func__ << "] Sampled spectrum fraction " << _energyFraction << std::endl; + std::cout << "[" << __func__ << "] Sampled spectrum fraction (with cos(theta_z)) " << (_energyFraction)*((_czMax - _czMin)/2.) << std::endl; } std::vector generate() override; @@ -69,6 +74,7 @@ namespace mu2e { double _czMax; double _spectrumXMin; double _spectrumXMax; + double _energyFraction; bool _flatSpectrum; std::unique_ptr _randomPoissonQ; @@ -121,7 +127,7 @@ namespace mu2e { default: spectrumVarName = "energy"; break; } - config->add_var(SpectrumConfig::RestrictedVar(spectrumVarName, 1., _spectrumXMin, _spectrumXMax, + config->add_var(SpectrumConfig::RestrictedVar(spectrumVarName, _energyFraction, _spectrumXMin, _spectrumXMax, _flatSpectrum ? SpectrumConfig::Type::kFlat : SpectrumConfig::Type::kPhysical)); config->add_var(SpectrumConfig::RestrictedVar("cosz", (_czMax - _czMin)/2., _czMin, _czMax)); return config; diff --git a/MCDataProducts/src/classes_def.xml b/MCDataProducts/src/classes_def.xml index 1bae40fc57..350fc1d0d6 100644 --- a/MCDataProducts/src/classes_def.xml +++ b/MCDataProducts/src/classes_def.xml @@ -40,6 +40,8 @@ + + From 27cdacdd06969978b8b0acb18f84959249fd4f32 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Fri, 10 Jul 2026 15:38:02 -0500 Subject: [PATCH 2/7] Reduce rates for non-primary muon capture products if energy is restricted --- EventGenerator/src/DIOGenerator_tool.cc | 2 +- EventGenerator/src/MuCapNeutronGenerator_tool.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/EventGenerator/src/DIOGenerator_tool.cc b/EventGenerator/src/DIOGenerator_tool.cc index 9e6126820a..5a8326973c 100644 --- a/EventGenerator/src/DIOGenerator_tool.cc +++ b/EventGenerator/src/DIOGenerator_tool.cc @@ -88,7 +88,7 @@ namespace mu2e { std::vector DIOGenerator::generate() { std::vector res; - const double r = (_czmax - _czmin)/2.; + const double r = _energyFraction*(_czmax - _czmin)/2.; // reduce the rate by the spectrum restriction if(_isPrimary || _randFlat->fire() <= r) { double energy = _spectrum.sample(_randSpectrum->fire()); diff --git a/EventGenerator/src/MuCapNeutronGenerator_tool.cc b/EventGenerator/src/MuCapNeutronGenerator_tool.cc index a131b65555..8e80206db3 100644 --- a/EventGenerator/src/MuCapNeutronGenerator_tool.cc +++ b/EventGenerator/src/MuCapNeutronGenerator_tool.cc @@ -57,7 +57,7 @@ namespace mu2e { void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override { _isPrimary = isPrimary; _rate = GlobalConstantsHandle()->getCaptureNeutronRate(material); - const double rate = _rate * (_czMax - _czMin)/2.; // accouunt for potential cz selection in the produced rates + const double rate = _rate * _energyFraction * (_czMax - _czMin)/2.; // account for potential spectrum restriction in the produced rates _randomUnitSphere = std::make_unique(eng, _czMin, _czMax); _randomPoissonQ = std::make_unique(eng, rate); _randSpectrum = std::make_unique(eng, _spectrum.getPDF(), _spectrum.getNbins()); From 81db65d35954346dfe744236b3e79c2df311d65c Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Sat, 11 Jul 2026 12:00:08 -0500 Subject: [PATCH 3/7] Add option to "correct" spectrum full integral and use for DIO --- EventGenerator/inc/ParticleGeneratorTool.hh | 11 +++++++++-- EventGenerator/src/DIOGenerator_tool.cc | 2 +- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh index 6b8684459c..f772b5590f 100644 --- a/EventGenerator/inc/ParticleGeneratorTool.hh +++ b/EventGenerator/inc/ParticleGeneratorTool.hh @@ -45,6 +45,7 @@ namespace mu2e { virtual ~ParticleGeneratorTool() noexcept = default; double calculateBinnedSpectrumEnergyFraction(fhicl::ParameterSet pset, + const bool correct_full_integral = false, // correct for missing low tail to full integral std::string var_low = "elow", // default to energy spectrum variables std::string var_high = "ehi", double full_var_low = 0., // if elow == ehi, it does the full spectrum @@ -64,11 +65,17 @@ namespace mu2e { for(size_t ibin=0;ibin < spectrum.getNbins();++ibin){ integral += spectrum.getPDF(ibin); } - - // The "full" integral may omit the 0 - first bin center contribution for(size_t ibin=0;ibin < fullSpectrum.getNbins();++ibin){ fullIntegral += fullSpectrum.getPDF(ibin); } + // The "full" integral may be missing the lowest tail component in some cases (e.g. DIO) + // --> apply a linear interpolation correction if requested + if(correct_full_integral) { + const double xmin = fullSpectrum.getAbscissa(0); + const double pdfmin = fullSpectrum.getPDF(0); + const double binsize = fullSpectrum.getBinWidth(); + fullIntegral += 0.5*pdfmin*xmin/binsize; // 0 to from x = 0 to x = linear interpolation + } // Evaluate the fraction being sampled const double fraction = (fullIntegral > 0.) ? integral / fullIntegral : 0.; diff --git a/EventGenerator/src/DIOGenerator_tool.cc b/EventGenerator/src/DIOGenerator_tool.cc index 5a8326973c..fc1c5d1efb 100644 --- a/EventGenerator/src/DIOGenerator_tool.cc +++ b/EventGenerator/src/DIOGenerator_tool.cc @@ -48,7 +48,7 @@ namespace mu2e { auto fullconfig = conf().spectrum.get(); _emin = fullconfig.get("elow", _spectrum.getXMin()); _emax = fullconfig.get("ehi", _spectrum.getXMax()); - _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig); + _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig, true); // energy fraction, correcting for missing 0 - 1 MeV component std::cout << "[" << __func__ << "] Cos(theta_z) min " << _czmin << " max " << _czmax << std::endl; std::cout << "[" << __func__ << "] Restricted Spectrum min " << _spectrum.getAbscissa(0) << " max " << _spectrum.getAbscissa(_spectrum.getNbins()-1) << std::endl; From 977788c4b0c57fc7f17e1973c186ec77fbfef936 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Sun, 12 Jul 2026 11:38:28 -0500 Subject: [PATCH 4/7] Linear interpolation correction by default is missing PDF range --- EventGenerator/inc/ParticleGeneratorTool.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh index f772b5590f..daa10a3912 100644 --- a/EventGenerator/inc/ParticleGeneratorTool.hh +++ b/EventGenerator/inc/ParticleGeneratorTool.hh @@ -45,7 +45,7 @@ namespace mu2e { virtual ~ParticleGeneratorTool() noexcept = default; double calculateBinnedSpectrumEnergyFraction(fhicl::ParameterSet pset, - const bool correct_full_integral = false, // correct for missing low tail to full integral + const bool correct_full_integral = true, // correct for missing low tail to full integral std::string var_low = "elow", // default to energy spectrum variables std::string var_high = "ehi", double full_var_low = 0., // if elow == ehi, it does the full spectrum @@ -70,11 +70,11 @@ namespace mu2e { } // The "full" integral may be missing the lowest tail component in some cases (e.g. DIO) // --> apply a linear interpolation correction if requested - if(correct_full_integral) { - const double xmin = fullSpectrum.getAbscissa(0); + const double xmin = fullSpectrum.getAbscissa(0); + if(correct_full_integral && xmin > full_var_low) { const double pdfmin = fullSpectrum.getPDF(0); const double binsize = fullSpectrum.getBinWidth(); - fullIntegral += 0.5*pdfmin*xmin/binsize; // 0 to from x = 0 to x = linear interpolation + fullIntegral += 0.5*(pdfmin/binsize)*(xmin - full_var_low); // 0 to over x = to linear interpolation } // Evaluate the fraction being sampled From 8a0737ac611cc40c1cc62aa2274bbb8bd3e9f874 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 16 Jul 2026 10:10:19 -0500 Subject: [PATCH 5/7] Preserve previous behavior --- EventGenerator/src/MuCapNeutronGenerator_tool.cc | 2 +- Mu2eG4/src/CompositeMaterialGenerator_module.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/EventGenerator/src/MuCapNeutronGenerator_tool.cc b/EventGenerator/src/MuCapNeutronGenerator_tool.cc index 8e80206db3..7fa15d811e 100644 --- a/EventGenerator/src/MuCapNeutronGenerator_tool.cc +++ b/EventGenerator/src/MuCapNeutronGenerator_tool.cc @@ -45,7 +45,7 @@ namespace mu2e { { if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax; auto fullconfig = conf().spectrum.get(); - _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig); + _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig, false); std::cout << "[" << __func__ << "] Sampled spectrum fraction " << _energyFraction << std::endl; std::cout << "[" << __func__ << "] Sampled spectrum fraction (with cos(theta_z)) " << (_energyFraction)*((_czMax - _czMin)/2.) << std::endl; } diff --git a/Mu2eG4/src/CompositeMaterialGenerator_module.cc b/Mu2eG4/src/CompositeMaterialGenerator_module.cc index 773f914a8b..3a0cc4220d 100644 --- a/Mu2eG4/src/CompositeMaterialGenerator_module.cc +++ b/Mu2eG4/src/CompositeMaterialGenerator_module.cc @@ -122,7 +122,7 @@ namespace mu2e{ auto generator_config = element.generator_tool.get(); auto generator_tool = art::make_tool(generator_config); // dummy argument because material is not actually needed in this context - generator_tool->finishInitialization(_engine, "", false); + generator_tool->finishInitialization(_engine, "", true); // enforce that keys to subgenerators are unique // would be better in a protected class, but eh From 1813c6adfd92d92a0e951ed65bc838be3ec33e19 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 16 Jul 2026 10:14:38 -0500 Subject: [PATCH 6/7] Fix xmin value in tail integral approx --- EventGenerator/inc/ParticleGeneratorTool.hh | 2 +- EventGenerator/src/MuCapNeutronGenerator_tool.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh index daa10a3912..41e3cb79be 100644 --- a/EventGenerator/inc/ParticleGeneratorTool.hh +++ b/EventGenerator/inc/ParticleGeneratorTool.hh @@ -70,7 +70,7 @@ namespace mu2e { } // The "full" integral may be missing the lowest tail component in some cases (e.g. DIO) // --> apply a linear interpolation correction if requested - const double xmin = fullSpectrum.getAbscissa(0); + const double xmin = fullSpectrum.getXMin(); if(correct_full_integral && xmin > full_var_low) { const double pdfmin = fullSpectrum.getPDF(0); const double binsize = fullSpectrum.getBinWidth(); diff --git a/EventGenerator/src/MuCapNeutronGenerator_tool.cc b/EventGenerator/src/MuCapNeutronGenerator_tool.cc index 7fa15d811e..8e80206db3 100644 --- a/EventGenerator/src/MuCapNeutronGenerator_tool.cc +++ b/EventGenerator/src/MuCapNeutronGenerator_tool.cc @@ -45,7 +45,7 @@ namespace mu2e { { if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax; auto fullconfig = conf().spectrum.get(); - _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig, false); + _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig); std::cout << "[" << __func__ << "] Sampled spectrum fraction " << _energyFraction << std::endl; std::cout << "[" << __func__ << "] Sampled spectrum fraction (with cos(theta_z)) " << (_energyFraction)*((_czMax - _czMin)/2.) << std::endl; } From 8653622fde0c9d34b3d36b586117a71033629b7c Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 16 Jul 2026 10:30:34 -0500 Subject: [PATCH 7/7] Correct integrals, address comments --- EventGenerator/inc/ParticleGeneratorTool.hh | 1 + EventGenerator/src/DIOGenerator_tool.cc | 2 + EventGenerator/src/Mu2eXGenerator_tool.cc | 40 +++++-------------- .../src/MuCapNeutronGenerator_tool.cc | 4 +- 4 files changed, 17 insertions(+), 30 deletions(-) diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh index 41e3cb79be..44a03facce 100644 --- a/EventGenerator/inc/ParticleGeneratorTool.hh +++ b/EventGenerator/inc/ParticleGeneratorTool.hh @@ -9,6 +9,7 @@ #include "art/Utilities/ToolConfigTable.h" #include "art/Framework/Services/Optional/RandomNumberGenerator.h" +#include "fhiclcpp/ParameterSet.h" #include "Offline/DataProducts/inc/PDGCode.hh" #include "Offline/MCDataProducts/inc/ProcessCode.hh" diff --git a/EventGenerator/src/DIOGenerator_tool.cc b/EventGenerator/src/DIOGenerator_tool.cc index fc1c5d1efb..7ea67201da 100644 --- a/EventGenerator/src/DIOGenerator_tool.cc +++ b/EventGenerator/src/DIOGenerator_tool.cc @@ -16,6 +16,8 @@ #include "fhiclcpp/types/DelegatedParameter.h" +#include + namespace mu2e { class DIOGenerator : public ParticleGeneratorTool { public: diff --git a/EventGenerator/src/Mu2eXGenerator_tool.cc b/EventGenerator/src/Mu2eXGenerator_tool.cc index fc92ff2259..8300986087 100644 --- a/EventGenerator/src/Mu2eXGenerator_tool.cc +++ b/EventGenerator/src/Mu2eXGenerator_tool.cc @@ -16,6 +16,8 @@ #include "fhiclcpp/types/DelegatedParameter.h" +#include + namespace mu2e { class Mu2eXGenerator : public ParticleGeneratorTool { public: @@ -33,39 +35,19 @@ namespace mu2e { _spectrum(BinnedSpectrum(conf().spectrum.get())), _emin(0.), _emax(0.), - _energy_fraction(1.), + _energyFraction(1.), _flatSpectrum(conf().spectrum.get().get("spectrumShape", "") == "flat") { - // compute normalization - double integral(0.0); - for(size_t ibin=0;ibin < _spectrum.getNbins();++ibin){ - integral += _spectrum.getPDF(ibin); - } - _emin = _spectrum.getXMin(); _emax = _spectrum.getXMax(); auto fullconfig = conf().spectrum.get(); - fullconfig.erase(std::string("elow")); - fullconfig.erase(std::string("ehi")); - fullconfig.put(std::string("elow"),double(0.0)); - fullconfig.put(std::string("ehi"),double(0.0)); - BinnedSpectrum fullspect(fullconfig); - double fullintegral(0.0); - for(size_t ibin=0;ibin < fullspect.getNbins();++ibin){ - fullintegral += fullspect.getPDF(ibin); - } - // correct for the missing prediction near threshold, assuming the rate falls to 0 linearly. - double pmin = _spectrum.getAbscissa(0); - double pdfmin = _spectrum.getPDF(0); - double binsize = _spectrum.getBinWidth(); - fullintegral += 0.5*pdfmin*pmin/binsize; - _energy_fraction = (fullintegral > 0.) ? integral/fullintegral : 0.; - std::cout << "Restricted Spectrum min " << _spectrum.getAbscissa(0) << " max " << _spectrum.getAbscissa(_spectrum.getNbins()-1) << std::endl; - std::cout << "Full Spectrum min " << fullspect.getAbscissa(0) << " max " << fullspect.getAbscissa(fullspect.getNbins()-1) << std::endl; - std::cout << "Restricted Spectrum integral " << integral << std::endl; - std::cout << "Full Spectrum integral " << fullintegral << std::endl; - std::cout << "Sampled spectrum fraction " << _energy_fraction << std::endl; + _emin = fullconfig.get("elow", _spectrum.getXMin()); + _emax = fullconfig.get("ehi", _spectrum.getXMax()); + _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig, true); // energy fraction evaluation + + std::cout << "[" << __func__ << "] Restricted Spectrum min " << _spectrum.getAbscissa(0) << " max " << _spectrum.getAbscissa(_spectrum.getNbins()-1) << std::endl; + std::cout << "[" << __func__ << "] Sampled spectrum fraction " << _energyFraction << std::endl; } @@ -86,7 +68,7 @@ namespace mu2e { BinnedSpectrum _spectrum; double _emin; double _emax; - double _energy_fraction; + double _energyFraction; bool _flatSpectrum; std::unique_ptr _randomUnitSphere; @@ -123,7 +105,7 @@ namespace mu2e { std::unique_ptr Mu2eXGenerator::spectrumConfig() { auto config = std::make_unique(); - config->add_var(SpectrumConfig::RestrictedVar("energy", _energy_fraction, _emin, _emax, + config->add_var(SpectrumConfig::RestrictedVar("energy", _energyFraction, _emin, _emax, _flatSpectrum ? SpectrumConfig::Type::kFlat : SpectrumConfig::Type::kPhysical)); return config; } diff --git a/EventGenerator/src/MuCapNeutronGenerator_tool.cc b/EventGenerator/src/MuCapNeutronGenerator_tool.cc index 8e80206db3..8659691b6c 100644 --- a/EventGenerator/src/MuCapNeutronGenerator_tool.cc +++ b/EventGenerator/src/MuCapNeutronGenerator_tool.cc @@ -17,6 +17,8 @@ #include "fhiclcpp/types/DelegatedParameter.h" +#include + namespace mu2e { class MuCapNeutronGenerator : public ParticleGeneratorTool { public: @@ -45,7 +47,7 @@ namespace mu2e { { if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax; auto fullconfig = conf().spectrum.get(); - _energyFraction = calculateBinnedSpectrumEnergyFraction(fullconfig); + _energyFraction = (_flatSpectrum) ? 1. : calculateBinnedSpectrumEnergyFraction(fullconfig, false); std::cout << "[" << __func__ << "] Sampled spectrum fraction " << _energyFraction << std::endl; std::cout << "[" << __func__ << "] Sampled spectrum fraction (with cos(theta_z)) " << (_energyFraction)*((_czMax - _czMin)/2.) << std::endl; }