Jungfraujoch/image_analysis/SensorAbsorption.h

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

#ifndef JUNGFRAUJOCH_SENSORABSORPTION_H
#define JUNGFRAUJOCH_SENSORABSORPTION_H

#include <cmath>
#include <string>

// How much of the beam a flat sensor stops, as a function of the angle at which the beam enters it.
//
// A photon arriving at incidence angle alpha to the sensor normal crosses t/cos(alpha) of sensor
// instead of t, so more of it is absorbed: the detector is MORE efficient at high angle than at
// normal incidence, and an uncorrected reflection out at the detector edge reads high. Dividing by
// that ratio is the whole correction.
//
//   a(alpha) = t / (L cos alpha)          crossing length in attenuation lengths
//   QE(alpha) = 1 - exp(-a)               absorbed fraction
//   correction to I = QE(0) / QE(alpha)   normalised so normal incidence is untouched
//
// Two limits, both physical and both reached in this corpus:
//   - L << t (long wavelength): every photon stops in the entrance skin, QE = 1 at every angle and
//     the correction is exactly 1. This is why the low-energy case is inert rather than gated off.
//   - L >> t (thin sensor, hard X-rays): absorption is proportional to path length and the
//     correction tends to cos(alpha) - a factor of 2 at 2theta = 60 degrees. This is the regime
//     small-molecule work at 18-25 keV lives in.
//
// alpha is the angle to the DETECTOR NORMAL, not the scattering angle. The two coincide only on an
// untilted detector. On a tilted one the correction acquires an azimuthal dependence at fixed
// resolution, which is the only part of it that is not degenerate with an overall Wilson B.
//
// Attenuation coefficients: NIST X-Ray Mass Attenuation Coefficients (Hubbell & Seltzer), total
// mu/rho with coherent scattering. Photoelectric absorption dominates over the whole range used
// here, and because the correction is a RATIO of two absorbed fractions the photoelectric branching
// ratio cancels exactly; only the attenuation length enters. Coherent scattering, which attenuates
// without converting, is a few per cent of mu/rho in Si below 20 keV and is kept in the total - it
// biases L slightly short, and the correction depends on L only through t/L.
//
// NOT modelled, and negligible where this is used: K-fluorescence escape from the sensor (in CdTe
// above the Cd K edge at 26.71 keV a fluorescence photon can leave the sensor, so the charge is
// recorded in the wrong pixel or not at all); charge sharing between pixels; and the obliquity of
// the entrance window. Above the Cd/Te K edges the escape term is large and this model should not
// be trusted as it stands.
namespace sensor_absorption {

// NIST mass attenuation coefficient mu/rho [cm^2/g] against photon energy [keV]. Absorption edges
// are the repeated energies - the table is walked in order and interpolated log-log, which is how
// these tables are meant to be read. Si has only the K edge at 1.839 keV, far below anything used.
struct MuRow { double E_keV, mu_rho; };

inline constexpr MuRow SI_TABLE[] = {
    {1.0, 1570.0}, {1.5, 535.5}, {1.8389, 309.2}, {1.8389, 3192.0}, {2.0, 2777.0},
    {3.0, 978.4}, {4.0, 452.9}, {5.0, 245.0}, {6.0, 147.0}, {8.0, 64.68}, {10.0, 33.89},
    {15.0, 10.34}, {20.0, 4.464}, {30.0, 1.436}, {40.0, 0.7012}, {50.0, 0.4385}};

inline constexpr MuRow CD_TABLE[] = {
    {1.0, 7350.0}, {1.5, 2931.0}, {2.0, 1473.0}, {3.0, 541.4}, {3.5375, 357.5},
    {3.5375, 1152.0}, {3.63101, 1083.0}, {3.727, 1013.0}, {3.727, 1389.0}, {4.0, 1170.0},
    {4.018, 1157.0}, {4.018, 1324.0}, {5.0, 768.5}, {6.0, 479.3}, {8.0, 225.4}, {10.0, 124.4},
    {15.0, 41.78}, {20.0, 19.20}, {26.7112, 8.809}, {26.7112, 50.65}, {30.0, 37.65},
    {40.0, 17.78}, {50.0, 9.779}};

inline constexpr MuRow TE_TABLE[] = {
    {1.0, 8434.0}, {1.5, 3608.0}, {2.0, 1832.0}, {3.0, 679.2}, {4.0, 329.7}, {4.3414, 267.8},
    {4.3414, 788.2}, {4.47465, 750.4}, {4.612, 699.5}, {4.612, 944.5}, {4.7728, 878.2},
    {4.9392, 806.2}, {4.9392, 929.2}, {5.0, 901.4}, {6.0, 572.1}, {8.0, 270.2}, {10.0, 150.1},
    {15.0, 50.78}, {20.0, 23.41}, {30.0, 7.878}, {31.8138, 6.738}, {31.8138, 37.19},
    {40.0, 20.64}, {50.0, 11.45}};

// Log-log interpolation, clamped at both ends of the table.
template <std::size_t N>
double InterpMuRho(const MuRow (&tab)[N], double E_keV) {
    if (E_keV <= tab[0].E_keV)
        return tab[0].mu_rho;
    if (E_keV >= tab[N - 1].E_keV)
        return tab[N - 1].mu_rho;
    std::size_t i = 1;
    while (i < N - 1 && tab[i].E_keV < E_keV)
        i++;
    const double e0 = tab[i - 1].E_keV, e1 = tab[i].E_keV;
    if (!(e1 > e0))                    // the two rows of an absorption edge: take the upper side
        return tab[i].mu_rho;
    const double f = (std::log(E_keV) - std::log(e0)) / (std::log(e1) - std::log(e0));
    return std::exp(std::log(tab[i - 1].mu_rho) + f * (std::log(tab[i].mu_rho) - std::log(tab[i - 1].mu_rho)));
}

// Attenuation length 1/mu [um] of a sensor material at a given wavelength. Si and CdTe are the
// sensors in use; an unrecognised material is treated as silicon, which is what DetectorSetup
// defaults to anyway.
inline double AttenuationLength_um(const std::string &material, double lambda_A) {
    if (!(lambda_A > 0.0))
        return 0.0;
    const double E_keV = 12.39842 / lambda_A;
    double mu_rho_cm2_g, rho_g_cm3;
    if (material == "CdTe") {
        // Mass fractions from the atomic weights, Cd 112.414 and Te 127.60.
        constexpr double w_cd = 112.414 / 240.014, w_te = 127.60 / 240.014;
        mu_rho_cm2_g = w_cd * InterpMuRho(CD_TABLE, E_keV) + w_te * InterpMuRho(TE_TABLE, E_keV);
        rho_g_cm3 = 5.85;
    } else {
        mu_rho_cm2_g = InterpMuRho(SI_TABLE, E_keV);
        rho_g_cm3 = 2.3290;
    }
    const double mu_cm = mu_rho_cm2_g * rho_g_cm3;
    return mu_cm > 0.0 ? 1e4 / mu_cm : 0.0;
}

// Everything the per-reflection correction needs, reduced to the two numbers that are constant for
// a dataset. Built once on the host; the GPU predictor takes the two floats.
struct SensorQE {
    float a0 = 0.0f;        // t / L, the optical thickness at normal incidence
    float qe0 = 1.0f;       // absorbed fraction at normal incidence
    bool active = false;    // false where the sensor is opaque and the correction is exactly 1

