diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh index e0995f47a4..44a03facce 100644 --- a/EventGenerator/inc/ParticleGeneratorTool.hh +++ b/EventGenerator/inc/ParticleGeneratorTool.hh @@ -9,11 +9,13 @@ #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" #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 +45,45 @@ namespace mu2e { virtual ~ParticleGeneratorTool() noexcept = default; + double calculateBinnedSpectrumEnergyFraction(fhicl::ParameterSet pset, + 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 + 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); + } + 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 + const double xmin = fullSpectrum.getXMin(); + if(correct_full_integral && xmin > full_var_low) { + const double pdfmin = fullSpectrum.getPDF(0); + const double binsize = fullSpectrum.getBinWidth(); + fullIntegral += 0.5*(pdfmin/binsize)*(xmin - full_var_low); // 0 to over x = to linear interpolation + } + + // 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..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: @@ -48,29 +50,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, 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; + 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 +79,7 @@ namespace mu2e { BinnedSpectrum _spectrum; double _emin; double _emax; - double _energy_fraction; + double _energyFraction; bool _flatSpectrum; std::unique_ptr _randomUnitSphere; @@ -105,7 +90,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()); @@ -135,7 +120,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/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 042d93c977..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: @@ -40,9 +42,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 = (_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; } std::vector generate() override; @@ -52,7 +59,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()); @@ -69,6 +76,7 @@ namespace mu2e { double _czMax; double _spectrumXMin; double _spectrumXMax; + double _energyFraction; bool _flatSpectrum; std::unique_ptr _randomPoissonQ; @@ -121,7 +129,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 @@ + + 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