diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h
index 3f35ddf43..9de58357f 100644
--- a/artisoptions_kilonova_lte.h
+++ b/artisoptions_kilonova_lte.h
@@ -6,7 +6,7 @@
 
 #include "constants.h"
 
-constexpr int MPKTS = 15000;
+constexpr int MPKTS = 100000;
 
 constexpr auto GRID_TYPE = GridType::CARTESIAN3D;
 constexpr int CUBOID_NCOORDGRID_X = 50;
@@ -151,5 +151,9 @@ constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
 constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
+constexpr float BETAMINUS_ENERGY_GAMMA_SPLITUP = 0.45;
+
+constexpr float BETAMINUS_ENERGY_ELECTRON_SPLITUP = 0.2;
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H
diff --git a/decay.cc b/decay.cc
index 6326a684e..8761dd0de 100644
--- a/decay.cc
+++ b/decay.cc
@@ -999,6 +999,12 @@ void init_nuclides(const std::vector<int> &custom_zlist, const std::vector<int>
   // Read in data for gamma ray lines and make a list of them in energy order.
   gammapkt::init_gamma_data();
 
+  // manipulate betaminus decay splitup ratios here
+  for (int i = 0; i < static_cast<int>(nuclides.size()); i++) {
+    nuclides[i].endecay_gamma = BETAMINUS_ENERGY_GAMMA_SPLITUP * nuclides[i].endecay_q[DECAYTYPE_BETAMINUS];
+    nuclides[i].endecay_electron = BETAMINUS_ENERGY_ELECTRON_SPLITUP * nuclides[i].endecay_q[DECAYTYPE_BETAMINUS];
+  }
+
   // TODO: generalise this to all included nuclides
   printout("decayenergy(NI56), decayenergy(CO56), decayenergy_gamma(CO56): %g, %g, %g\n",
            nucdecayenergytotal(28, 56) / MEV, nucdecayenergytotal(27, 56) / MEV, nucdecayenergygamma(27, 56) / MEV);
diff --git a/gammapkt.cc b/gammapkt.cc
index 56893e361..65b7a3e8a 100644
--- a/gammapkt.cc
+++ b/gammapkt.cc
@@ -236,7 +236,10 @@ void init_xcom_photoion_data() {
 __host__ __device__ auto choose_gamma_ray(const int nucindex) -> double {
   // Routine to choose which gamma ray line it'll be.
 
-  const double E_gamma = decay::nucdecayenergygamma(nucindex);  // Average energy per gamma line of a decay
+  double E_gamma = 0.;  // Average energy per gamma line of a decay
+  for (ptrdiff_t n = 0; n < std::ssize(gamma_spectra[nucindex]); n++) {
+    E_gamma += gamma_spectra[nucindex][n].probability * gamma_spectra[nucindex][n].energy;
+  }
 
   const double zrand = rng_uniform();
   double runtot = 0.;