fjebaker / SpectralFitting.jl

✨🛰 Fast and flexible spectral fitting in Julia.
https://fjebaker.github.io/SpectralFitting.jl
MIT License
9 stars 5 forks source link

Unfolded data #82

Open fjebaker opened 4 months ago

fjebaker commented 4 months ago

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

phajy commented 1 month ago

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 commented 1 month ago

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:

$$ 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.

$$ \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 commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

$$ 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

https://github.com/fjebaker/SpectralFitting.jl/blob/94bcb07d7e894716de0f94f75da795f2391dc0df/src/datasets/response.jl#L62-L72

(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 commented 1 month ago

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.