Skip to main content

fit_counts

Function fit_counts 

Source
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₁)/τ)) dt

with φ_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.