diff --git a/EventGenerator/inc/ParticleGeneratorTool.hh b/EventGenerator/inc/ParticleGeneratorTool.hh
index 3b28fd2c4c..d983500fd7 100644
--- a/EventGenerator/inc/ParticleGeneratorTool.hh
+++ b/EventGenerator/inc/ParticleGeneratorTool.hh
@@ -27,7 +27,8 @@ namespace mu2e {
       CLHEP::HepLorentzVector fourmom;
     };
 
-    virtual void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& materialName) = 0;
+    virtual void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& materialName,
+                                      const bool isPrimary) = 0;
 
     virtual std::vector<Kinematic> generate() = 0;
 
@@ -35,6 +36,8 @@ namespace mu2e {
     virtual void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) = 0;
 
     virtual ~ParticleGeneratorTool() noexcept = default;
+
+    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 d609ac5aa5..b63d4bf771 100644
--- a/EventGenerator/src/DIOGenerator_tool.cc
+++ b/EventGenerator/src/DIOGenerator_tool.cc
@@ -1,8 +1,8 @@
 #include "art/Utilities/ToolMacros.h"
 #include "cetlib_except/exception.h"
 
-#include "CLHEP/Random/RandPoissonQ.h"
 #include "CLHEP/Random/RandGeneral.h"
+#include "CLHEP/Random/RandFlat.h"
 
 #include "Offline/EventGenerator/inc/ParticleGeneratorTool.hh"
 
@@ -73,9 +73,11 @@ namespace mu2e {
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string&) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string&, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng, _czmin, _czmax);
       _randSpectrum = std::make_unique<CLHEP::RandGeneral>(eng, _spectrum.getPDF(), _spectrum.getNbins());
+      _randFlat = std::make_unique<CLHEP::RandFlat>(eng);
     }
 
   private:
@@ -88,20 +90,24 @@ namespace mu2e {
 
     std::unique_ptr<RandomUnitSphere>   _randomUnitSphere;
     std::unique_ptr<CLHEP::RandGeneral> _randSpectrum;
+    std::unique_ptr<CLHEP::RandFlat>    _randFlat;
   };
 
 
   std::vector<ParticleGeneratorTool::Kinematic> DIOGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  res;
+    const double r = (_czmax - _czmin)/2.;
+    if(_isPrimary || _randFlat->fire() <= r) {
 
-    double energy = _spectrum.sample(_randSpectrum->fire());
+      double energy = _spectrum.sample(_randSpectrum->fire());
 
-    const double p = energy * sqrt(1 - std::pow(_mass/energy,2));
-    CLHEP::Hep3Vector p3 = _randomUnitSphere->fire(p);
-    CLHEP::HepLorentzVector fourmom(p3, energy);
+      const double p = energy * sqrt(1 - std::pow(_mass/energy,2));
+      CLHEP::Hep3Vector p3 = _randomUnitSphere->fire(p);
+      CLHEP::HepLorentzVector fourmom(p3, energy);
 
-    ParticleGeneratorTool::Kinematic k{_pdgId, ProcessCode::mu2eMuonDecayAtRest, fourmom};
-    res.emplace_back(k);
+      ParticleGeneratorTool::Kinematic k{_pdgId, ProcessCode::mu2eMuonDecayAtRest, fourmom};
+      res.emplace_back(k);
+    }
 
     return res;
   }
diff --git a/EventGenerator/src/Mu2eXGenerator_tool.cc b/EventGenerator/src/Mu2eXGenerator_tool.cc
index 833017bd53..2e840a178b 100644
--- a/EventGenerator/src/Mu2eXGenerator_tool.cc
+++ b/EventGenerator/src/Mu2eXGenerator_tool.cc
@@ -64,7 +64,8 @@ namespace mu2e {
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string&) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string&, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng);
       _randSpectrum = std::make_unique<CLHEP::RandGeneral>(eng, _spectrum.getPDF(), _spectrum.getNbins());
     }
