Neutron generator updates

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 <pdfmin> over x = <spectrum low edge> to <xmin> 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
   };
 }

EventGenerator/src/DIOGenerator_tool.cc

@@ -16,6 +16,8 @@
 
 #include "fhiclcpp/types/DelegatedParameter.h"
 
+#include <iostream>
+
 namespace mu2e {
   class DIOGenerator : public ParticleGeneratorTool {
   public:
@@ -48,29 +50,12 @@ namespace mu2e {
       auto fullconfig = conf().spectrum.get<fhicl::ParameterSet>();
       _emin = fullconfig.get<double>("elow", _spectrum.getXMin());
       _emax = fullconfig.get<double>("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>   _randomUnitSphere;
@@ -105,7 +90,7 @@ namespace mu2e {
 
   std::vector<ParticleGeneratorTool::Kinematic> DIOGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  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<SpectrumConfig> DIOGenerator::spectrumConfig() {
     auto config = std::make_unique<SpectrumConfig>();
-    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;

EventGenerator/src/Mu2eXGenerator_tool.cc

@@ -16,6 +16,8 @@
 
 #include "fhiclcpp/types/DelegatedParameter.h"
 
+#include <iostream>
+
 namespace mu2e {
   class Mu2eXGenerator : public ParticleGeneratorTool {
   public:
@@ -33,39 +35,19 @@ namespace mu2e {
       _spectrum(BinnedSpectrum(conf().spectrum.get<fhicl::ParameterSet>())),
       _emin(0.),
       _emax(0.),
-      _energy_fraction(1.),
+      _energyFraction(1.),
       _flatSpectrum(conf().spectrum.get<fhicl::ParameterSet>().get<std::string>("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<fhicl::ParameterSet>();
-      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<double>("elow", _spectrum.getXMin());
+      _emax = fullconfig.get<double>("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>   _randomUnitSphere;
@@ -123,7 +105,7 @@ namespace mu2e {
 
   std::unique_ptr<SpectrumConfig> Mu2eXGenerator::spectrumConfig() {
     auto config = std::make_unique<SpectrumConfig>();
-    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;
   }

EventGenerator/src/MuCapNeutronGenerator_tool.cc

@@ -17,6 +17,8 @@
 
 #include "fhiclcpp/types/DelegatedParameter.h"
 
+#include <iostream>
+
 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<fhicl::ParameterSet>().get<std::string>("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<fhicl::ParameterSet>();
+      _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<ParticleGeneratorTool::Kinematic> 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<PhysicsParams>()->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<RandomUnitSphere>(eng, _czMin, _czMax);
       _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, rate);
       _randSpectrum = std::make_unique<CLHEP::RandGeneral>(eng, _spectrum.getPDF(), _spectrum.getNbins());
@@ -69,6 +76,7 @@ namespace mu2e {
     double _czMax;
     double _spectrumXMin;
     double _spectrumXMax;
+    double _energyFraction;
     bool _flatSpectrum;
 
     std::unique_ptr<CLHEP::RandPoissonQ> _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;

MCDataProducts/src/classes_def.xml

@@ -40,6 +40,8 @@
  <class name="art::Wrapper<mu2e::SpectrumConfig::RestrictedVar>"/>
  <class name="std::map<std::string,mu2e::SpectrumConfig::RestrictedVar>"/>
  <class name="art::Wrapper<std::map<std::string,mu2e::SpectrumConfig::RestrictedVar>>"/>
+ <class name="std::pair<std::string,mu2e::SpectrumConfig::RestrictedVar>"/>
+ <class name="art::Wrapper<std::pair<std::string,mu2e::SpectrumConfig::RestrictedVar>>"/>
 
 <!--  ********* simulation  ********* -->
  <class name="mu2e::StatusG4"/>

Mu2eG4/src/CompositeMaterialGenerator_module.cc

@@ -122,7 +122,7 @@ namespace mu2e{
       auto generator_config = element.generator_tool.get<fhicl::ParameterSet>();
       auto generator_tool = art::make_tool<ParticleGeneratorTool>(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