Creating source models¶

Warning

While this is possibly one of the most important pieces of documentation, it is still work in progress. Do not hesitate to complain if something is not clear incomplete, or simply wrong.

Source models are specified in ixpeobssim through Python configuration files—effectively simple Python modules. As we shall see in a second, this is the key of the flexibility in terms of source specification provided by the simulation framework. All the source model configuration files shipped with ixpeobssim are included in the [github]/ixpeobssim/config folder, and that is a good starting point to get up and running with creating your own observation-simulation setup.

In a nutshell, the one thing that a source configuration file must contain is an instance of the ixpeobssim.srcmodel.roi.xROIModel class called ROI_MODEL (mind that the name and case of the variable are important, as when you run a simulation the relevant application will try and import that exact name from the Python configuration module—and fail miserably if that is not possible). In the rest of this section we shall walk through the relevant classes for the definition of a region of interest (ROI) to be simulated.

Spectra and polarization patterns¶

Understanding how photon spectra and polarization patterns are specified in ixpeobssim is key to being able to create a functional simulation setup. The two basic rules, here, are:

the photon spectrum can be an arbitrary function of energy and time, expressed in units of $\text{cm}^{-2}~\text{s}^{-1}~\text{keV}^{-1}$ ;
the polarization degree and angle can be arbitrary functions of energy, time, and position in the sky.

From the implementation standpoint, the photon spectrum and the polarization degree and angle are simply Python functions with the proper signature:

def spec(E, t):
    """Definition of the photon spectrum (cm^-2 s^-1 keV^-1).
    """
    # Function definition here.

def pol_deg(E, t, ra, dec):
    """Definition of the polarization degree (between 0 and 1).
    """
    # Function definition here.

def pol_ang(E, t, ra, dec):
    """Definition of the polarization angle (in radians).
    """
    # Function definition here.

These functions are then used internally by the simulation code to generate a proper photon list, and we shall see in the next section how they are passed to the relevant classes. We emphasize that this approach is fairly flexible, as the user can put in the function bodies literally anything: a constant, an analitic formula, an interpolator constructed from a numeric table, etc. The only important thing is to make sure that at runtime the function can be called with the proper signature (and, by the way, can operate on numpy arrays—so always make sure that you use the numpy version of the mathematical functions that you need in your function body).

Note

At this point you might wonder why the rules are comparatively restrictive for the photon spectra which, in contrast to the polarization patterns, cannot be dependent on the position on the sky. Well, the reason is that a fully arbitrary source model doesn’t really fit with the current implementation of the basic simulation strategy. But before you panic, a couple of things:

for extended sources, as we shall see in the next section, the intensity of the emission is determined by the spatial template and therefore it depends on the position—in other words it is really the spectral shape that is bound to be constant across the source extension;
there are ways to break an extended source in an arbitrary number of smaller components, each with its own independent spectral shape—which effectively allows to overcome the limitation, albeit in a possibly inefficient way;
the Chandra-to-IXPE conversion tool shipped with ixpeobssim provides an alternative observation-simluation strategy fully preserving the correlations between the sky position and the spectrak shape.

Say that you want, e.g., simulate a stationary source with a power-law energy spectrum; all you have to do is something along the lines of the following snippet.

def spec(E, t):
    """Definition of the photon spectrum (cm^-2 s^-1 keV^-1).

    Change the power-law normalization and index to the values that suit
    you. Also, note that t is not used, i.e., the source is stationary, but
    you are free to put in the function body whatever complicated expression
    of energy `and` time, e.g., a prefactor and an index that vary with time.
    """
    return 10. * (E**-2.)

On a related note, a polarization degree constant in energy, time and position in the sky is readily specified.

def pol_deg(E, t, ra, dec):
    """Definition of the polarization degree (between 0 and 1).
    """
    return 0.5

Note that in all cases the signature of the function must be the correct one, no matter whether any specific dynamical variable is effectively used in the function body or not.

Tip

Since we don’t necessarily want you to reinvent the wheel every time, ixpeobssim provides convenience wrappers to power laws and constants, among other things. The first and last snippets above would be more succintly coded as

from ixpeobssim.srcmodel.spectrum import power_law
from ixpeobssim.srcmodel.polarization import constant

spec = power_law(10., 2.)
pol_deg = constant(0.5)

And you should take a look at the documentation of the two functions for more details (note that the power-law wrapper support time-varying prefactor and index).

Using XSPEC models¶

