/*******************************************************************************
* Component: Masked_ESS_butterfly
*
* An ESS_butterfly source that emits only where a chopper train lets neutrons through.
*
* `choppers` is a `chopper_parameters` array reinterpreted as `double *`, which is the
* only pointer a McStas SETTING PARAMETER can carry. It was four doubles per chopper
* before chopper-lib 4.0.0 and is a structure holding a pointer now, so an instrument
* that really did hand over a flat array of numbers must build the structures instead:
*
*   double edges[] = {-1.8, 1.8};
*   chopper_parameters pars[] = {{14.0, 0.02, 0.0, 2, edges, 10.0, 4.0}};
*   chopper_ptr = (double *) pars;
*   chopper_cnt = sizeof(pars) / sizeof(chopper_parameters);
*
* `edges` is a flat, increasing list of angles in degrees, two per opening, in the disk's
* own frame; `beam` is the angle on the beam path at `delay`; and `aperture` is how wide
* the beam is on the disk, in degrees about its spindle, which widens every window in time
* and is what a mask for anything but a pencil beam needs. Leave it off and the row still
* compiles, describing the point beam it described before. See the chopper-lib README.
*
* <b>What the mask is spent on: `resample`</b>
*
* By default a ray the mask excludes is ABSORBed. That is correct -- the discs would have
* stopped it -- but a train that passes a percent of the plane then spends ninety-nine
* percent of `ncount` on rays that die where they were born, and the run has the statistics
* of one a hundredth the size.
*
* `resample = 1` draws the excluded ray again, from inside the allowed region, instead of
* throwing it away. This is the same measurement, not an approximation, because the source
* samples the two coordinates the mask is drawn in uniformly and independently:
* `lambda = Lmin + range * rand01()` with `1/v = lambda V2K / 2 pi`, and
* `t = rand01() * tmax_multiplier * ESS_SOURCE_DURATION`. Restricting a uniform draw to a
* region and multiplying the weight by that region's share of the whole -- chopper-lib's
* `acceptance` -- leaves every downstream estimator unbiased. Every emitted ray is
* multiplied by it, the ones that landed in the allowed region unaided included.
*
* Two things this factor is not, both easy to reach for: it is not
* `chopper_unmasked_probability`, which is the allowed fraction of the *weighted* signal
* and so the transmission rather than the normalisation; and it is not recoverable from
* counting rejection attempts, since `E[1/k]` is not `1/E[k]`.
*
* The independence the argument rests on fails under time focusing, and INITIALIZE refuses
* `resample` there rather than quietly biasing the answer: the emission time is then drawn
* in a window centred on `tfocus_time - tfocus_dist / vz`, which moves with the neutron's
* own velocity. No single factor corrects that -- but none is needed, because time focusing
* already *is* this trick done exactly. The window traces a band of slope `-tfocus_dist` in
* the (1/v, t) plane, which is the shape of a chopper band, and `w_tfocus` is already the
* compensating weight. Set the three `tfocus_*` parameters from the pulse-shaping disc and
* the source restricts itself.
*
* <b>Under MPI</b>
*
* The three grids this component keeps are summed across the nodes in SAVE and written by
* master alone -- every node runs SAVE, and they all share one output directory, so
* ungated writes would put one interleaved copy per node into each file. The acceptance is
* not summed and must not be: it comes from the mask's geometry, so every node has the same
* number already. Nor is anything divided by the node count, because the ray weights carry
* one over the *whole* run's ncount: ESS_butterfly reads mcget_ncount() in INITIALIZE,
* which mccode_main runs before it slices ncount across the nodes.
*
* MPI is the parallelism this has been run under. The grids are accumulated atomically, so
* OpenACC threads are safe too.
*
* Resampling changes only how the surviving rays are drawn, never which region survives.
* It inherits the mask's own approximations -- bin quantisation, `mask_grow`, the disc
* `aperture` -- exactly, and adds none.
*
* <b>What the grid has to span</b>
*
* Neither of the two coordinates is read straight off the ray, and getting either wrong
* absorbs the beam rather than masking it.
*
* Under time focusing the emission time is drawn in a `tfocus_width` window whose centre
* slides with the neutron's inverse velocity, so the band it sweeps is
* `tfocus_dist * inverse_velocity_range` longer than the window; a grid sized by
* `tfocus_width` alone sits around the slowest neutron's window and nowhere near the rest.
* And ESS_butterfly adds the offset of a randomly chosen pulse to `t` before this component
* sees it, so with `n_pulses > 1` what arrives is not an emission time at all. The grid is
* built over the whole swept band clipped to the pulse, and the pulse offset is taken back
* off in TRACE -- which needs the emission window to be shorter than the gap between
* pulses, so INITIALIZE checks that too.
*
* %P
* INPUT PARAMETERS:
*
* choppers: [chopper_parameters,...]  information about the chopper train (see above)
* chopper_count: [1]                  the number of choppers described
* inverse_velocity_bin: [s/m]         the width of mask bins in 1/v
* time_bin : [s]                      the width of emission time bins
* filename: [str]                     the base name of any output files, replaced by NAME_CURRENT_COMP if missing
* noise_fraction: [1.]                accept out-of-mask rays with this probability
* mask_grow: [1]                      expand the chopper-accepted mask by this many bins in each direction
* use_mask: [1]                       use (1) or do not use (0) the calculated mask for limiting source emission
* resample: [1]                       draw (1/v, t) directly from the mask, 100% ray emission from source
* save_mask: [1]                      output the binary mask as {filename}.mask
* save_total: [1]                     output the probability drawn from the source before masking as {filename}.total
* save_emitted: [1]                   output the probability emitted after masking as {filename}.emitted
* save_count: [1]                     output the number of drawn rays before masking as {filename}.count
* verify_count: [1]                   report the fraction of drawn rays landing inside the mask, which measures
*                                     the acceptance that resample=1 corrects every ray weight by, and should
*                                     match it -- it is counted from the run, where the acceptance comes from
*                                     the mask geometry alone, so the two agreeing is a real check

*******************************************************************************/
DEFINE COMPONENT Masked_ESS_butterfly
INHERIT ESS_butterfly
DEFINITION PARAMETERS ()
SETTING PARAMETERS (
  double * choppers,
  int chopper_count,
  double inverse_velocity_bin,
  double time_bin,
  string filename=0,
  noise_fraction=0.0,
  int mask_grow=1,
  int use_mask=1,
  int resample=0,
  int save_mask=1,
  int save_total=1,
  int save_emitted=1,
  int save_count=1,
  int verify_count=1
)
OUTPUT PARAMETERS ()
/* The three grids below are written from TRACE with `#pragma acc atomic`, so the accumulation
 * is race-free as it stands, and everything they and the sampler point at is ordinary
 * calloc'd memory that `-gpu=mem:managed` -- the flag every McCode OpenACC toolchain sets --
 * makes reachable from a kernel. That covers the bytes; it does not cover the *pointer
 * fields* themselves in this component's device-resident twin, which McStas's own generated
 * `update device(...)` only shallow-copies -- so INITIALIZE attaches each one explicitly with
 * `acc_attach`, the same fix `off_init` elsewhere in the generated file uses for its own
 * malloc'd mesh arrays. */