diff --git a/EventGenerator/src/MuCap1809keVGammaGenerator_tool.cc b/EventGenerator/src/MuCap1809keVGammaGenerator_tool.cc
index d0e534c8ac..9f1b38bd30 100644
--- a/EventGenerator/src/MuCap1809keVGammaGenerator_tool.cc
+++ b/EventGenerator/src/MuCap1809keVGammaGenerator_tool.cc
@@ -20,14 +20,12 @@ namespace mu2e {
       using Comment=fhicl::Comment;
       fhicl::Atom<double> czMin  {Name("czmin")  , Comment("Restrict cos(theta_z) minimum"), -1.};
       fhicl::Atom<double> czMax  {Name("czmax")  , Comment("Restrict cos(theta_z) maximum"),  1.};
-      fhicl::Atom<bool>   fireAll{Name("fireAll"), Comment("Add a photon to all events, otherwise use the branching fraction"),  false};
     };
     typedef art::ToolConfigTable<PhysConfig> Parameters;
 
     explicit MuCap1809keVGammaGenerator(Parameters const& conf) :
       _czMin(conf().czMin()),
       _czMax(conf().czMax()),
-      _fireAll(conf().fireAll()),
       _pdgId(PDGCode::gamma),
       _mass(GlobalConstantsHandle<ParticleDataList>()->particle(_pdgId).mass()),
       _randomUnitSphere(nullptr),
@@ -39,7 +37,8 @@ namespace mu2e {
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _energy = GlobalConstantsHandle<PhysicsParams>()->get1809keVGammaEnergy(material);
       _intensity = GlobalConstantsHandle<PhysicsParams>()->get1809keVGammaIntensity(material);
       _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng, _czMin, _czMax);