Arbitrary XSPEC models can be fed into xpobssim, providing that PyXSPEC is up and running on your system. (If you are in a hurry, the section about Toy source models provides a working example that you can take inspiration from.)

In order to effectively use an underlying XSPEC model, you have to create a ixpeobssim.srcmodel.spectrum.xXspecModel instance. This is essentially an univariate, interpolated spline (coming with all the standard spline facilities), with the only exception that it can be called with a phony, second argument (in addition to the energy) so that it can be readily used to define a source component—in this case the second argument is acting as the time and has no effect.

The easiest way to create an XSPEC model is through a generic expression and a list specifying all the model parameters.

from ixpeobssim.srcmodel.spectrum import xXspecModel

spec = xXspecModel('powerlaw', [2., 0.1])

The ixpeobssim.srcmodel.spectrum.xXspecModel constructor mimics the underlying XSPEC bindings and allows passing the parameters in the form of a dictionary indexed by sequential identifiers, which is handy when one wants to only set a subset of them and use the XSPEC defaults for the others:

from ixpeobssim.srcmodel.spectrum import xXspecModel

spec = xXspecModel('powerlaw', {1: 2., 2: 0.1})
spec = xXspecModel('powerlaw', {2: 0.1})

Alternatively, the model can be read from a text file with a suitable syntax, see ixpeobssim.srcmodel.spectrum.xXspecModel.from_file()

from ixpeobssim.srcmodel.spectrum import xXspecModel

spec = xXspecModel.from_file('path/to/your/file')

In addition, XSPEC models provide a convenient facility to dump the data points to an ASCII file, see ixpeobssim.srcmodel.spectrum.xXspecModel.save_ascii().

Model components¶

The concept of a model component is central in the definition of an observation simulation setup, and is used throughout ixpeobssim to encapsulate the relevant properties of a few slightly different, but related, things:

a simple celestial source that cannot be further decomposed, such as that defined in [github]/ixpeobssim/config/toy_point_source.py;
different physical components of the same source—e.g., the thermal and non thermal emission in Cas A as modeled in [github]/ixpeobssim/config/toy_casa.py;
different sources in the same ROI—e.g., the Crab pulsar and nebula in [github]/ixpeobssim/config/crab_complex.py.

From a technical standpoint, a model component is basically an instance of the ixpeobssim.srcmodel.roi.xCelestialModelComponentBase class (or, to be more precise, an instance of a subclass of the base class) and encapsulates, among other things, the source name, photon spectrum, and polarization degree and angle, the three latter quantities being specified at construction time as explained in the previous section. (Take a look at the signature of the constructors in the module documentation for more details.)

The subclasses of ixpeobssim.srcmodel.roi.xModelComponentBase deal for the most part (with the exception of that describing periodic sources) with the source morphology, re-implementing the method for extracting vectors of positions in the sky (ra, dec). For completeness, the actual subclasses that can be instantiated and used in simulations are

Source type	Specialized subclass
Point source	`ixpeobssim.srcmodel.roi.xPointSource`
Periodic source	`ixpeobssim.srcmodel.roi.xPeriodicPointSource`
Binary system	`ixpeobssim.srcmodel.roi.xBinarySource`
Uniform disk	`ixpeobssim.srcmodel.roi.xUniformDisk`
Gaussian disk	`ixpeobssim.srcmodel.roi.xGaussianDisk`
Uniform annulus	`ixpeobssim.srcmodel.roi.xUniformAnnulus`
Extended source	`ixpeobssim.srcmodel.roi.xExtendedSource`
Chandra observation	`ixpeobssim.srcmodel.roi.xChandraObservation`
Instrumental background	`ixpeobssim.srcmodel.bkg.xInstrumentalBkg`
Instrumental background	`ixpeobssim.srcmodel.bkg.xPowerLawInstrumentalBkg`
Extragalactic background	`ixpeobssim.srcmodel.bkg.xExtragalacticBkg`
Galactic background	`ixpeobssim.srcmodel.bkg.xGalacticBkg`

and provide a fairly broad and diverse selection of morphological characteristics that should be sufficient for most of the applications.

Now, the simplest possible model component (which in this case is a fully-fledged source) can be defined with only a few lines of code.

import numpy

from ixpeobssim.srcmodel.roi import xPointSource
from ixpeobssim.srcmodel.spectrum import power_law
from ixpeobssim.srcmodel.polarization import constant

ra, dec = 45., 45.
spec = power_law(10., 2.)
pol_deg = constant(0.5)
pol_ang = constant(numpy.radians(65.))
src = xPointSource('Point source', ra, dec, spec, pol_deg, pol_ang)