SHARE INHERIT ESS_butterfly EXTEND %{
  %include "chopper-lib"

#if !defined(CHOPPER_LIB_VERSION) || CHOPPER_LIB_VERSION < 40100
#error "Masked_ESS_butterfly describes choppers by edges and resamples with chopper_mask_sampler; chopper-lib 4.1.0 or newer is required"
#endif
%}
DECLARE INHERIT ESS_butterfly EXTEND %{
  double * total;
  double * emitted;
  double * counted;
  int * mask;
  double * inverse_velocity_edges;
  double * time_edges;
  unsigned inverse_velocity_count;
  unsigned time_count;
  double minimum_inverse_velocity_edge;
  double inverse_velocity_range;
  double minimum_time_edge;
  double time_range;
  double pulse_period;
  chopper_mask_sampler mask_sampler;
  int save_or_verify_count;
%}
INITIALIZE INHERIT ESS_butterfly EXTEND %{
  // Use the same minimum inverse velocity as in ESS_butterfly, which draws its wavelength
  // uniformly and so its inverse velocity uniformly too: 1/v = lambda V2K / 2 pi exactly.
  minimum_inverse_velocity_edge = Lmin * V2K / 2 / PI;
  inverse_velocity_range = (Lmax - Lmin) * V2K / 2 / PI;
  if (inverse_velocity_bin <= 0) inverse_velocity_bin = inverse_velocity_range;
  inverse_velocity_count = (unsigned) ceil(inverse_velocity_range / inverse_velocity_bin);
  inverse_velocity_edges = (double *) calloc(inverse_velocity_count + 1, sizeof(double));

  /* The emission times the source can produce, which is what the grid has to span.
   *
   * Without time focusing ESS_butterfly draws t uniformly on [0, tmax_multiplier *
   * ESS_SOURCE_DURATION). With it, t is drawn in a window tfocus_width wide centred on
   * `tfocus_time - tfocus_dist / vz` -- a centre that slides with the neutron's own inverse
   * velocity, so the band that window sweeps out is tfocus_dist * inverse_velocity_range
   * longer than the window itself. Sizing the grid by tfocus_width alone puts it around the
   * slowest neutron's window and nowhere near the rest, and the bounds check in TRACE then
   * absorbs the entire beam as though the choppers had stopped it.
   *
   * vz rather than v, strictly, but the mask is drawn against the along-the-path inverse
   * velocity throughout and the two differ by the direction cosine, which is 1 to a part in
   * 1e4 for any sensible focusing rectangle.
   *
   * Either way the source absorbs whatever falls outside the pulse, so the reachable
   * interval is that band clipped to it.
   */
  double emission_end = tmax_multiplier * ESS_SOURCE_DURATION;
  if (tfocus_width > 0) {
    double earliest = tfocus_time - tfocus_dist * (minimum_inverse_velocity_edge + inverse_velocity_range) - tfocus_width / 2.0;
    double latest = tfocus_time - tfocus_dist * minimum_inverse_velocity_edge + tfocus_width / 2.0;
    if (earliest < 0.0) earliest = 0.0;
    if (latest > emission_end) latest = emission_end;
    minimum_time_edge = earliest;
    time_range = latest - earliest;
    if (time_range <= 0.0) {
      MPI_MASTER(
        fprintf(stderr, "%s: the time focusing window reaches [%g, %g] s, which does not "
                        "overlap the [0, %g] s pulse the source emits in, so nothing is "
                        "emitted at all\n",
                NAME_CURRENT_COMP,
                tfocus_time - tfocus_dist * (minimum_inverse_velocity_edge + inverse_velocity_range) - tfocus_width / 2.0,
                tfocus_time - tfocus_dist * minimum_inverse_velocity_edge + tfocus_width / 2.0,
                emission_end);
      );
      exit(1);
    }
  } else {
    minimum_time_edge = 0.0;
    time_range = emission_end;
  }
  if (time_bin <= 0) time_bin = time_range;
  time_count = (unsigned) ceil(time_range / time_bin);
  time_edges = (double *) calloc(time_count + 1, sizeof(double));

  /* ESS_butterfly picks a pulse per ray and adds its offset to t, so what arrives in TRACE
   * is not the emission time the mask is drawn against. The offset comes back off there --
   * unambiguously, so long as the window the source emits in is shorter than the gap
   * between pulses, which it is for any ordinary tmax_multiplier. */
  pulse_period = 1.0 / ESS_SOURCE_FREQUENCY;
  if (n_pulses > 1 && (minimum_time_edge < 0.0 || minimum_time_edge + time_range > pulse_period)) {
    MPI_MASTER(
      fprintf(stderr, "%s: with n_pulses=%d the emission window [%g, %g] s has to fit inside "
                      "one %g s pulse period, or the offset ESS_butterfly adds to t cannot be "
                      "told from the emission time it is added to\n",
              NAME_CURRENT_COMP, n_pulses, minimum_time_edge,
              minimum_time_edge + time_range, pulse_period);
    );
    exit(1);
  }

  save_or_verify_count = verify_count ? 1 : save_count ? 1 : 0;
  total = save_total ? (double *) calloc(inverse_velocity_count * time_count, sizeof(double)) : NULL;
  emitted = save_emitted ? (double *) calloc(inverse_velocity_count * time_count, sizeof(double)) : NULL;
  counted = save_or_verify_count ? (double *) calloc(inverse_velocity_count * time_count, sizeof(double)) : NULL;
  mask = (int *) calloc(inverse_velocity_count * time_count, sizeof(int));

  if (!inverse_velocity_edges || !time_edges || !mask
      || (save_total && !total) || (save_emitted && !emitted)
      || (save_or_verify_count && !counted)){
    /* Unlike every other failure here this one is a single node's own news, so it is printed
     * unreduced -- the node that ran out may not be master -- and it aborts the job rather
     * than calling exit, which under MPI means MPI_Finalize: finalizing alone while the
     * other nodes are still working leaves them waiting in the next collective forever. */
    fprintf(stderr, "Out of memory in %s!\n", NAME_CURRENT_COMP);
#ifdef USE_MPI
    MPI_Abort(MPI_COMM_WORLD, -1);
#endif
    exit(1);
  }
  inverse_velocity_edges[0] = minimum_inverse_velocity_edge;
  time_edges[0] = minimum_time_edge;
  for (unsigned it=0; it<time_count; ++it){
    time_edges[it+1] = time_edges[it] + time_bin;
  }
  for (unsigned iv=0; iv<inverse_velocity_count; ++iv){
    inverse_velocity_edges[iv+1] = inverse_velocity_edges[iv] + inverse_velocity_bin;
  }
  // The mask function overwrites every bin, but `total`, `emitted` and `counted` are
  // accumulated into, so they start empty. Index them the way everything else does:
  // `it * inverse_velocity_count + iv`, which is what chopper_inverse_velocity_time_mask
  // fills and what the TRACE below reads.
  for (unsigned it=0; it<time_count; ++it){
    for (unsigned iv=0; iv<inverse_velocity_count; ++iv){
      const unsigned index = it * inverse_velocity_count + iv;
      mask[index] = CHOPPER_MASK_INCLUDED;
      if (save_total) total[index] = 0;
      if (save_emitted) emitted[index] = 0;
      if (save_or_verify_count) counted[index] = 0;
    }
  }

  chopper_parameters * chop_pars = (chopper_parameters *) choppers;

  /* A caller still filling the pre-4.0.0 layout positionally compiles clean and lands a
   * flight path in `edge_count` and rubbish in `edges`, so check before dereferencing it:
   * a disk described by fewer than two edges is not a disk. */
  for (int ci = 0; ci < chopper_count; ++ci) {
    if (chop_pars[ci].edge_count < 2 || chop_pars[ci].edge_count % 2
        || chop_pars[ci].edges == NULL) {
      MPI_MASTER(
        fprintf(stderr, "%s: chopper %d has %u edges at %p; each opening needs two, and "
                        "chopper-lib 4.0.0 takes {speed, delay, beam, edge_count, edges, "
                        "path}\n",
                NAME_CURRENT_COMP, ci, chop_pars[ci].edge_count, (void *) chop_pars[ci].edges);
      );
      exit(1);
    }
  }

  int allowed_bin_count = chopper_inverse_velocity_time_mask(
    mask, inverse_velocity_count, time_count,
    inverse_velocity_edges, inverse_velocity_count + 1,
    time_edges, time_count + 1,
    chop_pars, chopper_count,
    mask_grow /* grow the region by this many bins in each direction */
  );

  if (allowed_bin_count < 1){
    MPI_MASTER(
      fprintf(stderr, "Choppers allow no transmission from %s!", NAME_CURRENT_COMP);
    );
    exit(1);
  }

  /* mask/total/emitted/counted are read or written from TRACE, which runs on the GPU; attach
   * their pointer fields so the device-resident twin of this component points at them, the
   * same fix `off_init` uses for its own malloc'd mesh arrays elsewhere in this file. */
  #ifdef OPENACC
  acc_attach((void *)&mask);
  if (save_total)          acc_attach((void *)&total);
  if (save_emitted)        acc_attach((void *)&emitted);
  if (save_or_verify_count) acc_attach((void *)&counted);
  #endif

  /* Drawing from the allowed region instead of absorbing outside it is only the same
   * measurement while the two mask coordinates are sampled uniformly and independently of
   * everything else. Refuse the two configurations where they are not, rather than apply a
   * correction that does not correct. */
  chopper_mask_sampler_empty(&mask_sampler);
  if (resample) {
    if (!use_mask) {
      MPI_MASTER(
        fprintf(stderr, "%s: resample=1 needs use_mask=1; there is nothing to resample "
                        "away from with the mask switched off\n", NAME_CURRENT_COMP);
      );
      exit(1);
    }
    if (noise_fraction != 0.0) {
      MPI_MASTER(
        fprintf(stderr, "%s: resample=1 and noise_fraction=%g are contradictory. Resampling "
                        "emits no masked-out ray at all, so there is no leak to set a rate "
                        "for\n", NAME_CURRENT_COMP, noise_fraction);
      );
      exit(1);
    }
    if (tfocus_width > 0) {
      MPI_MASTER(
        fprintf(stderr, "%s: resample=1 cannot be used with time focusing. The emission time "
                        "is drawn in a window centred on tfocus_time - tfocus_dist/vz, which "
                        "moves with the neutron's own velocity, so no single weight factor "
                        "puts the restriction back.\n"
                        "  Time focusing already does this exactly: the window is a band of "
                        "slope -tfocus_dist in the (1/v, t) plane, the shape of a chopper "
                        "band, and w_tfocus is the compensating weight. Set tfocus_dist, "
                        "tfocus_time and tfocus_width from the pulse-shaping disc and leave "
                        "resample=0.\n", NAME_CURRENT_COMP);
      );
      exit(1);
    }
    mask_sampler = chopper_mask_sampler_make(
      mask, inverse_velocity_count, time_count,
      inverse_velocity_edges, time_edges,
      minimum_inverse_velocity_edge, inverse_velocity_range,
      minimum_time_edge, time_range
    );
    if (mask_sampler.count == 0) {
      MPI_MASTER(
        fprintf(stderr, "%s: the mask allows %d bins but none of them overlap the region the "
                        "source samples; nothing can be resampled into\n",
                NAME_CURRENT_COMP, allowed_bin_count);
      );
      exit(1);
    }
    // Same reason as mask/total/etc above: chopper_mask_sampler_draw is `#pragma acc routine
    // seq` and runs on the GPU, so mask_sampler's five calloc'd arrays need their pointer
    // fields attached too.
    #ifdef OPENACC
    acc_attach((void *)&mask_sampler.cumulative);
    acc_attach((void *)&mask_sampler.inverse_velocity_low);
    acc_attach((void *)&mask_sampler.inverse_velocity_width);
    acc_attach((void *)&mask_sampler.time_low);
    acc_attach((void *)&mask_sampler.time_width);
    #endif
    MPI_MASTER(
      printf("%s: resampling into %u of %u mask cells; acceptance %.6g, so the ray weight "
             "carries that factor and the run keeps its whole ncount\n",
             NAME_CURRENT_COMP, mask_sampler.count,
             inverse_velocity_count * time_count, mask_sampler.acceptance);
    );
  }

  if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP);