@@ -49,7 +48,6 @@ namespace mu2e {
   private:
     double _czMin;
     double _czMax;
-    bool   _fireAll;
     PDGCode::type _pdgId;
     double _mass;
     double _energy = 0.;
@@ -62,7 +60,8 @@ namespace mu2e {
   std::vector<ParticleGeneratorTool::Kinematic> MuCap1809keVGammaGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  res;
 
-    if (_fireAll || _randFlat->fire() < _intensity) {
+    const double intensity = _intensity * (_czMax - _czMin)/2.; // account for potential cz selection in the intensity
+    if (_isPrimary || _randFlat->fire() < intensity) {
       const double momentum = _energy * sqrt(1 - std::pow(_mass/_energy,2));
       CLHEP::Hep3Vector p3 = _randomUnitSphere->fire(momentum);
       CLHEP::HepLorentzVector fourmom(p3, _energy);
diff --git a/EventGenerator/src/MuCapDeuteronGenerator_tool.cc b/EventGenerator/src/MuCapDeuteronGenerator_tool.cc
index b48081e583..d2ec2dfc07 100644
--- a/EventGenerator/src/MuCapDeuteronGenerator_tool.cc
+++ b/EventGenerator/src/MuCapDeuteronGenerator_tool.cc
@@ -1,4 +1,5 @@
 #include "art/Utilities/ToolMacros.h"
+#include "cetlib_except/exception.h"
 
 #include "CLHEP/Random/RandPoissonQ.h"
 #include "CLHEP/Random/RandGeneral.h"
@@ -25,6 +26,9 @@ namespace mu2e {
 
       fhicl::DelegatedParameter spectrum{Name("spectrum"), Comment("Parameters for BinnedSpectrum)")};
       fhicl::Atom<std::string> spectrumVariable{Name("spectrumVariable"), Comment("Variable that spectrum is defined in (\"totalEnergy\", \"kineticEnergy\", or \"momentum\")")};
+      fhicl::Atom<double> czMin{Name("czmin"), Comment("Restrict cos(theta_z) minimum"), -1.};
+      fhicl::Atom<double> czMax{Name("czmax"), Comment("Restrict cos(theta_z) maximum"),  1.};
+
     };
     typedef art::ToolConfigTable<PhysConfig> Parameters;
 
@@ -32,16 +36,22 @@ namespace mu2e {
       _pdgId(PDGCode::deuteron),
       _mass(GlobalConstantsHandle<ParticleDataList>()->particle(_pdgId).mass()),
       _spectrum(BinnedSpectrum(conf().spectrum.get<fhicl::ParameterSet>())),
-      _spectrumVariable(parseSpectrumVar(conf().spectrumVariable()))
-    {}
+      _spectrumVariable(parseSpectrumVar(conf().spectrumVariable())),
+      _czMin(conf().czMin()),
+      _czMax(conf().czMax())
+    {
+      if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax;
+    }
 
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _rate = GlobalConstantsHandle<PhysicsParams>()->getCaptureDeuteronRate(material);
-      _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng);
-      _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, _rate);
+      const double rate = _rate * (_czMax - _czMin)/2.; // accouunt for potential cz selection 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());
     }
 
@@ -52,6 +62,8 @@ namespace mu2e {
 
     BinnedSpectrum    _spectrum;
     SpectrumVar       _spectrumVariable;
+    double _czMin;
+    double _czMax;
 
     std::unique_ptr<CLHEP::RandPoissonQ> _randomPoissonQ;
     std::unique_ptr<RandomUnitSphere>    _randomUnitSphere;
@@ -62,7 +74,7 @@ namespace mu2e {
   std::vector<ParticleGeneratorTool::Kinematic> MuCapDeuteronGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  res;
 
-    int n_gen = _randomPoissonQ->fire();
+    const int n_gen = (_isPrimary) ? 1 : _randomPoissonQ->fire();
     for (int i_gen = 0; i_gen < n_gen; ++i_gen) {
       double energy = _spectrum.sample(_randSpectrum->fire());
 
diff --git a/EventGenerator/src/MuCapNeutronGenerator_tool.cc b/EventGenerator/src/MuCapNeutronGenerator_tool.cc
index 7ee90cff12..d5dcdd3b19 100644
--- a/EventGenerator/src/MuCapNeutronGenerator_tool.cc
+++ b/EventGenerator/src/MuCapNeutronGenerator_tool.cc
@@ -1,4 +1,5 @@
 #include "art/Utilities/ToolMacros.h"
+#include "cetlib_except/exception.h"
 
 #include "CLHEP/Random/RandPoissonQ.h"
 #include "CLHEP/Random/RandGeneral.h"
@@ -25,6 +26,8 @@ namespace mu2e {
 
       fhicl::DelegatedParameter spectrum{Name("spectrum"), Comment("Parameters for BinnedSpectrum)")};
       fhicl::Atom<std::string> spectrumVariable{Name("spectrumVariable"), Comment("Variable that spectrum is defined in (\"totalEnergy\", \"kineticEnergy\", or \"momentum\")")};
+      fhicl::Atom<double> czMin{Name("czmin"), Comment("Restrict cos(theta_z) minimum"), -1.};
+      fhicl::Atom<double> czMax{Name("czmax"), Comment("Restrict cos(theta_z) maximum"),  1.};
     };
     typedef art::ToolConfigTable<PhysConfig> Parameters;
 
@@ -32,16 +35,22 @@ namespace mu2e {
       _pdgId(PDGCode::n0),
       _mass(GlobalConstantsHandle<ParticleDataList>()->particle(_pdgId).mass()),
       _spectrum(BinnedSpectrum(conf().spectrum.get<fhicl::ParameterSet>())),
-      _spectrumVariable(parseSpectrumVar(conf().spectrumVariable()))
-    {}
+      _spectrumVariable(parseSpectrumVar(conf().spectrumVariable())),
+      _czMin(conf().czMin()),
+      _czMax(conf().czMax())
+    {
+      if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax;
+    }
 
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _rate = GlobalConstantsHandle<PhysicsParams>()->getCaptureNeutronRate(material);
-      _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng);
-      _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, _rate);
+      const double rate = _rate * (_czMax - _czMin)/2.; // accouunt for potential cz selection 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());
     }
 
@@ -52,6 +61,8 @@ namespace mu2e {
 
     BinnedSpectrum    _spectrum;
     SpectrumVar       _spectrumVariable;
+    double _czMin;
+    double _czMax;
 
     std::unique_ptr<CLHEP::RandPoissonQ> _randomPoissonQ;
     std::unique_ptr<RandomUnitSphere>    _randomUnitSphere;
@@ -61,7 +72,7 @@ namespace mu2e {
   std::vector<ParticleGeneratorTool::Kinematic> MuCapNeutronGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  res;
 
-    int n_gen = _randomPoissonQ->fire();
+    const int n_gen = (_isPrimary) ? 1 : _randomPoissonQ->fire();
     for (int i_gen = 0; i_gen < n_gen; ++i_gen) {
       double energy = _spectrum.sample(_randSpectrum->fire());
 
diff --git a/EventGenerator/src/MuCapPhotonGenerator_tool.cc b/EventGenerator/src/MuCapPhotonGenerator_tool.cc
index 3897583668..e9e925bcd2 100644
--- a/EventGenerator/src/MuCapPhotonGenerator_tool.cc
+++ b/EventGenerator/src/MuCapPhotonGenerator_tool.cc
@@ -1,4 +1,5 @@
 #include "art/Utilities/ToolMacros.h"
+#include "cetlib_except/exception.h"
 
 #include "CLHEP/Random/RandPoissonQ.h"
 #include "CLHEP/Random/RandGeneral.h"
@@ -22,20 +23,28 @@ namespace mu2e {
       using Comment=fhicl::Comment;
 
       fhicl::DelegatedParameter spectrum{Name("spectrum"), Comment("Parameters for BinnedSpectrum)")};
+      fhicl::Atom<double> czMin{Name("czmin"), Comment("Restrict cos(theta_z) minimum"), -1.};
+      fhicl::Atom<double> czMax{Name("czmax"), Comment("Restrict cos(theta_z) maximum"),  1.};
     };
     typedef art::ToolConfigTable<PhysConfig> Parameters;
 
     explicit MuCapPhotonGenerator(Parameters const& conf) :
-      _spectrum(BinnedSpectrum(conf().spectrum.get<fhicl::ParameterSet>()))
-    {}
+      _spectrum(BinnedSpectrum(conf().spectrum.get<fhicl::ParameterSet>())),
+      _czMin(conf().czMin()),
+      _czMax(conf().czMax())
+    {
+      if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax;
+    }
 
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _rate = GlobalConstantsHandle<PhysicsParams>()->getCapturePhotonRate(material);
-      _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng);
-      _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, _rate);
+      const double rate = _rate * (_czMax - _czMin)/2.; // accouunt for potential cz selection 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());
     }
 
@@ -43,6 +52,8 @@ namespace mu2e {
     double _rate = 0.;
 
     BinnedSpectrum    _spectrum;
+    double _czMin;
+    double _czMax;
 
     std::unique_ptr<CLHEP::RandPoissonQ> _randomPoissonQ;
     std::unique_ptr<RandomUnitSphere>    _randomUnitSphere;
@@ -53,7 +64,7 @@ namespace mu2e {
   std::vector<ParticleGeneratorTool::Kinematic> MuCapPhotonGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  res;
 
-    int n_gen = _randomPoissonQ->fire();
+    const int n_gen = (_isPrimary) ? 1 : _randomPoissonQ->fire();
     for (int i_gen = 0; i_gen < n_gen; ++i_gen) {
       double energy = _spectrum.sample(_randSpectrum->fire());
       const double p = energy;
diff --git a/EventGenerator/src/MuCapProtonGenerator_tool.cc b/EventGenerator/src/MuCapProtonGenerator_tool.cc
index b9d0b4f74a..f0f1ccd6c3 100644
--- a/EventGenerator/src/MuCapProtonGenerator_tool.cc
+++ b/EventGenerator/src/MuCapProtonGenerator_tool.cc
@@ -1,4 +1,5 @@
 #include "art/Utilities/ToolMacros.h"
+#include "cetlib_except/exception.h"
 
 #include "CLHEP/Random/RandPoissonQ.h"
 #include "CLHEP/Random/RandGeneral.h"
@@ -25,6 +26,8 @@ namespace mu2e {
 
       fhicl::DelegatedParameter spectrum{Name("spectrum"), Comment("Parameters for BinnedSpectrum)")};
       fhicl::Atom<std::string> spectrumVariable{Name("spectrumVariable"), Comment("Variable that spectrum is defined in (\"totalEnergy\", \"kineticEnergy\", or \"momentum\")")};
+      fhicl::Atom<double> czMin{Name("czmin"), Comment("Restrict cos(theta_z) minimum"), -1.};
+      fhicl::Atom<double> czMax{Name("czmax"), Comment("Restrict cos(theta_z) maximum"),  1.};
     };
     typedef art::ToolConfigTable<PhysConfig> Parameters;
 
@@ -32,16 +35,22 @@ namespace mu2e {
       _pdgId(PDGCode::proton),
       _mass(GlobalConstantsHandle<ParticleDataList>()->particle(_pdgId).mass()),
       _spectrum(BinnedSpectrum(conf().spectrum.get<fhicl::ParameterSet>())),
-      _spectrumVariable(parseSpectrumVar(conf().spectrumVariable()))
-    {}
+      _spectrumVariable(parseSpectrumVar(conf().spectrumVariable())),
+      _czMin(conf().czMin()),
+      _czMax(conf().czMax())
+    {
+      if(_czMin < -1. || _czMax > 1. || _czMin > _czMax) throw cet::exception("BADCONFIG") << "Cos(theta_z) range unphysical: " << _czMin << " - " << _czMax;
+    }
 
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _rate = GlobalConstantsHandle<PhysicsParams>()->getCaptureProtonRate(material);
-      _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng);
-      _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, _rate);
+      const double rate = _rate * (_czMax - _czMin)/2.; // accouunt for potential cz selection 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());
     }
 
@@ -52,6 +61,8 @@ namespace mu2e {
 
     BinnedSpectrum    _spectrum;
     SpectrumVar       _spectrumVariable;
+    double _czMin;
+    double _czMax;
 
     std::unique_ptr<CLHEP::RandPoissonQ> _randomPoissonQ;
     std::unique_ptr<RandomUnitSphere>    _randomUnitSphere;
@@ -61,7 +72,7 @@ namespace mu2e {
   std::vector<ParticleGeneratorTool::Kinematic> MuCapProtonGenerator::generate() {
     std::vector<ParticleGeneratorTool::Kinematic>  res;
 
-    int n_gen = _randomPoissonQ->fire();
+    const int n_gen = (_isPrimary) ? 1 : _randomPoissonQ->fire();
     for (int i_gen = 0; i_gen < n_gen; ++i_gen) {
       double energy = _spectrum.sample(_randSpectrum->fire());
 
diff --git a/EventGenerator/src/MuStopProductsGun_module.cc b/EventGenerator/src/MuStopProductsGun_module.cc
index a5521010e7..fd8b3c7c4a 100644
--- a/EventGenerator/src/MuStopProductsGun_module.cc
+++ b/EventGenerator/src/MuStopProductsGun_module.cc
@@ -107,12 +107,12 @@ namespace mu2e {
     const auto cap_psets = conf_.captureProducts.get<std::vector<fhicl::ParameterSet>>();
     for (const auto& i_cap_pset : cap_psets) {
       _muonCaptureGenerators.push_back(art::make_tool<ParticleGeneratorTool>(i_cap_pset));
-      _muonCaptureGenerators.back()->finishInitialization(eng_, material_);
+      _muonCaptureGenerators.back()->finishInitialization(eng_, material_, false);
     }
     const auto decay_psets = conf_.decayProducts.get<std::vector<fhicl::ParameterSet>>();
     for (const auto& i_decay_pset : decay_psets) {
       _muonDecayGenerators.push_back(art::make_tool<ParticleGeneratorTool>(i_decay_pset));
-      _muonDecayGenerators.back()->finishInitialization(eng_, material_);
+      _muonDecayGenerators.back()->finishInitialization(eng_, material_, false);
     }
   }
 
diff --git a/EventGenerator/src/MuplusMichelGenerator_tool.cc b/EventGenerator/src/MuplusMichelGenerator_tool.cc
index 0e43e5b83a..9a1883ff7c 100644
--- a/EventGenerator/src/MuplusMichelGenerator_tool.cc
+++ b/EventGenerator/src/MuplusMichelGenerator_tool.cc
@@ -36,7 +36,8 @@ namespace mu2e {
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string&) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string&, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _randomUnitSphere = std::make_unique<RandomUnitSphere>(eng);
       _randSpectrum = std::make_unique<CLHEP::RandGeneral>(eng, _spectrum.getPDF(), _spectrum.getNbins());
     }
diff --git a/EventGenerator/src/Pileup_module.cc b/EventGenerator/src/Pileup_module.cc
index 558c8eb488..dddb333cb1 100644
--- a/EventGenerator/src/Pileup_module.cc
+++ b/EventGenerator/src/Pileup_module.cc
@@ -111,13 +111,13 @@ namespace mu2e {
     const auto cap_psets = conf().captureProducts.get<std::vector<fhicl::ParameterSet>>();
     for (const auto& pset : cap_psets) {
       muonCaptureGenerators_.push_back(art::make_tool<ParticleGeneratorTool>(pset));
-      muonCaptureGenerators_.back()->finishInitialization(eng_, conf().stoppingTargetMaterial());
+      muonCaptureGenerators_.back()->finishInitialization(eng_, conf().stoppingTargetMaterial(), false);
     }
 
     const auto decay_psets = conf().decayProducts.get<std::vector<fhicl::ParameterSet>>();
     for (const auto& pset : decay_psets) {
       muonDecayGenerators_.push_back(art::make_tool<ParticleGeneratorTool>(pset));
-      muonDecayGenerators_.back()->finishInitialization(eng_, conf().stoppingTargetMaterial());
+      muonDecayGenerators_.back()->finishInitialization(eng_, conf().stoppingTargetMaterial(), false);
     }
   }
 
diff --git a/EventGenerator/src/RMCGenerator_tool.cc b/EventGenerator/src/RMCGenerator_tool.cc
index 63e804f977..d4174bcf08 100644
--- a/EventGenerator/src/RMCGenerator_tool.cc
+++ b/EventGenerator/src/RMCGenerator_tool.cc
@@ -98,14 +98,15 @@ namespace mu2e {
     std::vector<ParticleGeneratorTool::Kinematic> generate() override;
     void generate(std::unique_ptr<GenParticleCollection>& out, const IO::StoppedParticleF& stop) override;
 
-    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material) override {
+    void finishInitialization(art::RandomNumberGenerator::base_engine_t& eng, const std::string& material, const bool isPrimary) override {
+      _isPrimary = isPrimary;
       _randomUnitSphereExternal = std::make_unique<RandomUnitSphere>(eng, _czmin, _czmax);
       _randomUnitSphereInternal = std::make_unique<RandomUnitSphere>(eng);
       _randFlat = std::make_unique<CLHEP::RandFlat>(eng);
       _randSpectrum = std::make_unique<CLHEP::RandGeneral>(eng, _spectrum.getPDF(), _spectrum.getNbins());
       _muonCaptureSpectrum = std::make_unique<MuonCaptureSpectrum>(_randFlat.get(), _randomUnitSphereInternal.get());
       if(_useRate) {
-        _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, GlobalConstantsHandle<PhysicsParams>()->getCaptureRMCRate(material));
+        _randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, GlobalConstantsHandle<PhysicsParams>()->getCaptureRMCRate(material) * (_czmax - _czmin)/2.);
         _internalRate = GlobalConstantsHandle<PhysicsParams>()->getCaptureRMCInternalRate(material);
       }
     }
diff --git a/EventGenerator/src/SingleProcessGenerator_module.cc b/EventGenerator/src/SingleProcessGenerator_module.cc
index b30141d5bf..1192d245fb 100644
--- a/EventGenerator/src/SingleProcessGenerator_module.cc
+++ b/EventGenerator/src/SingleProcessGenerator_module.cc
@@ -99,7 +99,7 @@ namespace mu2e {
 
     const auto pset = conf().decayProducts.get<fhicl::ParameterSet>();
     Generator_ = art::make_tool<ParticleGeneratorTool>(pset);
-    Generator_->finishInitialization(eng_, conf().stoppingTargetMaterial());
+    Generator_->finishInitialization(eng_, conf().stoppingTargetMaterial(), true);
 
     if(pdgId_==PDGCode::e_plus) {
       muonLifeTime_=0; //decay time already included for stopped muon(+) FIXME!!!
diff --git a/Mu2eG4/src/CompositeMaterialGenerator_module.cc b/Mu2eG4/src/CompositeMaterialGenerator_module.cc
index 352f1db335..773f914a8b 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<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, "");
+      generator_tool->finishInitialization(_engine, "", false);
 
       // enforce that keys to subgenerators are unique
       // would be better in a protected class, but eh