Subtract backgrounds

[fjebaker]

With #20, there should be a further overhaul to let SpectralDataset hold onto most parts of OGIP, and flatten out some of the common meta fields (quality, grouping, etc.), and also hold onto the background spectrum, so this may be more easily subtracted.

[phajy]

Can you update spectral-datasets.jl and xmm-newton-mission.jl (and probably elsewhere) to include the background file. When I load an XMM-Newton data set it loads the spectrum, RMF, and ARF as specified in the spectrum header, but not the background file.

I think that simply subtracting the background from the source spectrum is the most straightforward approach in the first instance. But you might want to allow different backgrounds to be used in the future (e.g., an analytic model of the background or whatever).

[fjebaker]

@phajy Have been working on this today. I am in the progress of a larger parser rewrite since #45, so this should hopefully not take too long.

[phajy]

Sounds good. I was just writing this when your reply appeared!

Looking at ogip.jl I'm not sure what is the best way to handle the background. For now it might be good to assume/check the background spectrum is from the same telescope and instrument with the same number of channels so the background subtraction can be done channel by channel and then grouped. I think this is reasonable because it will cover almost all cases. The background exposure time will likely be the same (although it doesn't have to be) but source and background areas will be different (i.e., they will have different BACKSCAL keywords). The background could even be from a different observation (although probably not in the cases we'll be looking at initially).

[fjebaker]

Yes, that sounds good. From the XSPEC manual I have the method for background-subtracted counts $C(I)$ as

$$ C(I) = \frac{D(I)}{a_D(I) t_D} - \frac{b_D(I)}{b_B(I)} \frac{B(I)}{a_B(I) t_B}, $$

where the exposure times is $t_i$ of the data $D$ and background $B$, with area and background scaling factors $a_i$ and $b_i$ respectively.

From looking at the headers of the XMM and NuSTAR samples I am using in the test suite (see here), everything is defined, so that should be simple enough. On that note, if you have any other CC spectra that I could add to the test suite... that would help tremendously with avoiding future errors :)

I'm not going to check the instrument / telescope match currently, and will instead just assume that the user knows what they are doing with the background files. I can add a warning that will complain if they differ, but I am a little cautious about enforcing any rules given some of the edge cases I want to be able to address.

[fjebaker]

Actually looking at that equation, it must be the background-subtracted rates, since it is being divided by the exposure time. I must have written down the wrong symbol on the LHS.

[fjebaker]

After

• 48

then background subtracting can be implemented. I think the easiest way to do it is when we assemble the target vector for fitting.

Will add a background_subtracted_target_vector method for this. This can then be extended without too much difficulty to support other background models, depending on the type of the background in SpectralDataset. This is currently fixed to Union{Missing,Spectrum} but it could be extended for some analytic background type, in which case only background_subtracted_target_vector would need to updated for the fitting to still be able to work okay.