%}
TRACE INHERIT ESS_butterfly EXTEND %{
  // Since this is after the trace of ESS_butterfly, a neutron ray has already been selected.

  // ESS_butterfly has already chosen which pulse this ray belongs to and added that pulse's
  // offset to t. The mask is drawn against the emission time inside a single pulse, so take
  // the offset back off before looking anything up, and put it back on anything redrawn.
  // INITIALIZE has checked the emission window is shorter than the gap between pulses,
  // which is what makes the split unambiguous.
  double pulse_offset = n_pulses > 1 ? pulse_period * floor(t / pulse_period) : 0.0;
  double t_emission = t - pulse_offset;

  // Decide which bin the generated neutron ray should go into:
  double inv_v = 1.0 / sqrt(vx*vx + vy*vy + vz*vz);
  unsigned inverse_velocity_index = (unsigned) floor((inv_v - minimum_inverse_velocity_edge) / inverse_velocity_bin);
  unsigned time_index = (unsigned) floor((t_emission - minimum_time_edge) / time_bin);
  // Bounds check to prevent out-of-bounds array access
  if (inverse_velocity_index >= inverse_velocity_count || time_index >= time_count) {
    ABSORB;
  }
  // this indexing must match the internal working of 'chopper_inverse_velocity_time_mask'!
  unsigned linear_index = time_index * inverse_velocity_count + inverse_velocity_index;

  // Record the total probability distribution sampled, and the sampling itself. `counted`
  // carries no weight on purpose: the acceptance the resampling corrects by counts draws,
  // not intensity, so this is what SAVE checks it against.
  //
  // Atomic because every ray in flight shares these grids. Written the way McCode's own
  // monitors write theirs -- see PSD_monitor.comp -- as a read and a store rather than a
  // compound assignment.
  if (save_total) {
    #pragma acc atomic
    total[linear_index] = total[linear_index] + p;
  }
  if (save_or_verify_count) {
    #pragma acc atomic
    counted[linear_index] = counted[linear_index] + 1.0;
  }

  if (use_mask) {
    if (resample) {
      if (mask[linear_index] == CHOPPER_MASK_EXCLUDED) {
        // Draw a fresh (inverse velocity, emission time) from inside the allowed region and
        // rebuild everything downstream of the wavelength. The emission point, the surface
        // it came from and the direction it was focused into are all still good -- they are
        // sampled independently of these two coordinates -- so ESS_butterfly's own TRACE
        // locals are reused rather than redrawn, and only the wavelength-dependent tail of
        // its weight is repeated here.
        chopper_mask_sampler_draw(&mask_sampler, rand01(), rand01(), rand01(), &inv_v, &t_emission);
        lambda = inv_v * 2 * PI / V2K;  // the source's own 1/v = lambda V2K / 2 pi, inverted
        k = 2 * PI / lambda;
        v = K2V * k;
        vz = v * dz / r;
        vy = v * dy / r;
        vx = v * dx / r;

        // ESS_butterfly.comp:588-607, with dt = 0 and w_tfocus = 1 because time focusing is
        // refused above. The brilliance is a function of the emission time, which is why the
        // parent evaluates it before adding the pulse offset and why this does too. The
        // Schoenfeldt functions assign the weight rather than scaling it, so p needs no
        // initialisation.
        if (iscold) {
          ESS_2015_Schoenfeldt_cold(&t_emission, &p, lambda, tfocus_width, tfocus_time, dt, yheight,
                                    Mwidth_t, yheight, Mwidth_c, tmax_multiplier,
                                    beamportangle, modX, modY);
          p *= c_performance;
          p *= ColdScalars[beamline - 1];
        } else {
          ESS_2015_Schoenfeldt_thermal(&t_emission, &p, lambda, tfocus_width, tfocus_time, dt, yheight,
                                       Mwidth_t, yheight, Mwidth_c, tmax_multiplier,
                                       beamportangle, modX, modY);
          p *= t_performance;
          p *= ThermalScalars[beamline - 1];
        }
        p *= w_stat * w_focus * w_geom * w_mult * w_tfocus;
        p *= cos_factor;
        if (iscold) {
          p /= cold_frac;
        } else {
          p /= (1 - cold_frac);
        }
        // The pulse this ray was assigned to is untouched by the redraw: it is chosen
        // independently of the wavelength and the emission time, so conditioning those two
        // on the mask leaves it exactly as ESS_butterfly drew it.
        t = t_emission + pulse_offset;

        // Re-bin the redrawn ray from its velocity, the same way the original was binned.
        inv_v = 1.0 / sqrt(vx*vx + vy*vy + vz*vz);
        inverse_velocity_index = (unsigned) floor((inv_v - minimum_inverse_velocity_edge) / inverse_velocity_bin);
        time_index = (unsigned) floor((t_emission - minimum_time_edge) / time_bin);
        if (inverse_velocity_index >= inverse_velocity_count || time_index >= time_count) {
          ABSORB;  // the sampler draws inside the grid, so this is unreachable short of a bug
        }
        linear_index = time_index * inverse_velocity_count + inverse_velocity_index;
      }
      // Every emitted ray pays the acceptance, whether it needed redrawing or not: the ones
      // that landed inside the allowed region on their own are a draw from the same
      // restricted distribution, and carry the same correction.
      p *= mask_sampler.acceptance;
    } else if (mask[linear_index] == CHOPPER_MASK_EXCLUDED && rand01() > noise_fraction) {
      // Masked-out, and the noise leak did not save it
      ABSORB;
    }
  }
  // Otherwise, let it continue on its way -- the SCATTER call was made in ESS_butterfly
  if (save_emitted) {
    #pragma acc atomic
    emitted[linear_index] = emitted[linear_index] + p;
  }
%}
SAVE %{
  /* SAVE runs on every MPI node -- mccode_main calls finally() on all of them and finally()
   * calls save() unconditionally -- so the reduction below is reached everywhere, which it
   * has to be: mc_MPI_Sum is a collective MPI_Allreduce and a master-only call deadlocks.
   * Only the file writing is master's alone. It has to be: mcuse_dir hands every node the
   * same output directory, so without the gate every node writes the same three paths. */
  double * total_all = total;
  double * emitted_all = emitted;
  double * counted_all = counted;
  const unsigned mask_cells = inverse_velocity_count * time_count;

#ifdef USE_MPI
  if (mpi_node_count > 1) {
    /* Reduced into scratch rather than in place. mc_MPI_Sum copies its answer back over the
     * buffer it was handed, and SAVE is not called once per run: SIGUSR2 saves and *resumes*,
     * and finally() then saves again. Reducing the accumulators in place would leave the
     * second save reducing already-reduced data, once per node all over again. */
    total_all = save_total ? (double *) calloc(mask_cells, sizeof(double)) : NULL;
    emitted_all = save_emitted ? (double *) calloc(mask_cells, sizeof(double)) : NULL;
    counted_all = save_or_verify_count ? (double *) calloc(mask_cells, sizeof(double)) : NULL;
    if ((save_total && !total_all) || (save_emitted && !emitted_all) || (save_or_verify_count && !counted_all)) {
      /* Out of memory is this node's own news, so it says so itself rather than leaving it
       * to master, and aborts the job rather than calling exit -- which under MPI means
       * MPI_Finalize, and finalizing alone while the others are still working hangs them. */
      fprintf(stderr, "%s: out of memory reducing %u mask cells on node %i\n",
              NAME_CURRENT_COMP, mask_cells, mpi_node_rank);
      MPI_Abort(MPI_COMM_WORLD, -1);
    }
    if (save_total) memcpy(total_all, total, mask_cells * sizeof(double));
    if (save_emitted) memcpy(emitted_all, emitted, mask_cells * sizeof(double));
    if (save_or_verify_count) memcpy(counted_all, counted, mask_cells * sizeof(double));
    /* Summed and not averaged. The ray weights already carry w_stat = 1/mcget_ncount(), and
     * ESS_butterfly reads that in INITIALIZE -- which runs before mccode_main slices ncount
     * across the nodes -- so every node weights by one over the *whole* run's count while
     * emitting its own share of the rays. The sum over nodes is the whole-run answer, which
     * is why mcdetector_import divides by mpi_node_count only for a detector carrying no
     * count array, and ours carries one. */
    if (save_total) mc_MPI_Sum(total_all, mask_cells);
    if (save_emitted) mc_MPI_Sum(emitted_all, mask_cells);
    if (save_or_verify_count) mc_MPI_Sum(counted_all, mask_cells);
  }
#endif

  /* Three independent facts, and a run is no longer obliged to produce any of them, so each
   * clause is written only when the array behind it was filled. The transmission in
   * particular reads off `total`, which `save_total=0` leaves empty in a serial run and NULL
   * in a reduced one -- reporting it unconditionally is a zero in the first case and a
   * dereference of nothing in the second.
   *
   * The acceptance is the factor the resampling multiplies every ray weight by; it exists
   * only when there is a sampler to have computed it.
   *
   * The sampled fraction is that same acceptance measured, from `counted`, which carries no
   * weight on purpose: the factor counts draws rather than intensity. It comes from the mask
   * geometry and never from the run, so this is an independent check on it rather than a
   * restatement. The acceptance itself is not reduced and must not be -- every node computes
   * the same number, and summing it would scale every ray weight in the run by the node
   * count. Only the measured figure is a node's own share.
   *
   * The draw count rides along with it because it is the cheapest check that the reduction
   * happened at all: it is the whole run's ncount, never a single node's slice of it. */
  MPI_MASTER(
    char report[320];
    char clause[128];
    report[0] = '\0';
    if (resample) {
      snprintf(clause, sizeof(clause), "mask acceptance %.6g; ", mask_sampler.acceptance);
      strncat(report, clause, sizeof(report) - strlen(report) - 1);
    }
    if (verify_count) {
      double draws = 0.0;
      for (unsigned i = 0; i < mask_cells; ++i) draws += counted_all[i];
      snprintf(clause, sizeof(clause), "sampled %.6g over %g draws; ",
               chopper_unmasked_probability(counted_all, mask, inverse_velocity_count, time_count),
               draws);
      strncat(report, clause, sizeof(report) - strlen(report) - 1);
    }
    if (save_total) {
      snprintf(clause, sizeof(clause), "weighted transmission %.6g; ",
               chopper_unmasked_probability(total_all, mask, inverse_velocity_count, time_count));
      strncat(report, clause, sizeof(report) - strlen(report) - 1);
    }
    const size_t report_length = strlen(report);
    if (report_length > 2) {
      report[report_length - 2] = '\0';  /* the trailing separator of the last clause */
      printf("%s: %s\n", NAME_CURRENT_COMP, report);
    }

    // dirname is a McCode defined static global variable, and MC_PATHSEP_S is a McCode defined macro for the path separator as a string literal
    if (save_mask) {
      chopper_write_mask_to_file(dirname, filename, ".mask", MC_PATHSEP_S, mask, inverse_velocity_count, time_count, inverse_velocity_edges, time_edges);
    }
    if (save_total) {
      chopper_write_total_to_file(dirname, filename, ".total", MC_PATHSEP_S, total_all, inverse_velocity_count, time_count, inverse_velocity_edges, time_edges);
    }
    // What actually left the moderator, on the same grid: the same picture as `.total` with
    // `resample=0`, and the allowed region alone with `resample=1`.
    if (save_emitted) {
      chopper_write_total_to_file(dirname, filename, ".emitted", MC_PATHSEP_S, emitted_all, inverse_velocity_count, time_count, inverse_velocity_edges, time_edges);
    }
    if (save_count) {
      chopper_write_total_to_file(dirname, filename, ".count", MC_PATHSEP_S, counted_all, inverse_velocity_count, time_count, inverse_velocity_edges, time_edges);
    }
  );

#ifdef USE_MPI
  if (mpi_node_count > 1) {
    if (save_total) free(total_all);
    if (save_emitted) free(emitted_all);
    if (save_or_verify_count) free(counted_all);
  }
#endif
%}
FINALLY %{
  if (total) free(total);
  if (emitted) free(emitted);
  if (counted) free(counted);
  if (mask) free(mask);
  if (inverse_velocity_edges) free(inverse_velocity_edges);
  if (time_edges) free(time_edges);
  chopper_mask_sampler_free(&mask_sampler);
%}
END