diff --git a/gammapkt.cc b/gammapkt.cc
index 720815db4..278bc331d 100644
--- a/gammapkt.cc
+++ b/gammapkt.cc
@@ -647,6 +647,16 @@ constexpr auto meanf_sigma(const double x) -> double {
   return tot;
 }
 
+// get the gamma-ray opacity (with the expected energy loss per interaction factor included)
+auto get_kappa(const Packet &pkt) -> double {
+  const double xx = H * pkt.nu_cmf / ME / CLIGHT / CLIGHT;
+  const int mgi = grid::get_cell_modelgridindex(pkt.where);
+  const double chi = ((meanf_sigma(xx) * grid::get_nnetot(mgi)) + get_chi_photo_electric_rf(pkt) +
+                      (sigma_pair_prod_rf(pkt) * (1. - (2.46636e+20 / pkt.nu_cmf))));
+  const double kappa = 1 / grid::get_rho(mgi) * chi;
+  return kappa;
+}
+
 // Subroutine to record the heating rate in a cell due to gamma rays.
 // By heating rate I mean, for now, really the rate at which the code is making
 // k-packets in that cell which will then convert into r-packets. This is (going
@@ -1079,6 +1089,9 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) {
   // printout("initialise pol state of packet %g, %g, %g, %g,
   // %g\n",pkt.stokes_qu[0],pkt.stokes_qu[1],pkt.pol_dir[0],pkt.pol_dir[1],pkt.pol_dir[2]);
   // printout("pkt direction %g, %g, %g\n",pkt.dir[0],pkt.dir[1],pkt.dir[2]);
+
+  atomicadd(globals::estimator_gamma_kappa_decayspec, get_kappa(pkt) * pkt.e_cmf);
+  atomicadd(globals::estimator_gamma_nu_cmf_decayspec, pkt.nu_cmf * pkt.e_cmf);
 }
 
 __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {
@@ -1096,6 +1109,8 @@ __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {
 
   if (pkt.type != TYPE_GAMMA && pkt.type != TYPE_ESCAPE) {
     atomicadd(globals::timesteps[nts].gamma_dep_discrete, pkt.e_cmf);
+    atomicadd(globals::estimator_gamma_kappa_absorbedspec, get_kappa(pkt) * pkt.e_cmf);
+    atomicadd(globals::estimator_gamma_nu_cmf_absorbedspec, pkt.nu_cmf * pkt.e_cmf);
 
     if constexpr (GAMMA_THERMALISATION_SCHEME != ThermalisationScheme::DETAILED) {
       // no transport, so the path-based gamma deposition estimator won't get updated unless we do it here
diff --git a/globals.h b/globals.h
index d46106e78..dd52d73ce 100644
--- a/globals.h
+++ b/globals.h
@@ -212,6 +212,11 @@ inline std::vector<double> dep_estimator_positron;
 inline std::vector<double> dep_estimator_electron;
 inline std::vector<double> dep_estimator_alpha;
 
+inline double estimator_gamma_kappa_decayspec = 0.;
+inline double estimator_gamma_nu_cmf_decayspec = 0.;
+inline double estimator_gamma_kappa_absorbedspec = 0.;
+inline double estimator_gamma_nu_cmf_absorbedspec = 0.;
+
 inline int bfestimcount{0};
 
 // for USE_LUT_PHOTOION = true
diff --git a/sn3d.cc b/sn3d.cc
index 568fe2733..6d5ecf4b0 100644
--- a/sn3d.cc
+++ b/sn3d.cc
@@ -617,6 +617,11 @@ void zero_estimators() {
   std::ranges::fill(globals::dep_estimator_electron, 0.);
   std::ranges::fill(globals::dep_estimator_alpha, 0.);
 
+  globals::estimator_gamma_nu_cmf_decayspec = 0.;
+  globals::estimator_gamma_nu_cmf_absorbedspec = 0.;
+  globals::estimator_gamma_kappa_decayspec = 0.;
+  globals::estimator_gamma_kappa_absorbedspec = 0.;
+
   if constexpr (USE_LUT_PHOTOION) {
     if (globals::nbfcontinua_ground > 0) {
       std::ranges::fill(globals::gammaestimator, 0.);
@@ -643,6 +648,20 @@ void normalise_deposition_estimators(int nts) {
   globals::timesteps[nts].electron_dep = 0.;
   globals::timesteps[nts].alpha_dep = 0.;
 
+  globals::estimator_gamma_kappa_decayspec /= globals::timesteps[nts].gamma_emission;
+  globals::estimator_gamma_nu_cmf_decayspec /= globals::timesteps[nts].gamma_emission;
+  const double mean_en_decay_mev = H * globals::estimator_gamma_nu_cmf_decayspec / MEV;
+
+  globals::estimator_gamma_kappa_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete;
+  globals::estimator_gamma_nu_cmf_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete;
+  const double mean_en_absorb_mev = H * globals::estimator_gamma_nu_cmf_absorbedspec / MEV;
+
+  printout(
+      "timestep %d %.2f days kappa_eff [cm2/g]: decay_spec %.3f deposit_spec %.3f. mean gamma "
+      "E[MeV]: decay %.2f deposit %.2f\n",
+      nts, globals::timesteps[nts].mid / DAY, globals::estimator_gamma_kappa_decayspec,
+      globals::estimator_gamma_kappa_absorbedspec, mean_en_decay_mev, mean_en_absorb_mev);
+
   for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) {
     const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);