print(src)

And, for completeness, the output of the last print function should be something along the lines of the following snippet. We’re almost there.

...
xPointSource "Point source" (id = None)
    Galactic column density: 0.000e+00 cm^{-2}
    Redshift: 0.000
    Unabsorbed flux @ t = 0: 2.221e-08 erg/cm2/s (1110.55 mcrab)
    Position: RA = 45.0 deg, Dec = 45.0 deg

Note

ixpeobssim has a concept of Calibration sources and associated region of interest. The reader is referred to the toy_flat_field.py and fcw_calC.py configuration files for usage examples.

(Note that calibration configuration files are supposed to be run through xpcalib rather than xpobssim.)

Regions of interest¶

A complete model of a region of interest is basically a collection of an arbitrary number of model components—and it is the basic simulation unit in the context of ixpeobssim.

From an implementation standpoint, it is an instance of the ixpeobssim.srcmodel.roi.xROIModel class, and is specified passing the coordinates of its center; components are added after the fact. A small tweak to the last snippet in the previous section

import numpy

from ixpeobssim.srcmodel.roi import xPointSource, xROIModel
from ixpeobssim.srcmodel.spectrum import power_law
from ixpeobssim.srcmodel.polarization import constant

# Define the source properties and create the actual source.
ra, dec = 45., 45.
spec = power_law(10., 2.)
pol_deg = constant(0.5)
pol_ang = constant(numpy.radians(65.))
src = xPointSource('Point source', ra, dec, spec, pol_deg, pol_ang)

# Create the complete region of interest.
ROI_MODEL = xROIModel(ra, dec)
ROI_MODEL.add_source(src)

print(ROI_MODEL)

make it for a valid ixpeobsim configuration file that can be used to run an actual simulation. The print output will look like this.

...
ROI centered at (45.0000, 45.0000):
- xPointSource "Point source" (id = 0)
    Galactic column density: 0.000e+00 cm^{-2}
    Redshift: 0.000
    Unabsorbed flux @ t = 0: 2.221e-08 erg/cm2/s (1110.55 mcrab)
    Position: RA = 45.0 deg, Dec = 45.0 deg

Chandra region of interest¶

The region of interest for a Chandra observation follows the same concept as described above, with the difference that the spectral and morphological information of the source is directly taken from the provided Chandra event list and does not need to be specified. Then the xpobssim tool will take care to convert the Chandra observation into a correponding IXPE event list, opportunely down-sampling and smearing the events with the instrument response functions. This simulation technique has the advantage to preserve the full correlation between the source morphology and the spectrum, especially important for extended sources.

From tecnical standpoint, it is an instance of the ixpeobssim.srcmodel.roi.xChandraROIModel class, derived from the ixpeobssim.srcmodel.roi.xROIModel base class, and is initialized by passing the Chandra event file and choosing the ACIS detector used for the observation (I or S). Sub-regions of the ROI can be defined (with individual polarization models) as instances of the ixpeobssim.srcmodel.roi.xChandraObservation class by specifying the id, polarization degree and angle and the region definition (via a ds9 region file).

import os

from ixpeobssim.srcmodel.roi import xChandraObservation, xChandraROIModel
from ixpeobssim.srcmodel.polarization import constant
from ixpeobssim.utils.astro import read_ds9
from ixpeobssim import IXPEOBSSIM_CONFIG

# This is the input event file taken from the Chandra database (obs id 8489)
# http://cda.harvard.edu/pop/mainEntry.do
# ftp://cdaftp.cfa.harvard.edu//pub/science/ao08/cat7/8489/primary/acisf08489N003_evt2.fits.gz
EVENT_FILE_PATH = os.path.join(IXPEOBSSIM_CONFIG, 'fits', 'cena_evt.fits')

# cena_jet+core contains the definition of jet and core in the Chandra map
REG_SOURCE_FILE_PATH = os.path.join(IXPEOBSSIM_CONFIG, 'fits',
                                    'cena_jet+core.reg')

regs = read_ds9(REG_SOURCE_FILE_PATH)

ROI_MODEL = xChandraROIModel(EVENT_FILE_PATH, acis='I')

polarization_degree = constant(0.)
polarization_angle = constant(0.)

# jet region
jet = xChandraObservation('Jet', polarization_degree, polarization_angle, regs[0])
ROI_MODEL.add_source(jet)
# core region
core = xChandraObservation('Core', polarization_degree, polarization_angle, regs[1])
ROI_MODEL.add_source(core)
# the remaining part
core = xChandraObservation('Others', polarization_degree, polarization_angle)
ROI_MODEL.add_source(core)

