Jungfraujoch/common/Reflection.h

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <cstdint>
#include <optional>
#include <cmath>

#include "SpotToSave.h"

struct Reflection {
    int32_t h;
    int32_t k;
    int32_t l;
    float image_number; // Can be in-between for 3D integration
    float delta_phi_deg; // phi angle from XDS - difference from middle of current frame (NOT an absolute angle)
    float predicted_x;
    float predicted_y;
    float observed_x;
    float observed_y;
    float d;
    float I;
    float bkg;
    float var_bkg;        // non-signal (background) part of sigma^2, carried to the merge
    float sigma;
    float dist_ewald;     // stills: signed excitation error |S| - 1/lambda (1/A), positive outside the sphere
    // The reciprocal Lorentz factor (rotation only - a still's Lorentz factor is one) times the
    // reciprocal polarization factor, and nothing else. This is what LP means everywhere in the
    // field, and it is what the CBOR key "rlp" and the HDF5 dataset "lp" store the reciprocal of.
    // What it is not is a scale: the fitted per-image scale and the partiality stay out of it and
    // are divided in separately below. (Named after DIALS's prescaling_correction.)
    float prescaling_corr;
    // The sensor's angle-dependent efficiency, QE(0)/QE(alpha): always <= 1, and exactly 1 where the
    // sensor is opaque or its material and thickness are unknown. It is carried BESIDE
    // prescaling_corr rather than inside it, because LP and detector response are two different
    // things and every file the field reads keeps them apart. It is one of the three factors whose
    // product is the total deterministic correction, with prescaling_corr above and flight_corr
    // below. Defaulted to 1 so a reflection read from a file written before this existed is a no-op
    // rather than a zero.
    float qe_corr = 1.0f;
    // The medium in the sample-to-pixel flight path, exp(D/L * (1/cos(alpha) - 1)) with D the
    // normal-incidence distance, L the medium's attenuation length and alpha the angle of incidence
    // on the detector: always >= 1, because an oblique reflection crossed more of the medium than
    // one arriving head-on. Exactly 1 under --flight-path vacuum. Carried beside the two above for the same
    // reason they are carried apart - it is neither beam geometry nor detector response but the
    // medium in between, and unlike either of them it is set by the flight distance. The total
    // deterministic correction on a reflection is prescaling_corr * qe_corr * flight_corr, and every
    // site that corrects an intensity multiplies all three. Defaulted to 1 so a reflection read
    // from a file written before this existed is a no-op rather than a zero.
    float flight_corr = 1.0f;
    float partiality;     // fraction of the reflection recorded in the sampled (rocking) slice
    float zeta;
    float image_scale_corr; // I_true = image_scale_corr * I; = prescaling_corr * qe_corr * flight_corr / (partiality * image_scale)
    bool observed = false;
    bool on_ice_ring = false; // sits on a hexagonal-ice powder ring: excluded from scaling; merged on stills, ring by ring on rotation (DecideIceRings)
    // The signal disk lost pixels to the mask or to saturation, so the intensity is the profile's
    // estimate over what remained. Such a measurement may be low, and the merge does not let it
    // testify against a larger observation of the same reflection (WilsonOutliers.h).
    bool clipped = false;
    // A pixel of the signal disk was saturated. The reflection is then not a measurement of its
    // intensity at all - its brightest pixels are missing - and a rotation merge drops the whole
    // rocking event it belongs to, as XDS does with an overloaded reflection.
    bool overloaded = false;
};

// One full reflection of a rotation sweep as the merge used it - its partials summed into one
// measurement, every correction and the per-frame scale applied, its sigma under the error model the
// merge weighted it with - but NOT averaged with its symmetry mates. h k l are the index it was
// measured at, not its ASU representative. Written as the SHELX HKLF 4 file, so that the program
// reading it sees the equivalents and computes Rint itself, and as the intensities of the unmerged MTZ.
struct ScaledFull {
    int32_t h = 0;
    int32_t k = 0;
    int32_t l = 0;
    float I = NAN;
    float sigma = NAN;
    float image_number = NAN;   // the frame of the event's peak partial, which tells two events of one index apart
    // The full's counting variance as the merge rebuilds it at an intensity I, before its error model:
    // var_bkg + var_per_I * I + (capture * I)^2, all at the intensity's scale (what a later merge rebuilds
    // the full from). NaN where not exported.
    float var_bkg = NAN;
    float var_per_I = NAN;
    float capture = NAN;
    float d = NAN;              // resolution, A
    // The diffracted beam's direction in the crystal: its components along the real-space cell vectors
    // (a.s1, b.s1, c.s1), so that s1 = sum_i (a_i.s1) a*_i - reindexed like h, k, l. What a joint merge of
    // several sweeps fits one absorption surface over. NaN where not known.
    float s1_cell[3] = {NAN, NAN, NAN};
};

struct MergedReflection {
    int32_t h = 0;
    int32_t k = 0;
    int32_t l = 0;
    float I = NAN;
    float sigma = NAN;
    float I_half[2] = {NAN, NAN};
    float sigma_half[2] = {NAN, NAN};
    // Half-set means with every frame's observations kept whole in one half (stills: R_split and the
    // CrystFEL-style CC1/2). NaN where a half saw none, and on rotation.
    float I_split[2] = {NAN, NAN};
    // Weight of this reflection in a CC1/2: the precision its half-sets would have had with every
    // observation at the run's typical frame scale, over the precision they have. 1 on a sweep
    // without a weak stretch; small for a reflection measured only where the crystal barely
    // diffracted, whose scaled-up noise would otherwise count as much as a well-measured pair.
    // Internal, not reported: it enters the CC1/2 the processing decides on (the automatic resolution
    // cut, MergeStatisticsShell::cc_half_internal), never the reported CC1/2, which is the plain Pearson.
    float cc_weight = 1.0f;
    float d = 0.0;
    // Any observation of this reflection sat on an ice ring. The intensity is still merged and
    // written - the ring contaminates it, it does not make it absent - and this only marks it so
    // a consumer that must not read the ring as crystal signal can leave it out.
    bool on_ice_ring = false;
    bool rfree_flag = false;
    float F = NAN;       // French-Wilson amplitude |F| (filled by ApplyFrenchWilson at end of merge)
    float sigmaF = NAN;  // its sigma
    // Anomalous (Bijvoet) split of this reflection's own observations, kept even when the merge is
    // Friedel-averaged (I above is the Friedel mean). Lets I(+)/I(-) be written and CCano reported by
    // default without scaling anomalously; NaN when a hand was not measured or for centrics.
    float I_plus = NAN;
    float sigma_plus = NAN;
    float I_minus = NAN;
    float sigma_minus = NAN;
    // French-Wilson amplitudes of the two hands (filled by ApplyFrenchWilson from I_plus/I_minus).
    float F_plus = NAN;
    float sigmaF_plus = NAN;
    float F_minus = NAN;
    float sigmaF_minus = NAN;
};