Unfolded data

[fjebaker]

the data points plotted are calculated by D*(unfolded model)/(folded model), where D is the observed data,

• XSPEC manual

[phajy]

See also Mike Nowak's conference proceeding, X-ray and Radio Monitoring of GX 339-4 and Cyg X-1, especially §3, "A Rant on the Nature of Evil". I think we should implement two ways of unfolding spectra. A more "proper" way (see Mike's conference proceeding) and the XSPEC way to allow comparison with the literature.

@fjebaker - since we won't use unfolded spectra for fitting (although perhaps that should just be ill advised rather than forbidden?) I was thinking of putting this in as a plot recipe. What do you think? It does require access to the data, responses, and model, but I think that is now all bundled up in the fitting result?

[fjebaker]

Sounds good, I hadn't seen that Nowak proceeding before. Just for my own comparison, I write here how XSPEC and ISIS do their unfolded spectra and model:

• XSPEC

$$ D \times \frac{M}{R \times M}, $$

where $D$ is the data, $M$ is the model, and $R$ is the response (with ancillary folded in). That is, the data multiplied by the model over the folded model.

• ISIS

$$ \frac{C - B}{\Delta t \int R(E) \text{d} E}, $$

where $C$ are the counts, $B$ the background, $\Delta t$ is the integrated exposure time, and $R$ is now a unit normalized response matrix (with ancillary). We have a function for normalizing $R$, but AFAIK I haven't found a response matrix yet that hadn't already been normalized).

I read the integral as saying the response matrix summed (integrated) over the photon energy $E$, as in, the total probability of detecting in that channel? Which is why the response matrix must be unit normalized.

[phajy]

The response matrices are typically unit normalised but this is not always the case and I have run into two different cases: 1) response matrices that are nearly unit normalised (not sure why such things should exist), and 2) (much more commonly) combined RMF and ARF (sometimes given the extension .rsp) which can be created with, e.g., the FTOOL marfrmf.

ISIS has factor_rsp command that returns a normalised RMF and creates an appropriate ARF, designed to split apart response files that are RMF $\times$ ARF.

Typo: the integrand in the denominator of the ISIS expression should contain the term $A(E)$.

[fjebaker]

I don't understand how ISIS can factor out the ARF. I thought merging the RMF and ARF was a matrix multiplication

$$ \tilde{R} = A^{\text{T}} R, $$

so how can you invert this to find both $R$ and $A$ without additional information about $R$ or $A$?

[phajy]

Ah, right, a key point here is that the ARF is not a general matrix - it is vector (or diagonal matrix). Take a look at The Calibration Requirements for Spectral Analysis. So this means you can, at a given energy (or channel), normalise the combined response (ARF and RMF in one matrix) $R$, so that

$$ \begin{aligned} A(E_i) &= \sum_j R(E_i, E_j) \cr R^\prime(E_i, E_j) &= \frac{R(E_i, E_j)}{A(E_i)} \end{aligned} $$

where $A(E_i)$ is the effective area at energy $E_i$, $R^\prime(E_i, E_j)$ is the probability that a photon of energy $E_i$ is detected in channel $E_j$ and is unit normalised at each energy $E_i$.

This procedure will take $R$ and give you an ARF, $A$, and unit normalised RMF, $R^\prime$.

P.S. I note that the documentation says "This format is very similar to that used originally by the XSPEC .rsp standard-format response files" which suggests that XSPEC pre-dates the establishment of the ARF and RMF as response files.

P.P.S. I also just spotted the ftools chkarf, chkpha, and chkrmf - although the later has been replaced with ftchkrmf. The instructions do give some clues but you might need to look at the source code to see what these really do in full. Might be good consistency checks to perform when reading in PHA, RMF, and ARF files (something for the future!)

[fjebaker]

$$ A(E_i) = \sum_j R(E_i, E_j) $$

This was the relation I was missing! It all clicks into place now. Yes we can quite easily add those implementations into SpectralFitting.jl. I have this currently unused function in the source

SpectralFitting.jl/src/datasets/response.jl

Lines 62 to 72 in 94bcb07

	# todo: currently unused; do we need this?
	function normalise_rows!(matrix)
	# returns weights
	@views map(1:size(matrix, 2)) do
	w = sum(matrix[i, :])
	if ΣR > 0.0
	matrix[i, :] .= matrix[i, :] / w
	end
	w
	end
	end

(Edit: this code has an undefined symbol so doesn't actually run correctly, but it's a trivial fix)

which normalises each row of the response matrix to 1, which is in essence your $R^\prime$ I believe? Pulling this apart to also then calculate the ARF is simple.

We can (and probably should) also add functions that check that the matrix is unit normalized when loading in and issue a warning to the user if not, with instructions on how to normalize the matrix.

[fjebaker]

Infact the weights in the above function are already the ARF as in your equation, so we actually just need a wrapper function with a better name and expose it publicly.