It is also possible to remove sources in the region of interest by passing the exclude option in the instance of the ixpeobssim.srcmodel.roi.xChandraObservation class and to add sources using the standard model components described above. These features are particularly useful for treating the pileup phenomenon occurring in Chandra CCD detectors for bright sources, whenever two or more photons are detected as a single event. In these cases, the simplest strategy is to remove the bright source suffering of pileup from the conversion ROI and replace it with a standard model component with an appropriate source model description. The converter will then take care of it, simulating the source in the standard way and merging the photon lists at the end of the process. For reference, here are the additional lines that need to be included to the source model example above to remove a region and add a point source to the chandra region of interest:

# core region, exclude the core region by passing the exclude=True
core = xChandraObservation('Core', polarization_degree, polarization_angle,
                           regions[1], exclude=True)
ROI_MODEL.add_source(core)
# the remaining part
core = xChandraObservation('Others', polarization_degree, polarization_angle)
ROI_MODEL.add_source(core)

#Add a point source at a given location in the roi.
RA = 201.43
DEC = -43.00
PL_NORM = 0.0004,
PL_INDEX = 2.
spec = power_law(PL_NORM, PL_INDEX)
src = xPointSource('Point source', RA, DEC, spec, polarization_degree, 
                   polarization_angle)
ROI_MODEL.add_source(src)

Source model etiquette¶

Although the ROI_MODEL top-level object is all ixpeobssim is looking for in a configuration file (and everything will run happily as long as you provide that in a valid form) there is undoubtedly value in trying and keep a consistent style when creating source models. People will get acquainted to the conventions and will understand new models more easily. This section is essentially a list of suggestions on how to go about creating your own models.

We note, in passing, that the section about Toy source models is a good read, along with the material in this section.

Give your model file a sensible name (which, in most cases, should be obvious). It is conceivable that we might have multiple, alternative, configuration files for a given source; in this case, use your best judgement to try and convey with the file name the intent of any of these alternative models. (Having crab.py, crab2.py, crab_new.py, crab_newer.py stops being funny after a while.) Even more important, try and factor out as much as possible of the code that is shared across multiple models—resist the temptation to cut and paste code around!

Try and document your model in a reasonable fashion. This is typically achieved by putting arbitrary reStructuredText that can be then parsed in sphinx at the beginning of the Python configuration module, right after the license notice and before the __future__ Python imports.

Warning

For the documentation to be parsed by sphinx, the reStructuredText has to be before any actual Python code. Keep this in mind.

When you run a simulation and you want to compare the ixpeobssim output with the input model, you will need to have all the input parameters at hand. Keep this in mind when writing configuration files:

spec = power_law(1., 2.)

is more succinct than

pl_norm = 1.
pl_index = 2.
spec = power_law(pl_norm, pl_index)

but the second formulation makes it easier for you to retrieve the input spectral index when you fit your binned PHA1 spectrum in xspec and want to compare the best-fit parameters with the models.

Take a second to read our Coding guidelines and try and stick to them. In the long run people will be grateful to you for doing it.

Create a display() method in each configuration file that, when called, will plot anything that might be interesting in you model. (Each source model is different and it is hard to achieve this through a common, completely automated interface.) This method will not partecipate in the observation simulation when the model is parsed by ixpeobssim, but is something that, when protected by a if __name__ == '__main__' can be useful as a quick-look interface and for creating figures for the documentation. To this end, we have implemented a small bootstrapping function that you are welcome to use by doing something along the lines of

if __name__ == '__main__':
    from ixpeobssim.config import bootstrap_display
    bootstrap_display()

This will create a custom option parse for you with a limited set of options (e.g., to save the plots in a specific folder) and will fire up your display() method right away. Try at any time

python path/to/config/file.py --help

to have an up-to date help output.

Being consistent in terms of axis labels and units, conventions for file names and all that is cool. Really. Take a seconds to acquaint yourself with the following small utility functions and classes:

ixpeobssim.utils.matplotlib_.save_gcf(): save the current figure
ixpeobssim.utils.matplotlib_.setup_gca(): setup the current plot
ixpeobssim.utils.fmtaxis.fmtaxis: list of axis labels and units

Creating source models¶

Spectra and polarization patterns¶

Using XSPEC models¶

Model components¶

Regions of interest¶

Chandra region of interest¶

Source model etiquette¶

Navigation

Related Topics