    // Opaque beyond this: exp(-20) = 2e-9, so QE(0)/QE(alpha) rounds to exactly 1.0f for every
    // incidence angle and the correction is bit-identical to not applying it. This is what makes
    // the long-wavelength case inert without a flag or a threshold anyone has to choose.
    static constexpr float OPAQUE_A0 = 20.0f;

    static SensorQE Build(const std::string &material, double thickness_um, double lambda_A) {
        SensorQE q;
        const double L = AttenuationLength_um(material, lambda_A);
        if (!(thickness_um > 0.0) || !(L > 0.0))
            return q;
        q.a0 = static_cast<float>(thickness_um / L);
        q.qe0 = static_cast<float>(1.0 - std::exp(-q.a0));
        q.active = q.a0 < OPAQUE_A0;
        return q;
    }

    // The multiplicative correction to an intensity recorded at incidence angle alpha: QE(0)/QE(alpha).
    // Always <= 1, because a sensor is more efficient off-normal than head-on.
    [[nodiscard]] float Factor(float cos_alpha) const {
        if (!active || !(cos_alpha > 1e-3f))
            return 1.0f;
        const float qe = 1.0f - std::exp(-a0 / cos_alpha);
        return qe > 0.0f ? qe0 / qe : 1.0f;
    }
};

}   // namespace sensor_absorption

#endif   // JUNGFRAUJOCH_SENSORABSORPTION_H