add mean gamma energy for deposited gammas

gammapkt.cc

@@ -1112,6 +1112,7 @@ __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

globals.h

@@ -212,10 +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_nu_cmf_decayspec;
-inline double estimator_gamma_kappa_integrated;
-inline double estimator_gamma_kappa_decayspec;
-inline double estimator_gamma_kappa_absorbedspec;
+inline double estimator_gamma_kappa_integrated = 0.;
+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};
 

sn3d.cc

@@ -618,6 +618,7 @@ void zero_estimators() {
   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_integrated = 0.;
   globals::estimator_gamma_kappa_decayspec = 0.;
   globals::estimator_gamma_kappa_absorbedspec = 0.;
@@ -653,16 +654,21 @@ void normalise_deposition_estimators(int nts) {
     gamma_dep_erg += (globals::dep_estimator_gamma[nonemptymgi] / nprocs);
   }
   globals::estimator_gamma_kappa_integrated /= gamma_dep_erg;
+
   globals::estimator_gamma_kappa_decayspec /= globals::timesteps[nts].gamma_emission;
-  globals::estimator_gamma_kappa_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete;
   globals::estimator_gamma_nu_cmf_decayspec /= globals::timesteps[nts].gamma_emission;
-  const double mean_en_mev = H * globals::estimator_gamma_nu_cmf_decayspec / MEV;
+  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 %g days kappa_eff [cm2/g]: path-integrated %g decay_spec %g deposit_spec %g decay gamma mean E[MeV] "
-      "%g\n",
+      "timestep %d %.2f days kappa_eff [cm2/g]: path-integrated %.3f 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_integrated,
-      globals::estimator_gamma_kappa_decayspec, globals::estimator_gamma_kappa_absorbedspec, mean_en_mev);
+      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);