pub fn fit_counts(
measurement: &Measurement,
calibration: &Calibration,
) -> Result<CountsFit, PipelineError>Expand description
Fit the areal densities and temperatures of the regions’ materials, each
region’s background terms and gamma transmission, the normalization, the
time background, and the timing offset, flight path and pulse numbers of
calibration, that are not known, to the
raw counts of both runs in every region. On a uniform grid of flight
times u_i, energies E_i and step w, the counts of region r in bin
k are
O_rk = ℓ^O_rk · [w Σ_i φ_ri P_ki + p_r G_k + p_r ρ F^O w_k]
S_rk = ℓ^S_rk · [c_q · a · (w Σ_i φ_ri [T_ri + b_r(E_i)] P_ki + γ_r p_r G_k)
+ p_r ρ F^S w_k]
T_ri = (1 + v_r D_ri)^(−1/v_r), D_ri = Σ_m n_rm σ_m(E_i; T_r), or 1 without a material
b_r(E) = b0_r + b1_r/√E + b2_r·√E + d_r·exp(−f_r/√E)
G_k = ∫ over bin k of (C + A·exp(−(t − t₁)/τ)) dtwith φ_r the region’s beam per µs, P_ki the chance of a neutron at
u_i arriving in bin k, n_rm and σ_m each isotope’s density and
total cross section Doppler-broadened at the material’s temperature T_r,
c_q the charge ratio, a the normalization and ℓ^O_rk, ℓ^S_rk each
run’s live fraction. The background is beam neutrons that reach the
detector another way, so it passes through the pulse and scales with the
normalization; SAMMY’s BackA, BackB, BackC, BackD and BackF
(cro/mnrm1.f90) are a·b0, a·b1, a·b2, a·d and f, to within the
background’s change over the pulse’s delay. The time
background G is counts that follow the proton charge but reach the
detector without passing through the moderator pulse, such as the prompt
gammas of each proton pulse and their tail, at the time t of the edges’
clock, t₁ the first edge: it is not folded with the pulse, and the sample
run scales it with c_q·a, which correct the charge ratio and an
efficiency change alike in every pixel, and the fraction γ_r that reaches
the region’s pixels, p_r the region’s pixels. It turns
the ratio R of the neutrons alone into R + g·(γ_r − R), g its share of
the open-beam counts; SAMMY’s backgrounds in time (Bgfrpi) are added to
the transmission instead. The room background ρ (RoomBackground)
follows each run’s duration instead, F^O and F^S its frames and w_k
bin k’s width: it is neither folded with the pulse nor scaled by the
charge ratio, the normalization or the gamma transmission, so it separates
from the time background’s C only where c_q·a·γ_r differs from
F^S/F^O.
The regions share the normalization, the run pair’s scale; and t0, the
flight path and the pulse, which one clock, flight path and moderator set
for every pixel, and with them the grid and P_ki. Each region has its
own beam, from its own open-beam counts, with the intervals its
fit_open_beam chooses, fitted with
the rest to both runs; and its own
background, which depends on what surrounds its pixels. A region without
a material counts c_q·a·(1 + b_r) of its beam in the sample run, so it
pins the normalization when its b0 is known, as zero in the open-region
normalization of imaging, and nothing about it when its b0 is fitted.
The regions must be separate pixels: the covariance takes every bin as
independent.
The Calibration’s t0 and flight path L are fitted unless known.
The grid is built at a t0₀ and L₀, at first the starting ones. Grid
point i keeps its energy and arrives at t0 + (L/L₀)·u_i, with u_i
its flight time at t0₀ and L₀, standing for (L/L₀)·w of flight time;
the beams are functions of the arrival time less the starting t0. This
is SAMMY’s energy scale: its Tzero is t0, and its Elzero times its
flight path is L (dat/mdat0.f90, Mtzero).
The pulse’s numbers are fitted unless known, each point’s P_ki taken at
the pulse’s laws at its energy E_i. A trial pulse whose α or β is 0
at some energy is outside the model and rejected; a fit driven there ends
unconverged. The open-beam fits that start the beams and weight the
open-beam runs use the starting pulse and, at first, the time background’s
starting numbers; the fit is repeated with them at its fitted numbers, up
to twenty times, until the fitted numbers change the counts they took as
known by Σ_r Σ_k (ΔO_rk)²/O_rk at most BOUND over every region’s
open-beam bins, and from the third time on each region keeps the beam
intervals its second open-beam fit chose.
The grid’s first step is at most half the narrowest Doppler full width at
half maximum, in flight time, of any resonance of any region inside its
energy span, at that region’s starting temperature; it is then halved
until the counts of every run of every region, together, meet
BOUND. The step is uniform and shared, so a
wide window whose span holds a narrow resonance at high energy, or a low
fitted temperature in any region, can exceed the grid’s point cap; a
temperature a region’s counts barely determine can run to 1 K and refuse
the fit that way.
The fit minimizes each run’s half Poisson deviance over its
overdispersion, plus ½((x − value)/sd)² for each Value::Measured
quantity x, and ½(θ − m)ᵀC⁻¹(θ − m) over the pulse numbers θ a
calibration’s Pulse::prior covers,
with m and C that prior’s mean, which may lie past a number’s bound,
and covariance. Each region’s open-beam run is weighted with its
open-beam fit’s overdispersion, and its sample run first with the same.
The fit is repeated from its answer while any sample run’s
overdispersion, measured on the bins outside those left out that the first
fit predicts at least one count, changes by more than 1%, the coarser grid of the accepted pair is
wider than the rule at the fitted temperatures, its flight times miss
some that the fitted t0, L and pulse need, or the bins left out at
those differ from the ones it left out, which for the first fit are those
at the starting ones. The next first grid is the
finer of that pair’s coarser grid and the rule’s grid at the fitted
temperatures, or, when the flight times miss, the rule’s grid of a grid
built at the fitted t0, L and pulse. After twenty fits it is reported
unconverged.
The fitter finds a local minimum. Where a temperature is fitted and the
counts reject the fit, some region’s sample-run overdispersion exceeding 1
by more than three of its standard deviations under Poisson counts,
√(2/ν) over the run’s ν bins outside those left out, the fit is made
again with each fitted temperature started at 2,500 K, or at the bound
nearest it. The second fit is returned when it converged and the first
did not, or when both or neither converged and its half Poisson deviance
of every run’s counts, outside the sample bins either fit leaves out, is
lower by more than NEWTON_DECREMENT_TOL; a second fit that fails is
set aside. Where a black resonance
empties bins and the background is near zero, a fitted background with no
lower bound can drive a black bin’s prediction to zero, and the fit ends
unconverged; with b1 and b2 known, b0 bounded below by 0 ends on
that bound, converged.
§Errors
PipelineError::ShapeMismatch unless each run of each region has one
count, and one live fraction when given, per bin;
PipelineError::InvalidParameter if a count is not a whole non-negative
number, a run has no counts, a live fraction is not in (0, 1], the charge
ratio or a room background’s frames are not finite and positive, a
region’s pixels are not finite and
positive, a region without a material has capture gammas or in-scatter
other than known 0, a spread’s kernel has a term whose weight or length is
not finite and positive, its grid holds no cell or more than 2^20, it does
not give each region one or more cells on the grid, or a cell is in two
regions, a known, starting or measured value is
not finite and in its quantity’s range, a measured value’s sd is not finite
and positive, bounds are not lower < upper in that range
with the start between them, the pulse’s calibration covers a number that
is not fitted or is measured, no region holds a material, a material has
no isotopes or lists one twice, an isotope’s resonance data are not
finite, an isotope not known to be absent has a resonance between the
energies of the last and first time edges, at the starting t0 and flight
path or a converged fit’s, outside the pulse’s
line_span_ev, or the energies
its broadened cross section reads on the grid at the starting t0, flight
path and pulse, or at the fitted ones of any pass that rebuilds it, at the
region’s known temperature or at the upper bound of a fitted one, down to
zero for a window within the thermal spread of zero energy, are not inside
a single one of its evaluated resolved ranges,
or, in a region with in-scatter, an isotope has a resonance below that
grid’s lowest energy and above the lowest energy a neutron scattered forward
leaves with;
PipelineError::UnmodelledCounts if at the fit, converged or not, a bin
is predicted negative or non-finite counts, recorded or, under a spread,
arriving, or, outside the bins left
out, holds counts predicted below NEGLIGIBLE_PREDICTION: starting or known values the fitter cannot
leave;
PipelineError::InvalidParameter if a material’s covariance names no
isotope of it, or one twice, lists no resonance, one twice or one the
isotope does not have, or is not a finite symmetric positive
semi-definite matrix over their energies and widths, or if it is given for
an isotope another region also holds or with a pulse calibration’s prior;
everything fit_open_beam refuses for
any region;
PipelineError::FlightTimeGrid for the grid’s refusals at the starting
t0, flight path and pulse or the fitted ones of any pass, including more
points than it allows; PipelineError::Fitting if the fitter fails, a
region’s open-beam fit does not converge with the beam intervals it keeps,
the cross sections fail at the start or, for a listed resonance parameter
shifted, at the fit, or a decomposition of the pulse’s
consistency test fails; a failure at a trial temperature is a rejected
step.