Instrument cross calibration¶

Most of the time, when dealing with multiple X-ray instruments, one would want to account for the differences in calibrations between all observatories. The straightforward way to do this with jaxspec is to use an InstrumentModel. Let's first load observations from PN, MOS1 and MOS2.

In [ ]:

%%capture
# Hide the output of this cell
import numpyro

numpyro.enable_x64()
numpyro.set_platform("cpu")
numpyro.set_host_device_count(4)

import numpyro.distributions as dist
from jaxspec.fit import VIFitter
from jaxspec.data.util import load_example_obsconf
from jaxspec.model.additive import Blackbodyrad, Powerlaw
from jaxspec.model.multiplicative import Tbabs
from jaxspec.model.background import BackgroundWithError

spectral_model = Tbabs()*(Powerlaw() + Blackbodyrad())

obsconfs = load_example_obsconf("NGC7793_ULX4_ALL")

prior = {
    "spectrum.powerlaw_1.alpha": dist.Uniform(1, 3),
    "spectrum.powerlaw_1.norm": dist.LogUniform(1e-5, 1e-3),
    "spectrum.blackbodyrad_1.kT": dist.Uniform(0, 2),
    "spectrum.blackbodyrad_1.norm": dist.LogUniform(1e-2, 1),
    "spectrum.tbabs_1.nh": dist.Uniform(0, 1)
}

%%capture # Hide the output of this cell import numpyro numpyro.enable_x64() numpyro.set_platform("cpu") numpyro.set_host_device_count(4) import numpyro.distributions as dist from jaxspec.fit import VIFitter from jaxspec.data.util import load_example_obsconf from jaxspec.model.additive import Blackbodyrad, Powerlaw from jaxspec.model.multiplicative import Tbabs from jaxspec.model.background import BackgroundWithError spectral_model = Tbabs()*(Powerlaw() + Blackbodyrad()) obsconfs = load_example_obsconf("NGC7793_ULX4_ALL") prior = { "spectrum.powerlaw_1.alpha": dist.Uniform(1, 3), "spectrum.powerlaw_1.norm": dist.LogUniform(1e-5, 1e-3), "spectrum.blackbodyrad_1.kT": dist.Uniform(0, 2), "spectrum.blackbodyrad_1.norm": dist.LogUniform(1e-2, 1), "spectrum.tbabs_1.nh": dist.Uniform(0, 1) }

An InstrumentModel can be constructed from a ConstantGain and a ConstantShift. Their priors live in the unified prior dict under instrument.gain.factor and instrument.shift.offset. For instance, in the following observations, we will use the PN one as the reference, and specify it using "PN" because it is the key corresponding to this observation in the dictionary of observations we use.

In [ ]:

import numpyro.distributions as dist
from jaxspec.model.background import BackgroundWithError
from jaxspec.model.instrument import InstrumentModel, ConstantGain, ConstantShift

instrument_prior = {
    **prior,
    # `[*]` expands to one independent draw per applicable observation
    # (= each instrumented obs; here MOS1 and MOS2).
    "instrument.gain.factor[*]":  dist.Uniform(0.5, 1.5),   # 50–150% of nominal calibration
    "instrument.shift.offset[*]": dist.Uniform(-0.3, +0.3), # ±300 eV shift
}

fitter = VIFitter(
    spectral_model, instrument_prior, obsconfs,
    # PN is the explicit reference: ``None`` skips instrument calibration on PN.
    # Each MOS observation gets its own InstrumentModel; the per-obs gain/shift
    # values are sampled independently via the ``[*]``-suffixed priors.
    instrument_model={
        "PN":   None,
        "MOS1": InstrumentModel(gain=ConstantGain(), shift=ConstantShift()),
        "MOS2": InstrumentModel(gain=ConstantGain(), shift=ConstantShift()),
    },
    background_model=BackgroundWithError(),
)

import numpyro.distributions as dist from jaxspec.model.background import BackgroundWithError from jaxspec.model.instrument import InstrumentModel, ConstantGain, ConstantShift instrument_prior = { **prior, # `[*]` expands to one independent draw per applicable observation # (= each instrumented obs; here MOS1 and MOS2). "instrument.gain.factor[*]": dist.Uniform(0.5, 1.5), # 50–150% of nominal calibration "instrument.shift.offset[*]": dist.Uniform(-0.3, +0.3), # ±300 eV shift } fitter = VIFitter( spectral_model, instrument_prior, obsconfs, # PN is the explicit reference: ``None`` skips instrument calibration on PN. # Each MOS observation gets its own InstrumentModel; the per-obs gain/shift # values are sampled independently via the ``[*]``-suffixed priors. instrument_model={ "PN": None, "MOS1": InstrumentModel(gain=ConstantGain(), shift=ConstantShift()), "MOS2": InstrumentModel(gain=ConstantGain(), shift=ConstantShift()), }, background_model=BackgroundWithError(), )

Now lets fit the data as usual

In [ ]:

result = fitter.fit(num_steps=30_000, num_samples=10_000)

result = fitter.fit(num_steps=30_000, num_samples=10_000)

In [ ]:

result.plot_ppc();

result.plot_ppc();

In [ ]:

result.plot_corner()

result.plot_corner()

The associated calibrations can be found in the result object and accessed as follows

In [ ]:

import numpy as np
import matplotlib.pyplot as plt

plt.close('all')
plt.subplot(121)
plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS1.gain.factor']), label='Gain MOS1', alpha=0.5, bins=50, density=True)
plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS2.gain.factor']), label='Gain MOS2', alpha=0.5, bins=50, density=True)
plt.title('Fitted Gain')
plt.legend()

plt.subplot(122)
plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS1.shift.offset']), label='Shift MOS1', alpha=0.5, bins=50, density=True)
plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS2.shift.offset']), label='Shift MOS2', alpha=0.5, bins=50, density=True)
plt.title('Fitted Shift')
plt.legend()
plt.show()

import numpy as np import matplotlib.pyplot as plt plt.close('all') plt.subplot(121) plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS1.gain.factor']), label='Gain MOS1', alpha=0.5, bins=50, density=True) plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS2.gain.factor']), label='Gain MOS2', alpha=0.5, bins=50, density=True) plt.title('Fitted Gain') plt.legend() plt.subplot(122) plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS1.shift.offset']), label='Shift MOS1', alpha=0.5, bins=50, density=True) plt.hist(np.ravel(result.inference_data.posterior['forward.instrument.MOS2.shift.offset']), label='Shift MOS2', alpha=0.5, bins=50, density=True) plt.title('Fitted Shift') plt.legend() plt.show()