API

Top-level indices

• Index

• Module Index

• Search Page

Core Modules

Facilities related to FITS I/O.

ixpeobssim.core.fitsio.copy_hdu_list(hdu_list)

Create a copy in memory of a HDU list.

Note that the HDUList class provides a copy() method, but this is a shallow copy that is still tied to the underlying file the HDU list is read from: https://github.com/astropy/astropy/issues/6395

As a workaround we copy the HDUs one by one and assemble them back into a HDUList object.

ixpeobssim.core.fitsio.find_column_index(hdu, col_name)

Return the index corresponding to a given column name within a binary table, and None when the column does not exists.

Parameters:

• hdu (fits.BinTableHDU instance) – The underlying HDU

• col_name (str) – The name of the target column.

ixpeobssim.core.fitsio.read_hdu_list_in_memory(file_path, extension='fits')

Read a fits file and store the HDU list in memory.

This is a convenience function that is meant to open a file within a context manager (so that it is properly closed after the fact) and create a copy of the content in memory. This is less than trivial, as astropy typically keeps some kind of pointer back to the original file when dealing with FITS files, and it is not obvious how to make sure that all the resources are properly closed at runtime.

ixpeobssim.core.fitsio.set_tlbounds(hdu, col_name, min_, max_)

Set the proper TLMIN and TLMAX header keywords for a given column in a given HDU.

class ixpeobssim.core.fitsio.xBinTableHDUBase(data=None, keywords=None, comments=None)

Binary table HDU class.

This is a small wrapper around a standard binary table to facilitate customizations.

set_ext_name(name)

Set the extension name for the binary table.

classmethod spec_names()

Return the name of the data fields specified in the SPEC class member.

classmethod spec_names_and_types()

Return the name of the data fields specified in the SPEC class member.

class ixpeobssim.core.fitsio.xFITSImageBase(file_path)

Base class describing a FITS image.

Parameters:

file_path (string) – The path to the FITS file containing the image.

add_circle(x, y, radius, mode='radec', **kwargs)

Add a circle to the image.

Note that the x and y coordinates can be either RA and Dec in decimal degrees (if mode=’radec’) or relative to canvas coordinates (if mode=’canvas’).

Parameters:

• x (float) – The x position of the circle center (Right Ascension in decimal degrees if relative=False, relative position in canvas coordinates if relative=True).

• y (float) – The y position of the circle center (Declination in decimal degrees if relative=False, relative position in canvas coordinates if relative=True).

• radius (float) – The radius of the circle in arcseconds.

• mode ({'radec', 'canvas'}) – Optional flag defining how the center of the circle is expresses: in right ascension and declination (‘radec’) or in relative canvas coordinates (‘canvas’).

• kwargs (dict) – All the keyword arguments passed to the matplotlib Circle patch.

add_circle_canvas(ra, dec, radius, **kwargs)

Convenience function, see add_circle().

add_circle_radec(ra, dec, radius, **kwargs)

Convenience function, see add_circle().

static add_label(text, x=0.1, y=0.9, **kwargs)

Add a label to an image.

This is a shortcut to have all the formatting defined.

Parameters:

• text (string) – The label text

• x (float) – The x position of the label in relative coordinates

• y (float) – The y position of the label in relative coordinates

• **kwargs – All the keyword arguments passed to plt.text()

cdelt1()

Return the value of the CRVAL1 header keyword.

cdelt2()

Return the value of the CRVAL2 header keyword.

center()

Return the sky coordinates of the image center.

crpix1()

Return the value of the CRPIX1 header keyword.

crpix2()

Return the value of the CRPIX2 header keyword.

crval1()

Return the value of the CRVAL1 header keyword.

crval2()

Return the value of the CRVAL2 header keyword.

delta()

Return the pixel increment in the ra and dec coorinates.

get(keyword, default=None)

Retrieve the value of a primary header keyword.

inner_radius()

Return the inner radius of the image in degrees.

static make_plot(data, wcs_, slices=None, zlabel=None, stretch='linear', ticks=None, colorbar_format=None, **kwargs)

Underlying plotting routine.

This is implemented as a staticmethod in such a way it can be called from external methods dealing with multi-dimensional data arrays, which are not readily supported in the base class. This is achieved by passing the slices argument, aling with the proper data slice.

Parameters:

• data (array_like) – The underlying data to be plotted.

• wcs (astropy.wcs object) – The WCS object defining the transformation from pixel to sky coordinates.

• slices (tuple, optional) – Verbatim from the astropy.visualization documentation: For WCS transformations with more than two dimensions, we need to choose which dimensions are being shown in the 2D image. The slice should contain one x entry, one y entry, and the rest of the values should be integers indicating the slice through the data. The order of the items in the slice should be the same as the order of the dimensions in the WCS, and the opposite of the order of the dimensions in Numpy. For example, (50, ‘x’, ‘y’) means that the first WCS dimension (last Numpy dimension) will be sliced at an index of 50, the second WCS and Numpy dimension will be shown on the x axis, and the final WCS dimension (first Numpy dimension) will be shown on the y-axis (and therefore the data will be plotted using data[:, :, 50].transpose())

• zlabel (str) – The label for the z axis (e.g., the colorbar).

• stretch ({‘linear’, ‘sqrt’, ‘power’, log’, ‘asinh’}, optional) – The stretch function to apply to the image.

• ticks (array_like or ticker, optional) – The explicit setting for the colorbar ticks.

• colorbar_format ({None, 'log', 'scilog'}, optional) – Optional string for the colorbar formnatting.

• kwargs (dict) – All the keyword arguments passed to plt.imshow().

plot(zlabel=None, stretch='linear', ticks=None, colorbar_format=None, **kwargs)

Plot the underlying FITS image.

Since ixpeobssim version 8.9.0 this not using aplpy under the hood anymore, see https://bitbucket.org/ixpesw/ixpeobssim/issues/272

Parameters:

• zlabel (string) – The text label for the colorbar (use None for no colorbar)

• **kwargs – All the keyword arguments passed to plt.imshow()

plot_arrows(nside, model, threshold=0.0, **kwargs)

Overlay arrows on an image.

recenter(ra, dec, radius)

Recenter the image.

Parameters:

• ra (float) – The RA position of the circle center in decimal degrees.

• dec (float) – The DEC position of the circle center in decimal degrees.

• radius (float) – The radius of the region in arcseconds.

shape()

Return the shape of the image.

sky_bounding_box()

Return the bounding box of the FITS image in sky coordinates, i.e., a 4-element tuple of the form (ra_min, dec_min, ra_max, dec_max).

square_sky_grid(nside, ra0=None, dec0=None, half_size=None)

Create a square, regular grid in ra and dec.

This is essentially an overloaded version of the utils.astro.square_sky_grid() method, where the arguments are inferred, when possible, from the underlying image structure.

class ixpeobssim.core.fitsio.xHDUBase

Base class for FITS HDU.

add_comment(comment)

Add a comment to the table header.

add_keyword(key, value, comment='')

Add a keyword to the table header.

set_keyword(key, value)

Set a keyword into the table header.

set_keyword_comment(keyword, comment)

Set the comment for a header keyword.

setup_header(keywords=None, comments=None)

Update the table header with arbitrary additional information.

class ixpeobssim.core.fitsio.xPrimaryHDU(data=None, header=None, creator='ixpeobssim', keywords=None, comments=None)

Class describing a primary HDU to be written in a FITS file.

This is initializing a standard astropy.io.fits.PrimaryHDU object and adding the creator and creation time fields.

Parameters:

creator (str) – The application that created the header (defaults to ixpeobssim, and the ixpeobssim tag is automatically added.)

ixpeobssim.core.fitting.fit(model, xdata, ydata, p0=None, sigma=None, xmin=-inf, xmax=inf, absolute_sigma=True, check_finite=True, method=None, verbose=True, **kwargs)

Lightweight wrapper over the scipy.optimize.curve_fit() function to take advantage of the modeling facilities. More specifically, in addition to performing the actual fit, we update all the model parameters so that, after the fact, we do have a complete picture of the fit outcome.

Parameters:

• model (ixpeobssim.core.modeling.FitModelBase instance callable) – The model function, f(x, …). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments.

• xdata (array_like) – The independent variable where the data is measured.

• ydata (array_like) – The dependent data — nominally f(xdata, …)

• p0 (None, scalar, or sequence, optional) – Initial guess for the parameters. If None, then the initial values will all be 1.

• sigma (None or array_like, optional) – Uncertainties in ydata. If None, all the uncertainties are set to 1 and the fit becomes effectively unweighted.

• xmin (float) – The minimum value for the input x-values.

• xmax (float) – The maximum value for the input x-values.

• absolute_sigma (bool, optional) – If True, sigma is used in an absolute sense and the estimated parameter covariance pcov reflects these absolute values. If False, only the relative magnitudes of the sigma values matter. The returned parameter covariance matrix pcov is based on scaling sigma by a constant factor. This constant is set by demanding that the reduced chisq for the optimal parameters popt when using the scaled sigma equals unity.

• method ({'lm', 'trf', 'dogbox'}, optional) – Method to use for optimization. See least_squares for more details. Default is ‘lm’ for unconstrained problems and ‘trf’ if bounds are provided. The method ‘lm’ won’t work when the number of observations is less than the number of variables, use ‘trf’ or ‘dogbox’ in this case.

• verbose (bool) – Print the model if True.

• kwargs – Keyword arguments passed to leastsq for method='lm' or least_squares otherwise.

ixpeobssim.core.fitting.fit_gaussian_iterative(histogram, p0=None, sigma=None, xmin=-inf, xmax=inf, absolute_sigma=True, check_finite=True, method=None, verbose=True, num_sigma_left=2.0, num_sigma_right=2.0, num_iterations=2, **kwargs)

Fit the core of a gaussian histogram within a given number of sigma around the peak.

This function performs a first round of fit to the data and then repeats the fit iteratively limiting the fit range to a specified interval defined in terms of deviations (in sigma) around the peak.

For additional parameters look at the documentation of the ixpeobssim.core.fitting.fit_histogram()

Parameters:

• num_sigma_left (float) – The number of sigma on the left of the peak to be used to define the fitting range.

• num_sigma_right (float) – The number of sigma on the right of the peak to be used to define the fitting range.

• num_iterations (int) – The number of iterations of the fit.

ixpeobssim.core.fitting.fit_histogram(model, histogram, p0=None, sigma=None, xmin=-inf, xmax=inf, absolute_sigma=True, check_finite=True, method=None, verbose=True, **kwargs)

Fit a histogram to a given model.

This is basically calling ixpeobssim.core.fitting.fit() with some pre-processing to turn the histogram bin edges and content into x-y data. Particularly, the bin centers are taken as the independent data series, the bin contents are taken as the dependent data saries, and the square root of the counts as the Poisson error.

For additional parameters look at the documentation of the ixpeobssim.core.fitting.fit()

Parameters:

• model (ixpeobssim.core.modeling.FitModelBase instance or) – callable The fit model.

• histogram (ixpeHistogram1d instance) – The histogram to be fitted.

Warning

We’re not quite doing the right thing, here, as we should integrate the model within each histogram bin and compare that to the counts, but this is not an unreasonable first-order approximation. We might want to revise this, especially since we can probably provide an analytic integral for most of the model we need.

ixpeobssim.core.fitting.fit_modulation_curve(histogram, p0=None, sigma=None, xmin=-inf, xmax=inf, absolute_sigma=True, check_finite=True, method=None, verbose=True, degrees=True, **kwargs)

Fit a modulation curve.

For all the other parameters look at the documentation of the ixpeobssim.core.fitting.fit_histogram()

ixpeobssim.core.fitting.linear_analytical_fit(x, y, dy=None)

Simple vectorized, analytical linear fit.

See https://mathworld.wolfram.com/LeastSquaresFitting.html

ixpeobssim.core.fitting.power_law_analytical_fit(x, y)

Simple vectorized, analytical linear fit.

See https://mathworld.wolfram.com/LeastSquaresFitting.html

Geometry module.

class ixpeobssim.core.geometry.xLine(begin, end)

Class representing a line.

length()

Return the length of the line.

class ixpeobssim.core.geometry.xPoint(x, y, z)

Class representing a point in the three-dimensional space.

distance_to(other)

Return the distance to another point.

move(length, ray)

Move the point by a given length along a given ray.

norm()

Return the norm of the three-dimensional vector corresponding to the point.

unphysical()

Return True if the point is unphysical.

classmethod unphysical_point()

Return an unphysical point.

class ixpeobssim.core.geometry.xRay(origin, theta, phi)

Class representing a ray in the three-dimensional space.

xintersect(x)

Return the point where the ray intersects the plane at a given x.

yintersect(y)

Return the point where the ray intersects the plane at a given y.

zintersect(z)

Return the point where the ray intersects the plane at a given z.

Histogram facilities.

class ixpeobssim.core.hist.xGpdMap2d(nside=10, zlabel='Entries/bin')

2-dimensional GPD map.

class ixpeobssim.core.hist.xGpdMap3d(nside, zbins, zlabel='', wlabel='Entries/bin')

Three-dimensional histogram where the first two axes represent the GPD active area in Physical coordinates.

class ixpeobssim.core.hist.xHistogram1d(xbins, xlabel='', ylabel='Entries/bin')

Container class for one-dimensional histograms.

errorbar(**kwargs)

Plot the histogram as a scatter plot.

errorbar_data()

Return the x, y, dy arrays that can be used to build a scatter plot (with errors) from the histogram.

fit(model, p0=None, sigma=None, xmin=-inf, xmax=inf, absolute_sigma=True, check_finite=True, method=None, verbose=True, **kwargs)

Fit the histogram with a model.

gaussian_kde_smooth(bandwidth=2.0)

Create a copy of the histogram where the weights are smoothed with a gaussian kernel density estimatore.

Parameters:

bandwidth (float) – The sigma of the gaussian kernel, in (fractional) number of bins.

scatter_plot()

Turn the histogram into a scatter plot.

Warning

This is to be removed.

class ixpeobssim.core.hist.xHistogram2d(xbins, ybins, xlabel='', ylabel='', zlabel='Entries/bin')

Container class for two-dimensional histograms.

gaussian_kde_smooth(bandwidth=(2.0, 2.0))

Create a copy of the histogram where the weights are smoothed with a gaussian kernel density estimatore.

Warning

This is essentially untested.

Parameters:

bandwidth (2-element tuple float) – The sigma of the gaussian kernel, in (fractional) number of bins.

hbisect(y)

Return the horizontal slice corresponding to a given y value.

hslice(bin_)

Return the horizontal slice for a given bin.

hslices()

Return a list of all the horizontal slices.

vbisect(x)

Return the vertical slice corresponding to a given y value.

vslice(bin_)

Return the vertical slice for a given bin.

vslices()

Return a list of all the vertical slices.

class ixpeobssim.core.hist.xHistogram3d(xbins, ybins, zbins, xlabel='', ylabel='', zlabel='', wlabel='Entries/bin')

Container class for three-dimensional histograms.

class ixpeobssim.core.hist.xHistogramBase(binning, labels)

Base class for an n-dimensional histogram.

This interface to histograms is profoundly different for the minimal numpy/matplotlib approach, where histogramming methods return bare vectors of bin edges and counts. The main underlying ideas are

• we keep track of the bin contents, the bin entries and the sum of the weights squared (the latter for the purpose of calculating the errors);

• we control the axis label and the plotting styles;

• we provide two separate interfaces, fill() and set_content(), to fill the histogram from either unbinned or binned data;

• we support the basic arithmetics (addition, subtraction and multiplication by a scalar);

• we support full data persistence (I/O) in FITS format.

Note that this base class is not meant to be instantiated directly, and the interfaces to concrete histograms of specific dimensionality are defined in the sub-classes.

Parameters:

• binning (n-tuple of array) – the bin edges on the different axes.

• labels (n-tuple of strings) – the text labels for the different axes.

bin_centers(axis=0)

Return the bin centers for a specific axis.

bin_widths(axis=0)

Return the bin widths for a specific axis.

static binning_col_name()

Column name for the axis binnings.

static binning_hdu_name(axis)

Extension name for the axis binnings.

static bisect(binning, values, side='left')

Return the indices corresponding to a given array of values for a given binning.

copy()

Create a full copy of a histogram.

empty_copy()

Create an empty copy of a histogram.

static entries_hdu_name()

Extension name for the bin entries.

errors()

Return the bin errors.

fill(*data, weights=None)

Fill the histogram from unbinned data.

Note this method is returning the histogram instance, so that the function call can be chained.

find_bin(*coords)

Find the bin corresponding to a given set of “physical” coordinates on the histogram axes.

This returns a tuple of integer indices that can be used to address the histogram content.

find_bin_value(*coords)

Find the histogram content corresponding to a given set of “physical” coordinates on the histogram axes.

classmethod from_file(file_path)

Load the histogram from a FITS file.

Note that we transpose the image data back at read time—see the comment in the save() method.

static label_keyword(axis)

Header keyword for the axis labels.

mean(axis=0)

Calculate the (binned) mean eastimate along a given axis.

num_entries()

Return the number of entries in the histogram.

plot(**kwargs)

Plot the histogram.

rms(axis=0)

Calculate the (binned) rms eastimate along a given axis.

Warning

Mind this is not using anything clever (e.g. the Welford algorithm) and is not particularly numerically stable.

save(file_path, overwrite=True, **header_keywords)

Save the histogram (with edges) to a FITS file.

Note that all the image data are transposed so that the thing can be correctly visualized with standard FITS viewers—at least in two dimensions.

set_axis_label(axis, label)

Set the label for a given axis.

set_content(content, entries=None, errors=None)

Set the bin contents programmatically from binned data.

Note this method is returning the histogram instance, so that the function call can be chained.

set_errors(errors)

sum(axis=None)

return the sum of weights in the histogram.

static sumw2_hdu_name()

Extension name for the bin entries.

class ixpeobssim.core.hist.xModulationCube2d(xbins, ybins)

Specialized class for modulations cubes.

class ixpeobssim.core.hist.xScatterPlot(x, y, dy=None, dx=None, xlabel=None, ylabel=None)

Small class encapsulating a scatter plot.

Technically speaking, this would not belong here, as a scatter plot is not strictly related to any of the histogram classes, but a 1-dimensional histogram with errors can techincally be turned into a scatter plot, so we introduce the concept here for completeness.

Warning

Consider removing this class.

plot(**kwargs)

Plot the scatter plot.

class ixpeobssim.core.modeling.xConstant

Constant model.

cdf(x)

Overloaded cdf() method.

init_parameters(xdata, ydata, sigma)

Overloaded init_parameters() method.

static jacobian(x, constant)

Overloaded jacobian() method.

ppf(q)

Overloaded ppf() method.

static value(x, constant)

Overloaded value() method.

class ixpeobssim.core.modeling.xExponential

Exponential model.

static jacobian(x, normalization, exponent)

Overloaded jacobian() method.

static value(x, normalization, exponent)

Overloaded value() method.

class ixpeobssim.core.modeling.xExponentialOffset

Exponential model.

static jacobian(x, offset, normalization, exponent)

Overloaded jacobian() method.

static value(x, offset, normalization, exponent)

Overloaded value() method.

class ixpeobssim.core.modeling.xFe55

One-dimensional double gaussian model for Fe55 with 2 lines

init_parameters(xdata, ydata, sigma)

static jacobian(x, amplitude0, amplitude1, peak, sigma)

Overloaded jacobian() method.

plot(*parameters, **kwargs)

Overloaded plot() method.

static value(x, amplitude0, amplitude1, peak, sigma)

Overloaded value() method.

class ixpeobssim.core.modeling.xFitModelBase

Base class for a fittable model.

This base class isn’t really doing anything useful, the idea being that actual models that can be instantiated subclass xFitModelBase overloading the relevant class members.

The class features a number of static members that derived class should redefine as needed:

• PARAMETER_NAMES is a list containing the names of the model parameters. (It goes without saying thet its length should match the number of parameters in the model.)

• PARAMETER_DEFAULT_VALUES is a list containing the default values of the model parameters. When a concrete model object is instantiated these are the values being attached to it at creation time.

• PARAMETER_DEFAULT_BOUNDS is a tuple containing the default values of the parameter bounds to be used for the fitting. The values in the tuple are attached to each model object at creation time and are intended to be passed as the bounds argument of the scipy.optimize.curve_fit() function. From the scipy documentation: Lower and upper bounds on independent variables. Defaults to no bounds. Each element of the tuple must be either an array with the length equal to the number of parameters, or a scalar (in which case the bound is taken to be the same for all parameters.) Use np.inf with an appropriate sign to disable bounds on all or some parameters. By default models have no built-in bounds.

• DEFAULT_PLOTTING_RANGE is a two-element list with the default support (x-axis range) for the model. This is automatically updated at runtime depending on the input data when the model is used in a fit.

• DEFAULT_STAT_BOX_POSITION is the default location of the stat box for the model, see the gpdswpy.matplotlib_.stat_box() function for all the details.

In addition, each derived class should override the following things:

• the value(x, *args) static method: this should return the value of the model at a given point for a given set of values of the underlying parameters;

• (optionally) the jacobian(x, *args) static method. (If defined, this is passed to the underlying fit engine allowing to reduce the number of function calls in the fit; otherwise the jacobian is calculated numerically.)

Finally, if there is a sensible way to initialize the model parameters based on a set of input data, derived classes should overload the init_parameters(xdata, ydata, sigma) method of the base class, as the latter is called by fitting routines if no explicit array of initial values are passed as an argument. The default behavior of the class method defined in the base class is to do nothing.

See gpdswpy.modeling.xGaussian for a working example.

init_parameters(xdata, ydata, sigma)

Assign a sensible set of values to the model parameters, based on a data set to be fitted.

Note that in the base class the method is not doing anything, but it can be reimplemented in derived classes to help make sure the fit converges without too much manual intervention.

integral(edges)

Calculate the integral of the model within pre-defined edges.

Note that this assumes that the derived class provides a suitable cdf() method.

name()

Return the model name.

parameter_error(name)

Return the parameter error by name.

parameter_errors()

Return the vector of parameter errors.

parameter_status()

Return the complete status of the model in the form of a tuple of tuples (parameter_name, parameter_value, parameter_error).

Note this can be overloaded by derived classes if more information needs to be added.

parameter_value(name)

Return the parameter value by name.

parameter_values()

Return the vector of parameter values.

plot(*parameters, **kwargs)

Plot the model.

Note that when this is called with a full set of parameters, the self.parameters class member is overwritten so that the right values can then be picked up if the stat box is plotted.

reduced_chisquare()

Return the reduced chisquare.

reset()

Reset all the necessary stuff.

This method initializes all the things that are necessry to keep track of a parametric fit.

• the parameter values are set to what it specified in PARAMETER_DEFAULT_VALUES;

• the covatiance matrix is initialized to a matrix of the proper dimension filled with zeroes;

• the minimum and maximum values of the independent variable (relevant for plotting) are set to the values specified in DEFAULT_PLOTTING_RANGE;

• the model bounds are set to the values specified in PARAMETER_DEFAULT_BOUNDS;

• the chisquare and number of degrees of freedom are initialized to -1.

rvs(size=1)

Return random variates from the model.

Note that this assumes that the derived class provides a suitable ppf() method.

set_parameter(name, value)

Set a parameter value.

set_parameters(*pars)

Set all the parameter values.

Note that the arguments must be passed in the right order.

set_plotting_range(xmin, xmax)

Set the plotting range.

stat_box(position=None, plot=True, **kwargs)

Plot a ROOT-style stat box for the model.

static value(x, *parameters)

Eval the model at a given point and a given set of parameter values.

Warning

This needs to be overloaded in any derived class for the thing to do something sensible.

class ixpeobssim.core.modeling.xGaussian

One-dimensional Gaussian model.

static der_amplitude(x, amplitude, peak, sigma)

Return the amplitude derivative of the function, to be used in the calculation of the Jacobian.

static der_peak(x, amplitude, d_amplitude, peak, sigma)

Return the peak derivative of the function, to be used in the calculation of the Jacobian.

Note that we pass the pre-calculated values of the amplitude derivatives in order not to repeat the same calculation more times than strictly necessary.

static der_sigma(x, amplitude, d_amplitude, peak, sigma)

Return the sigma derivative of the function, to be used in the calculation of the Jacobian.

Note that we pass the pre-calculated values of the amplitude derivatives in order not to repeat the same calculation more times than strictly necessary.

fwhm()

Return the absolute FWHM of the model.

init_parameters(xdata, ydata, sigma)

Overloaded init_parameters() method.

static jacobian(x, amplitude, peak, sigma)

Overloaded jacobian() method.

resolution()

Return the resolution of the model, i.e., the FWHM divided by the peak value.

resolution_error()

Return the error on the resolution.

static value(x, amplitude, peak, sigma)

Overloaded value() method.

class ixpeobssim.core.modeling.xGeneralizedGaussian

Generalized gaussian fitting model.

See https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gennorm.html for implementation details.

static value(x, norm, mean, sigma, beta)

Overloaded value() method.

class ixpeobssim.core.modeling.xHat

Fit model representing a flat distribution truncated with an error function at both ends.

The ASYMMETRY class member controls the slope of the central part of the model.

static value(x, norm, x1, sigma1, x2, sigma2)

Overloaded value() method.

class ixpeobssim.core.modeling.xLine

Straight-line model.

static jacobian(x, intercept, slope)

Overloaded jacobian() method.

static value(x, intercept, slope)

Overloaded value() method.

class ixpeobssim.core.modeling.xLogNormal

Log-normal fitting model.

See https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html for more details on the implementation.

Note that our parametrization is done in terms of the mean and sigma of the distribution, in addition to the shape parameter.

static value(x, norm, mu, sigma, shape)

Overloaded value() method.

class ixpeobssim.core.modeling.xLorentzian

Lorentzian function

static jacobian(x, normalization, peak, hwhm)

Overloaded jacobian() method.

static value(x, normalization, peak, hwhm)

Overloaded value() method.

class ixpeobssim.core.modeling.xModulationCurveDeg

Modulation curve model (for fitting azimuthal distributions expressed in degrees).

init_parameters(xdata, ydata, sigma)

Overloaded init_parameters() method.

static jacobian(x, amplitude, modulation, phase)

Overloaded jacobian() method.

In addition to the degree-to-radians conversion of the inputs, here we have to multiply the last column of the jacobian (i.e., the phase derivatives) by pi / 180.

static value(x, amplitude, modulation, phase)

Overloaded value() method.

Here we are essentially calling the modulation curve model expressed in radians converting to radians the input values of the independent variable and the phase parameter.

class ixpeobssim.core.modeling.xModulationCurveRad

Modulation curve model (for fitting azimuthal distributions expressed in radians).

init_parameters(xdata, ydata, sigma)

Overloaded init_parameters() method.

static jacobian(x, amplitude, modulation, phase)

Overloaded jacobian() method.

static value(x, amplitude, modulation, phase)

Overloaded value() method.

class ixpeobssim.core.modeling.xPixelPha

Composite model for fitting the single-pixel pulse-height distributions.

This is essentially the sum of a gaussian centered at zero (and modeling the noise) and a power law with exponential cutoff (modeling the signal).

static jacobian(x, noise_amplitude, noise_sigma, signal_normalization, signal_index, signal_cutoff)

Overloaded jacobian() method.

plot(*parameters, **kwargs)

Overloaded plot() method.

static value(x, noise_amplitude, noise_sigma, signal_normalization, signal_index, signal_cutoff)

Overloaded value() method.

class ixpeobssim.core.modeling.xPowerLaw

Power law model.

static jacobian(x, normalization, index)

Overloaded jacobian() method.

static value(x, normalization, index)

Overloaded value() method.

class ixpeobssim.core.modeling.xPowerLawExpCutoff

Power law with an exponential cutoff.

static jacobian(x, normalization, index, cutoff)

Overloaded jacobian() method.

static value(x, normalization, index, cutoff)

Overloaded value() method.

ixpeobssim.core.pipeline.batch()

Return the current value of the ‘batch’ rc parameter.

ixpeobssim.core.pipeline.bootstrap_pipeline(model_name)

Convenience bootstrap function.

ixpeobssim.core.pipeline.figure(name)

Return a matplotlib figure with the model name prepended to the actual name.

ixpeobssim.core.pipeline.figure_name(name)

Small convience function to enfore a minimum uniformity in the naming scheme for the output file.

This will prepeng the model name to any name the use passes as an argument, and change the latter to all lower case.

ixpeobssim.core.pipeline.file_list(*args, **kwargs)

Create a file list from a series of patterns.

This method is very handy to retrieve the files generated by any given step of a pipeline in order to process them in the following step, and takes care automagically of looping over the three detector units.

The arguments can be either a string or a (str, int) tuple, in which case the string is a label and the integer is a uid attached to it (e.g., the identifier for a particular phase selection for a periodic source). The arguments are concatenated in the file name with an “_” character.

Assuming that we have a pipieline bootstrapped with a model called “mymodel”, the resolution rules will yield: pipeline.file_list() -> ‘$IXPEOBSSIM_DATA/mymodel_du1/2/3.fits’ pipeline.file_list((‘sel’, 1)) -> ‘$IXPEOBSSIM_DATA/mymodel_du1/2/3_sel0001.fits’ pipeline.file_list((‘sel’, 1), ‘cmap’) -> ‘$IXPEOBSSIM_DATA/mymodel_du1/2/3_sel0001_cmap.fits’

The label “pha1*” is peculiar in that, if encountered it is automatically expanded to include all the flavors of the Stokes parameter spectra.

Note that, in order to maintain compatibility with Python 2, we opted for the current function signature over what would have been the most natural way to do things in Python 3, i.e., file_list(*args, folder_path=IXPEOBSSIM_DATA, check_files=True)

ixpeobssim.core.pipeline.fit_ensamble_stokes_spectra(seed, normalized=False, **kwargs)

Fit an ensamble file with XSPEC.

This is a small convenience function that, under the hood, collects the relevant file list, runs xpxspec, and return a list with the best-fit parameters and errors.

ixpeobssim.core.pipeline.generate_ensamble(size=10, start_seed=0, processing_function=<function standard_ensamble_processing>, cleanup=True, **kwargs)

Generate an ensable of simulations to study the a-posteriory statistical properties of the measurement setup.

Parameters:

• size (int) – The number of independent realizations to be generated

• start_seed (int) – The first random seed to be used (will be incremented by one at each step)

• cleanup (bool) – If True, remove the (potentially large) event files after the processing.

• **Kwargs – All the keyword arguments to be passes to xpobssim.

ixpeobssim.core.pipeline.glob_ensamble(seed, *patterns)

Glob an ensamble folder searching for the files matching a particular seed and pattern.

Note this is completely different in spirit wrt file_list(), as the list is not assembled a priori, but globbed directly from the file system (i.e., you are guaranteed that the files in the list do exist).

ixpeobssim.core.pipeline.glob_ensamble_stokes_spectra(seed, normalized=False)

Specialized function to get hold of a consistent set of Stokes spectra for a given seed.

If the normalized argument is False, this returns the PHA1, PHA1Q and PHA1U files for the three detector units, otherwise the PHA1QN and PHA1UN spectra are returned.

ixpeobssim.core.pipeline.model()

Return the current model name.

ixpeobssim.core.pipeline.output_folder()

Return the current value of the ‘outfolder’ rc parameter.

ixpeobssim.core.pipeline.overwrite()

Return the current value of the ‘overwrite’ rc parameter.

ixpeobssim.core.pipeline.param(key, default=None)

Retrive a given global configuration parameter.

ixpeobssim.core.pipeline.params()

Return the underlying rc param dictionary.

ixpeobssim.core.pipeline.post_process_ensamble_pcubes(seed, label='pcube')

Post-process the polarization cubes for a given ensamble run.

ixpeobssim.core.pipeline.reset(model_name, **kwargs)

ixpeobssim.core.pipeline.residual_figure(name)

Create a figure for residual plots.

ixpeobssim.core.pipeline.save()

Return the current value of the ‘save’ rc parameter.

ixpeobssim.core.pipeline.set_gcf_name(name)

Set the canvas title for the current figure, following the same rules of the figure_name() method.

ixpeobssim.core.pipeline.set_model(model_name)

Set the pipeline model name.

By default this is the string that will be used, e.g., to create all the file path concatenations.

ixpeobssim.core.pipeline.setup(**kwargs)

Setup the global pipeline params.

ixpeobssim.core.pipeline.standard_ensamble_processing(file_list, mc=True)

Standard processing routine.

ixpeobssim.core.pipeline.suffix(label=None, index=None)

Create a suffix to be appended to file names when the pipeline tools are created.

ixpeobssim.core.pipeline.target()

Return the function target for a particular pipeline run.

ixpeobssim.core.pipeline.xpancrkey(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpbin(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpbinview(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpchrgmap(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpexposure(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpgrppha(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpmdp(**kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpobssim(**kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpophase(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpphase(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpphotonlist(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xppicorr(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xppimms(**kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpradialprofile(**kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpselect(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpsimfmt(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpsonify(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpstokesalign(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpstokesrandom(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpstokesshuffle(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpstokessmear(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpstripmc(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpvisibility(*args, **kwargs)

App wrapper.

ixpeobssim.core.pipeline.xpxspec(*args, **kwargs)

App wrapper.

Custom random generation classes, mostly based on splines.

class ixpeobssim.core.rand.xUnivariateAuxGenerator(rv, aux, pdf, bbox=None, kx=3, ky=3, xlabel='rv', ylabel='aux', zlabel='pdf')

Univariate random generator where the pdf of the random variable depends on an additional auxiliary variable.

Internally, a proper meshgrid of random and auxiliary variable values is created, and the pdf is evaluated on the grid to construct an interpolated bivariate spline, so that a horizontal slice of the bivariate spline at any given value of the auxiliary variable is the actual pdf of the random variable at that value of the auxiliary variable.

Note that the z function argument can either be an array-like object with shape (rv.size, aux.size), or a function that can be evaluated in the bi-dimensional phase space—see the base class xInterpolatedBivariateSpline for the implementation details.

Parameters:

• rv (array_like) – Input values for the random variable (assumed to be sorted).

• aux (array_like) – Input values for the auxiliary variable (assumed to be sorted).

• pdf (array_like of shape (x.size, y.size) or callable) – Input pdf values, with shape (x.size,y.size).

• bbox (array_like, optional) – The boundary of the approximation domain.

• kx (int) – The spline order on the two axes.

• ky (int) – The spline order on the two axes.

• xlabel (str, optional) – A text label for the random variable.

• ylabel (str, optional) – A text label for the auxiliary variable.

• zlabel (str, optional) – A text label for the pdf values.

Note

A major backward-incompatible change was introduced here in version 3.0.0, in the context of https://bitbucket.org/ixpesw/ixpeobssim/issues/196, and the constructor signature of the class now reads as xUnivariateAuxGenerator(rv, aux, pdf…), instead of xUnivariateAuxGenerator(aux, rv, pdf…) as it used to be historically in all the previous versions of ixpeobssim.

(On a related note, since some of the arguments are passed by name in the base spline classes—mostly the axis labels—we decided to keep the argument names in synch with those of the base classes.)

The main driver for making this breaking change was that one the main use cases for this class is to generate time-dependent spectra, and the basic signature that we picked in this case, as specified in the documentation, is spec(E, t)—where E (the energy) is the random variable and t (the time) the auxiliary variable. This makes sense, as in most cases the spectra are time-independent, and keeping the energy in the first position allows to give the time a default value, simplifying all the signatures in this fairly common use case. This also makes the signature of the pdf(rv, aux) class method more reasonable.

After the fix to the creation of the underlying meshgrid implemented in https://bitbucket.org/ixpesw/ixpeobssim/commits/cdbd99a7545bb1b0fd09591b2c9b973fb56325d4 insisting on the signature as xUnivariateAuxGenerator(aux, rv, pdf…) would have forced us to swap the spectrum arguments spec(E, t) -> spec(t, E) all over the place (e.g., in srcmodel.spectrum).

(In brief: we have rotated the underlying spline representation by 90 degrees with respect to the historical implementation.)

rvs(aux)

Return random variates for a given array of values of the auxiliary variable.

Note that the sample size, here, is determined by the size of the aux argument.

slice(aux)

Return the one-dimensional pdf for a given value of the auxiliary variable (i.e., a horizontal slice of the underlying bivariate spline).

class ixpeobssim.core.rand.xUnivariateGenerator(rv, pdf, w=None, bbox=None, k=3, ext=0, xlabel='rv', ylabel='pdf')

Univariate random number generator based on a linear interpolated spline.

This class is adding very little on top of the xInterpolatedUnivariateSpline spline class: essentially we are building the ppf in the constructor and we are adding the pdf() method (simply an alias to the underlying __call__() operator, as well as an rvs() method to generate actual random numbers).

Parameters:

• rv ((N,) array_like) – Array of points sampling the values of the random variable.

• pdf ((N,) array_like) – pdf values at the array rv.

• w ((N,) array_like, optional) – Weights for spline fitting.

• bbox ((2,) array_like, optional) – Bounduary of the approximation interval. .

• k (int) – The order of the spline (<=5, default 3, i.e., a cubic spline).

• ext (int or string, optional) – Controls the extrapolation mode for elements not in the interval defined by the knot sequence.

• xlabel (str, optional) – A text label of the random variable.

• ylabel (str, optional) – A text label for the pdf.

rvs(size=1)

Return random variates of arbitrary size.

Parameters:

size (int) – The number of random numbers to be generated.

rvs_bounded(size=1, rvmin=None, rvmax=None)

Return random variates of arbitrary size between the specified bounds.

This is an alternative to rvs(), where the uniformly-distributed random array passed to the ppf is spanning the appropriate subset of the [0, 1] interval.

Although the default behavior of this method is identical to the standard rvs(), we prefer to have a separate function to avoid unnecessary complexity in what is by far the most common use case.

class ixpeobssim.core.rand.xUnivariateGeneratorLinear(rv, pdf, ext=0, xlabel='rv', ylabel='pdf')

Subclass of xUnivariateGenerator restricted to a linear spline.

Spline utility module, building on top of the scipy.interpolate modules.

ixpeobssim.core.spline.interpolate(xa, ya, xb, yb, x)

Simple two-point linear interpolation/extrapolation.

This is a convenience function that is pretty much only used in the optimize_grid_linear() below.

ixpeobssim.core.spline.optimize_grid_linear(x, y, tolerance=0.0001)

Optimize a pair of (x, y) arrays for the corresponding spline definition.

This loops over the input arrays and removes unnecessary data points to minimize the length of the arrays necessary to the spline definition.

Parameters:

• x (array) – The input x-array.

• y (array) – The input y-array.

• tolerance (float) – The maximum relative difference between the generic yi value and the estrapolation of the two previous optimized data points for the point i to be removed.

class ixpeobssim.core.spline.xInterpolatedBivariateSpline(x, y, z, bbox=None, kx=3, ky=3, xlabel=None, ylabel=None, zlabel=None)

Bivariate interpolated spline on a rectangular grid.

This is somewhat similar in spirit to the corresponding univariate base class, except that the additional functionalities are, for the moment, limited to book-keeping and plotting facilities.

One handy facility that this class provides is that of allowing both an array like of shape (x.size, y.size) and a function as the z argument. This makes it possible to construct slides by passing either (x, y, z), as the signature of the underlying scipy’s interpolator, and (x, y, f(x, y)), which is easier and more intuitive in most cases.

Parameters:

• x (array_like) – Input x values (assumed to be sorted).

• y (array_like) – Input y values (assumed to be sorted).

• z (array_like of shape (x.size, y.size) or callable) – Input z values, with shape (x.size,y.size).

• bbox (array_like, optional) – The boundary of the approximation domain.

• kx (int) – The spline order on the two axes.

• ky (int) – The spline order on the two axes.

• xlabel (str, optional) – A text label for the quantity on the x-axis.

• ylabel (str, optional) – A text label for the quantity on the y-axis.

• zlabel (str, optional) – A text label for the quantity on the z-axis.

build_horizontal_ppf()

Create the horizontal percent point function (or inverse of cdf).

Warning

This really, really need to be fixed. Instead of grabbing a vertical slice at xmean, we should pass an argument to the function so that the subclasses can implement whatever is right for them.

hslice(y, k=3)

Return an horizontal slice at a given y of the bivariate spline.

Parameters:

• y (float) – The y value at which the horizontal slice should be calculated.

• k (int, optional) – The degree of the resulting spline.

plot(num_contours=75, logz=False, **kwargs)

Plot the spline.

plot_color(num_contours=75, logz=False, **kwargs)

Color plot

plot_contours(num_contours=10, colors='black', fontsize='small', logz=False, cfmt='%1.3f')

Contour plot.

scale(scale_factor, zlabel=None)

Scale the spline z values.

static transposed_meshgrid(x, y)

Return a suitable meshgrid to evaluate a function in the spline reference system.

This is needed because apparently meshgrid is using the matlab convention for the indices, which is inconsistent with most of the other interfaces in scipy. See: https://github.com/scipy/scipy/issues/3164

vslice(x, k=3)

Return a vertical slice at a given x of the bivariate spline.

Parameters:

• x (float) – The x value at which the vertical slice should be calculated.

• k (int, optional) – The degree of the resulting spline.

xmax()

Return the maximum of the underlying x-array.

xmin()

Return the minimum of the underlying x-array.

ymax()

Return the maximum of the underlying y-array.

ymin()

Return the minimum of the underlying y-array.

class ixpeobssim.core.spline.xInterpolatedBivariateSplineLinear(x, y, z, xlabel=None, ylabel=None, zlabel=None)

Bivariate linear interpolated spline on a rectangular grid.

class ixpeobssim.core.spline.xInterpolatedPiecewiseUnivariateSpline(x, y, breakpoints, w=None, bbox=None, k=3, ext=0, xlabel=None, ylabel=None)

Pieacewise interpolated spline.

In addition to all the regular xInterpolatedUnivariateSpline parameters, this takes an additional breakpoints argument that allows to split the spline creation and evaluation into independent pieces. This is achieved by stories a list of independent splines under the hood.

class ixpeobssim.core.spline.xInterpolatedUnivariateLogSpline(x, y, w=None, bbox=None, k=3, ext=0, xlabel=None, ylabel=None)

Poor man’s attempt at a spline in logarithmic space.

The class inherits from both xUnivariateSplineBase and scipy’s InterpolatedUnivariateSpline, and the basic strategy, here, is twofold: we initialize xUnivariateSplineBase from the physical x and y values, but we construct an actual interpolator in logarithmic space under the hood, by passing log10(x) and log10(y) to the InterpolatedUnivariateSpline constructor. We then overload the __call__() method by tranforming the argument into its logarithm, performing the interpolation in logarithmic space, and raising the result to the power 10.

The class comes with some pretty obvious shortcomings, among which the difficulty of implementing sensible derivatives and integrals. As far as the latter is concerned, we proceed by brute force and create a second spline, this time in linear space, with a relatively large number of points to guarantee the accuracy of the integral. Admittedly, this is not very elegant nor very roboust. For the derivative issue, we simply do not provide the nu argument in the overloaded __call__ method. For all the other scipy’s spline method that is not implemented in this context, we just overload it and raise a NotImplementedError exception.

At some point we even considered removing this class, but the truth is it is handy when one needs a quick way to extrapolate in log-log space. So the class is remaning, with all the caveates above.

antiderivative(n=1)

Overloaded method.

derivative(n=1)

Overloaded method.

derivatives(x)

Overloaded method.

integral(a, b)

Overloaded integral method.

The integral spline is calculated and cached the first time this method is called.

roots()

Overloaded method.

class ixpeobssim.core.spline.xInterpolatedUnivariateLogSplineLinear(x, y, w=None, bbox=None, ext=0, xlabel=None, ylabel=None)

Subclass of xInterpolatedUnivariateLogSpline with k=1.

antiderivative(n=1)

Overloaded method.

derivative(n=1)

Overloaded method.

derivatives(x)

Overloaded method.

roots()

Overloaded method.

class ixpeobssim.core.spline.xInterpolatedUnivariateSpline(x, y, w=None, bbox=None, k=3, ext=0, xlabel=None, ylabel=None)

Interpolated spline (i.e., passing through all the input points).

class ixpeobssim.core.spline.xInterpolatedUnivariateSplineLinear(x, y, w=None, bbox=None, ext=0, xlabel=None, ylabel=None)

Interpolate linear (k=1) spline.

class ixpeobssim.core.spline.xStepFunction(x, y, xlabel=None, ylabel=None)

Small convenience class describing a step function.

The basic rule, here, is that the value of the function between x[i] and x[i + 1] is exactly y[i].

Parameters:

• x (array_like) – The array of x values.

• y (array_like) – The array of y values, with shape (N, ), assuming that (N + 1, ) is the x shape.

• xlabel (str) – The text label for the x axis.

• ylabel (str) – The text label for the y axis.

plot(annotate=False, fmt='%.2f')

Plot the step function.

See the interesting discussion of a related problem at https://stackoverflow.com/questions/5347065 (For the reference, this is Will’s solution.)

class ixpeobssim.core.spline.xUnivariateSpline(x, y, w=None, bbox=None, k=3, s=None, ext=0, xlabel=None, ylabel=None, log=False)

Wrapper around scipy’s UnivariateSplineClass.

The basic idea is to keep track of the original arrays passed to the interpolator and to support arithmetic operations and plotting. We also allow the user to supply optional arguments to control the ranges and specify labels (e.g., names and units) for the quantities involved. We essentially maintain the original signature of the function, adding two optional arguments (xlabel and ylabel) at the end.

Note

Note that the interface to the base class has changed from numpy 0.14. An ext argument can be passed to the constructor starting with scipy 0.15 to control the extrapolation behavior and a check_finite argument is available in 0.16 to avoid nans in the input data. We currently do not expose either one and rely on the default parameters for backward compatibility.

Parameters:

• x ((N,) array_like) – Input x values for the spline.

• y ((N,) array_like) – Input y values for the spline.

• w ((N,) array_like, optional) – Weights for spline fitting.

• bbox ((2,) array_like, optional) – Bounduary of the approximation interval. .

• k (int) – The order of the spline (<=5, default 3, i.e., a cubic spline).

• s (float or None) – The spline smoothing factor (o means no smoothing).

• ext (int or string, optional) – Controls the extrapolation mode for elements not in the interval defined by the knot sequence.

• xlabel (str, optional) – A text label for the quantity on the x-axis.

• ylabel (str, optional) – A text label for the quantity on the y-axis.

• log (bool, default False) – Flag to initialize the underline scipy spline in logarithmic space.

build_cdf(ylabel='cdf', spline_class=None)

Create the cumulative distribution function.

Parameters:

• spline_class (class) – The specific spline class (e.g., xInterpolatedUnivariateSpline) used to build the cdf

• ylabel (str) – The (optional) label for the y axis

Note

At this point the dafault behavior is to use a linear spline for the cdf by default. This is mainly for historical reasons, but maybe a better option would be to use a cubic spline, or to delegate this to call self.__class__?

Warning

Re-enable the check on y > 0!

build_ppf(xlabel='q', spline_class=None)

Create the percent point function (or inverse of the cdf).

See the docstring for the build_cdf() method for more details.

inverse(xmin=None, xmax=None, **kwargs)

Calculate the inverse of the spline.

Note that the xmin and xmax arguments allow to invert the spline in a subset of its values.

Parameters:

• xmin (float) – The minimum value for the x-axis of the inverse spline (or y-axis of the original one).

• xmax (float) – The maximum value for the x-axis of the inverse spline (or y-axis of the original one).

norm()

Return the integral over the entire spline domain.

plot(num_points=1000, overlay=False, logx=False, logy=False, scale=1.0, offset=0.0, grids=False, **kwargs)

Plot the spline.

Parameters:

• num_points (int, optional) – The number of sampling points to be used to draw the spline.

• overlay (bool, optional) – If True, the original arrays passed to the spline are overlaid.

• logx (bool, optional) – If True, the spline is sampled and plotted with the log scale on the x axis.

• logy (bool, optional) – If True, the spline is plotted with the log scale on the y axis.

• scale (float, optional) – Optional scale factor for plotting (useful for overlaying splines with different ranges on the y axis).

• offset (float, optional) – Optional offset for plotting (useful for overlaying splines with different ranges on the y axis).

• **kwargs (dict) – Keyword arguments passed directly to the matplotlib.plot() method.

scale(scale_factor, **kwargs)

Return a different spline whose y values are scaled by a given factor wrt the original one.

xmax()

Return the maximum of the underlying x-array.

Note this works because the input array is sorted.

xmin()

Return the minimum of the underlying x-array.

Note this works because the input array is sorted.

class ixpeobssim.core.stokes.xDataStokesParameters

Small utility class to deal with the Stokes parameters in a data analysis context.

Warning

This is work in progress, and the class methods are known to fail by ZeroDivisionError in some circumstances.

static polarization_angle(Q, U, dQ, dU)

Return the polarization angle and proparate the uncertainties.

static polarization_degree(I, Q, U, dI, dQ, dU)

Return the polarization degree and proparate the uncertainties.

class ixpeobssim.core.stokes.xModelStokesParameters

Small utility class to deal with the Stokes parameters in a source model context.

Basically we provide conversion functions from Stokes parameters to polarizarion degree and angle and vice-versa.

Note that all the algebra, here, is coded in terms of the reduced Stokes parameters q = Q / I and u = U / I. The reason is twofold:

• when we simulate a model we typically decouple the definition of the spectrum from that of the polarization pattern;

• when converting polarization degree and angle to Stokes parameters we can only calculate, by definition, q and u—not Q and U.

For completeness: be aware that all the angles are measured in radians, and if you want to operate with degrees it is the user’s responsibility to do the conversion outside this class.

static normalize(QU, I)

Calcualate the Stokes normalized Q parameter.

static pdpa_to_xy(pol_deg, pol_ang, degrees=False)

Convert polarization degree and angle into the x and y components of the polarization vector.

This is assuming that the position angle is measured starting from the celestial North, see https://bitbucket.org/ixpesw/ixpeobssim/issues/597 for more discussion about this.

Warning

This nd the following function are encapsulating our convention for measuring position angles, and should be probably better suited in a different module?

static polarization_angle(q, u)

Convert q and u to the corresponding polarization angle (in radians).

static polarization_degree(q, u)

Convert q and u to the corresponding polarization degree.

static q(polarization_degree, polarization_angle)

Convert a polarization degree and angle (in radians) into the corresponding q reduced Stokes parameter.

static qu_to_xy(q, u)

Conver the Stokes parameters into the x and y components of the polarization vector.

static u(polarization_degree, polarization_angle)

Convert a polarization degree and angle (in radians) into the corresponding u reduced Stokes parameter.

Binned data products

Base classes for binned data products.

ixpeobssim.binning.base.peek_binning_algorithm(file_path)

Open a binned file and peek at the underlying algorithm

class ixpeobssim.binning.base.xBinnedFileBase(file_path)

Base class for binned files.

binning_algorithm()

Return the binning algorithm used to create the file.

du_id()

Return the DU ID for the binned file.

Added in version 12.0.0.

Warning

This was added in version 12 of ixpeobssim after https://bitbucket.org/ixpesw/ixpeobssim/issues/327

classmethod from_file_list(file_list)

Method to instance the binned class from a list of files.

ontime(scale=0.001)

Return the nominal ontime for the binned file.

Parameters:

scale (float) – A scale factor, by default 1.e-3 (i.e., ontime in ks).

plot(*args, **kwargs)

Do nothing method.

set_data(name, value)

Add a (key, value) pair to the underlying __data_dict class member.

Warning

This is yet another side effect of the perverse __getattr__ / __setattr__ mechanism we have put in place for this class—if we calculate a quantity dinamically for a binned object and we still want to be able to access it with the same rules, we have to manually add it to the dictionary.

write(file_path, overwrite=True)

Write the binned data to the output file.

class ixpeobssim.binning.base.xEventBinningBase(file_path, **kwargs)

Base class for the event binning.

This is essentially opening an event file and keeping track of the xEventFile object, along all the keyword arguments passed to the constructor.

bin_()

Do-nothing method to be reimplemented in the derived classes.

static bin_centers(bin_edges)

Return an array of bin centers given an array of bin edges.

Parameters:

bin_edges (1-d array of length (n + 1).) – The array with the bin edges.

Returns:

The array with the values of the bin centers.

Return type:

1-d array of length n.

static bin_widths(bin_edges)

Return an array of bin widths given an array of bin edges.

Parameters:

bin_edges (1-d array of length (n + 1).) – The array with the bin edges.

Returns:

The array with the values of the bin widths.

Return type:

1-d array of length n.

build_primary_hdu(data=None, header=None)

Build the primary HDU for the output file.

check_pcube_weighting_scheme(aeff)

Simple check on the weighting scheme being used.

For more background information, see https://bitbucket.org/ixpesw/ixpeobssim/issues/573 and https://bitbucket.org/ixpesw/ixpeobssim/issues/613

static equipopulated_binning(num_bins, data, min_val=None, max_val=None)

Create an equipopulated binning based on the values of a generic data column.

Parameters:

• num_bins (int) – The number of bins

• data (array) – The underlying data to be used for the binning

• min_val (float (optional)) – The minimum data value to be used for the binning

• max_val (float (optional)) – The maximum data value to be used for the binning

get(key, default=None)

Convenience method to address the keyword aguments.

load_aeff_for_polarization_analysis()

Load the proper arf file for a model-independent polarization analysis (i.e., to be used for polarization cubes, maps and map cubes).

This is loading the proper arf file making sure that, when weights are used, the SIMPLE weighting prescription is picked.

static make_binning(bin_alg, min_val=None, max_val=None, num_bins=None, bin_list=None, bin_data=None, bin_file=None)

Generic function to define a binning with several possible different algorithms.

static make_energy_binning(bin_data=None, **kwargs)

Small convenience function for energy binning.

process_image_ref_kwargs()

Set the xref and yref values for the output image (in sky) coordinates.

If either xref or yref are not specifified, the center of the ROI specified in the proper extension is assumed.

process_kwargs()

Check the keyword arguments and if the output file is not set, create a path on the fly, based on the event file and the binning algorithm.

static read_binning_from_file(file_path)

Read a custom binning from file and return a numpy array.

set(key, value)

Convenience method to set keyword arguments.

weight_data()

Retrieve the weights from the underlying event file.

This encapsulates the simple logic used downstream by the binning algorithms supporting weights, i.e., if the command-line argument weights is False the weights are all one, while if it is True returns they are the content of the column indicated by the weightcol command-line argument.

write_output_file(hdu_list, overwrite=True)

Write the binned data to the output file.

Binning data products in detector coordinates.

class ixpeobssim.binning.detector.xBinnedAreaEnergyFluxMap(file_path)

Display interface to binned EFLUX files.

class ixpeobssim.binning.detector.xBinnedAreaRateMap(file_path)

Display interface to binned ARMAP files.

plot()

Plot the data.

class ixpeobssim.binning.detector.xEventBinningARMAP(file_path, **kwargs)

Class for ARMAP binning.

bin_()

Overloaded method.

process_data()

Convenience function factoring out the code in common with the corresponding EFLUX class—see the overloaded method in there.

Here we are binning the event position in detector coordinates and dividing by the livetime and the bin area. The function returns a n x n array of area rate values to be written in the output file.

class ixpeobssim.binning.detector.xEventBinningEFLUX(file_path, **kwargs)

Class for EFLUX binning.

process_data()

Overloaded method.

Here we are just multiplying by the average measured event energy.

facilities for exposure calculation.

ixpeobssim.binning.exposure.create_psf_kernel(psf, pixsize, npix=21)

Generate a binned kernel for the PSF convolution.

Warning

There is potential overlap with the calculate_psf_kernel() function in the evt.deconvolution module, but the latter is doing something slightly different—a Monte Carlo integration over the pixels, rather than a direct evaluation at the pixel center. We should probably think about whether which one we want.

class ixpeobssim.binning.exposure.xBinnedLivetimeCube(file_path)

Read-mode interface to a LTCUBE FITS file.

TODO: Temporarily most of this class is just a copy of xBinnedMDPMapCube.

We should refactor the common methods in a separate class.

ds9_region_mask(region_list, region_slice=None)

Return the (spatial) array map corresponding to a given ds9 region list.

If there is more than one region in the region list, by default the output mask corresponds to the logical or of all the regions. The region_slice optional argument allows to restrict the mask to a subset of the regions.

Parameters:

• region_list – The region list defining the mask

• region_slice (int or slice (optional)) – An optional slice designator to select a subset of the region list.

map_shape()

Return the shape of the underlying sky-maps.

Mind the underlying arrays are all 3-dimensional cubes with the theta binning as the first dimension—so it’s the last two that we care about. (And, for completeness: since the arrays are all the same, we use I for convenience.)

num_theta_layers()

Return the number of theta layers in the cube.

pixel_size()

Return the pixel size of the underlying spatial map.

plot(prefix=None)

Plot the data.

plot_livetime_map(theta_layers=None, prefix=None, **kwargs)

Plot the MDP map.

sum_pixels(spatial_mask, theta_layer=0)

Sum the relevant quantities over a given spatial mask and theta slice.

theta_binning()

Return the underlying theta binning in the form of a numpy array.

Note the array has one more element than the underlying THETA_LO and THETA_HI arrays, i.e., if the cube is binned in two theta layers, say 0–2 arcmin and 2–4 arcmin, the theta binning is [0. 2. 4.].

theta_centers()

Return the centers of the binning over theta.

theta_label(theta_layer)

Return a text label corresponding to the theta layer (e.g, to identify the theta range on a plot).

theta_range(theta_layer)

Return the theta value for a given theta layer.

class ixpeobssim.binning.exposure.xEventBinningLTCUBE(file_path, **kwargs)

Class for LTCUBE binning.

bin_()

Overloaded method.

process_kwargs()

Overloaded method.

class ixpeobssim.binning.exposure.xExposureCube(exposure, ebounds=None, units=None)

Structure for writing/reading exposure cubes.

classmethod empty(shape)

Create an empty exposure cube.

classmethod from_file(file_path)

Read an exposure cube from file.

write(file_path, header, overwrite)

Write the exposure cube to file.

Warning

Need to add all the keywords from the LTCUBE file

Format definitions for binned data products.

class ixpeobssim.binning.fmt.xBinTableHDUEBOUNDS(data=None, keywords=None, comments=None)

Binary table for storing energy bounds.

Warning

Can we just reuse the same extension from the response matrix?

class ixpeobssim.binning.fmt.xBinTableHDULC(data=None, keywords=None, comments=None)

Binary table for binned LC data.

class ixpeobssim.binning.fmt.xBinTableHDUPCUBE(data=None, keywords=None, comments=None)

Binary table for binned PCUBE data.

class ixpeobssim.binning.fmt.xBinTableHDUPHA1(data=None, keywords=None, comments=None)

Binary table for binned PHA1 data.

class ixpeobssim.binning.fmt.xBinTableHDUPP(data=None, keywords=None, comments=None)

Binary table for binned PP data.

class ixpeobssim.binning.fmt.xBinTableHDUTHETABOUNDS(data=None, keywords=None, comments=None)

Binary table for storing off-axis angle bounds.

Miscellanea binned products.

class ixpeobssim.binning.misc.xBinnedLightCurve(file_path)

Binned light curve.

plot(mjd=False, **plot_opts)

Overloaded plot method.

rate()

rate_error()

class ixpeobssim.binning.misc.xBinnedMap(file_path)

Display interface to binned CMAP files.

While this is essentially an xFITSImage, we need to inherit from xBinnedFileBase, as the latter provides the implementation of the from_file_list() slot. Instead of inheriting from xFITSImage, too, we preferred composition, instead.

Warning

I am more and more convinced that this is poor design, and maybe the pylint super-init-not-called is there for a reason. We should probably refactor the class and make the implementation cleaner.

plot(**kwargs)

Plot the data. The kwargs passed to plt.imshow().

class ixpeobssim.binning.misc.xBinnedPulseProfile(file_path)

Binned pulse profile.

plot(**kwargs)

Overloaded plot method.

class ixpeobssim.binning.misc.xEventBinningCMAP(file_path, **kwargs)

Class for CMAP binning.

bin_()

Overloaded method.

process_kwargs()

Overloaded method.

class ixpeobssim.binning.misc.xEventBinningLC(file_path, **kwargs)

Class for LC binning.

bin_()

Overloaded method.

make_binning()

Build the light-curve binning.

process_kwargs()

Overloaded method.

class ixpeobssim.binning.misc.xEventBinningPP(file_path, **kwargs)

Class for pulse-profile binning.

bin_()

Overloaded method.

make_binning()

Build the light-curve binning.

Binned data structures pertaining to the (spectro-)polarimetric analysis.

class ixpeobssim.binning.polarization.xBinnedCountSpectrum(file_path)

Binned count spectrum.

plot(**kwargs)

Overloaded plot method.

respfile()

Return the value of the RESPFILE keyword in the SPECTRUM extension header, stripped from the first part of the pathe (i.e., the file name only).

class ixpeobssim.binning.polarization.xBinnedMDPMapCube(file_path)

Read-mode interface to a MDPMAPCUBE FITS file.

ds9_region_mask(region_list, region_slice=None)

Return the (spatial) array map corresponding to a given ds9 region list.

If there is more than one region in the region list, by default the output mask corresponds to the logical or of all the regions. The region_slice optional argument allows to restrict the mask to a subset of the regions.

Parameters:

• region_list – The region list defining the mask

• region_slice (int or slice (optional)) – An optional slice designator to select a subset of the region list.

energy_binning()

Return the underlying energy binning in the form of a numpy array.

Note the array has one more element than the underlying ENERG_LO and ENERG_HI arrays, i.e., if the cube is binned in two energy layers, say 2–4 keV and 4–8 keV, the energy binning is [2. 4. 8.].

energy_label(energy_layer)

Return a text label corresponding to the energy layer (e.g, to identify the energy range on a plot).

energy_range(energy_layer)

Return the energy range, i.e., (emin, emax) in keV, for a given energy layer.

map_shape()

Return the shape of the underlying sky-maps.

Mind the underlying arrays are all 3-dimensional cubes with the energy binning as the first dimension—so it’s the last two that we care about. (And, for completeness: since the arrays are all the same, we use I for convenience.)

num_energy_layers()

Return the number of energy layers in the cube.

pixel_size()

Return the pixel size of the underlying spatial map.

plot(prefix=None)

Plot the data.

plot_mdp_map(energy_layers=None, prefix=None, **kwargs)

Plot the MDP map.

sum_pixels(spatial_mask, energy_layer=0)

Sum the relevant quantities over a given spatial mask and energy slice.

class ixpeobssim.binning.polarization.xBinnedPolarizationCube(file_path)

Read-mode interface to a PCUBE FITS file.

Added in version 12.0.0.

as_table()

Return the polarization cube as an astropy table.

backscal()

Return the value of the BACKSCAL header keyword, if present.

plot(**kwargs)

Default plotting for the polarization cubes.

Warning

This should be considered in evolution, and the API might change.

This is plotting the polarization cube in normalized Stokes parameters space, with custom grids to help reading out the polarization degree and angle.

A brief explanation of the supporte keywords arduments.

Parameters:

• side (float, optional) – The absolute value of the maximum Q/I and U/I to be displayed (note the aspect ratio of the plot is set to ‘equal’, and the plot is assumed to be squared). By default this is driven by the maximum value of the polarization degree over the polarization cube.

• pd_grid (array_like) – The polarization degree radial grid in correspondence of which the gridding circumpherences are plotted to guide the eye.

• pd_grid_label_angle (float) – The angle at which the text labels for the PD grids are plotted.

• pa_grid_step (float) – The step of the tangential grid in the polarization angle.

• sigma_levels (array_like) – The sigma levels at which the ellipses repesenting the normalized Stokes parameter contours are plotted.

• sigma_ls (tuple of the same size as sigma_levels) – The line styles corresponding to sigma levels—by default the innermost ellipse is solid and the others are dashed.

• colors (array_like or str, optional) – The colors for the elliptical contours (by deafult the proper color for the specific DU in the polarization cube is picked).

• label (str, optional) – Optional label to be displayed in the legend—this is only used if the marker boolean flag is True and only associated to the first layer in the cube.

• marker (bool) – If True, a marker is plotted at the center of the elliptical contours.

• marker_size (float) – The size for the optional markers.

• annotate_energies (bool) – If true, the energy layers are annotated with nice arrows and text labels (note the positioning is still fragile).

• setup_axes (bool) – Call the underlying setup_gca_stokes() hook at the end. (This flag is provided to allow disingaging multiple grid plotting in case one wants to overlay multiple polarization cubes.)

plot_polarization_angle(**kwargs)

Plot the polarization angle as a function of energy.

Warning

This is deprecated in favor of the standard plots of the Stokes parameters.

plot_polarization_degree(min_num_sigmas=2.0, **kwargs)

Plot the polarization degree as a function of energy.

Warning

This is deprecated in favor of the standard plots of the Stokes parameters.

polarization()

Return the polarization information in the form of a four element tuple (PD, PD_err, PA, PA_err) of arrays.

class ixpeobssim.binning.polarization.xBinnedPolarizationMapCube(file_path)

Read-mode interface to a PMAPCUBE FITS file.

align()

Align the polarization direction to the radial direction on a bin by bin basis.

Warning

This will need to be reviewed when we fix issue #597.

calculate_significance_mask(num_sigma=2.0, intensity_percentile=0.0)

Utitlity function to calculate a mask on the underlying pixel cube where the polarization degree and angle can be reliably plotted.

The calculation is based on two different metrics:

• the value of I in each given pixel, compared with a given percentile of the overall I distribution across the cube;

• the number of sigma that the measured polarization degree differs from either 0 or 1 (if the measurement is too close to 0 or 1 and misses the physical bounds by less than a few error bars, then it is effectively not a measurement).

Note the logic is positive—i.e., we’re passing the mask identifying the good pixels. This is typically negated downstream to set to zero all the other pixels.

Parameters:

• num_sigma (float) – The number of standard deviations representing the significance of the measurement of the polarization degree.

• intensity_percentile (float) – the threshold on the I value for the pixels, based on the percentile of the I distribution across all the pixels.

convolve(kernel)

Convolve the polarization cube with a generic binned kernel.

normalized_stokes_parameters()

Calculate the normalized Stokes parameters, i.e., Q/I and U/I.

plot(prefix=None)

Plot everything.

plot_normalized_stokes_parameters(energy_layers=None, num_sigma=1.0, intensity_percentile=0.05, prefix=None, **kwargs)

Plot the normalized Stokes parameters.

plot_polarization_angle(energy_layers=None, num_sigma=2.0, prefix=None, **kwargs)

Plot the polarization angle.

plot_polarization_degree(energy_layers=None, num_sigma=2.0, arrows=True, prefix=None, **kwargs)

Plot the polarization degree.

plot_significance(energy_layers=None, prefix=None, **kwargs)

Plot the significance.

plot_stokes_parameters(energy_layers=None, prefix=None, **kwargs)

Plot the Stokes parameters.

radial_profile(n=10, layer=0)

Create a radial polarization profile of the source.

Warning

This needs to be cleaned up and properly documented.

save_to_ds9(output_file_path, num_sigma=2.0, intensity_percentile=0.0, line_color='white')

Save the polarization arrows to a ds9 region.

sum_pixels(spatial_mask, energy_layer=0)

Sum the relevant quantities over a given spatial mask and energy slice.

class ixpeobssim.binning.polarization.xEventBinningMDPMAP(file_path, **kwargs)

Class for MDPMAP binning.

process_kwargs()

Overloaded method.

Note we have to call the method of the base class closest in the inheritance hierarchy, to make sure that all the ingredients for the WCS object generation are properly setup.

class ixpeobssim.binning.polarization.xEventBinningMDPMAPCUBE(file_path, **kwargs)

Class for MDPMAPCUBE binning.

bin_()

Overloaded method.

process_kwargs()

Overloaded method.

class ixpeobssim.binning.polarization.xEventBinningPCUBE(file_path, **kwargs)

Class for PCUBE binning.

Added in version 12.0.0.

bin_()

Overloaded method.

class ixpeobssim.binning.polarization.xEventBinningPHA1(file_path, **kwargs)

Original algorithm for PHA1 files.

binned_data()

Overloaded binning algorithm.

class ixpeobssim.binning.polarization.xEventBinningPHA1Base(file_path, **kwargs)

Base class for PHA1 binning.

bin_()

Overloaded method.

binned_data()

Method to be overloaded by derived classes.

class ixpeobssim.binning.polarization.xEventBinningPHA1Q(file_path, **kwargs)

Class for PHA1 binning.

binned_data()

Overloaded binning algorithm.

class ixpeobssim.binning.polarization.xEventBinningPHA1QN(file_path, **kwargs)

Subclass for creating normalized Stokes Q spectra.

binned_data()

Overloaded binning algorithm.

class ixpeobssim.binning.polarization.xEventBinningPHA1U(file_path, **kwargs)

Subclass for creating Stokes U spectra.

binned_data()

Overloaded binning algorithm.

class ixpeobssim.binning.polarization.xEventBinningPHA1UN(file_path, **kwargs)

Subclass for creating normalized Stokes U spectra.

binned_data()

Overloaded binning algorithm.

class ixpeobssim.binning.polarization.xEventBinningPMAP(file_path, **kwargs)

Class for PMAP binning.

process_kwargs()

Overloaded method.

Note we have to call the method of the base class closest in the inheritance hierarchy, to make sure that all the ingredients for the WCS object generation are properly setup.

class ixpeobssim.binning.polarization.xEventBinningPMAPCUBE(file_path, **kwargs)

Class for PMAPCUBE binning.

Event-level analysis

Facilities for the rotation of the polarization angle/Stokes parameters given an input model.

This module provides a series of functions facilitating the search for large-scale polarization signatures, such as radial/tangential polarization in extended sources such SNRs.

ixpeobssim.evt.align.align_phi(phi, phi0)

Rotate the photoelectron azimuthal angle and recalculate the event-by-event Stokes parameters.

Parameters:

• phi (array_like) – The azimuthal angle for the input event list

• phi0 (array_like) – The model polarization direction, calculated at the positions of the input events.

• warning:: (..) – As we moved away from the used of azimuthal angles in the analysis, this function was added to the module for backward compatibility and testing purposes, but should no be used. Use align_stokes_parameters() below, instead.

ixpeobssim.evt.align.align_stokes_parameters(q, u, q0, u0)

Align the Stokes parameters according to an input polarization model.

Parameters:

• q (array_like) – The event-by-event Q Stokes parameters.

• u (array_like) – The event-by-event U Stokes parameters.

• q0 (array_like) – The input model Q Stokes parameters, calculated at the positions of the input events.

• u0 (array_like) – The input model U Stokes parameters, calculated at the positions of the input events.

Animation utilities.

class ixpeobssim.evt.animate.xMovingCircle(points, radius=25.0, kxy=1, kr=1)

Small convenience class representing a time-dependent circular patch.

The moving patch is initialized given a set of (t, x, y) points, referred to the center of the target image. The interpolation is done via interpolating splines of order 1.

Parameters:

• points (iterable of three-elements (t, x, y) tuples) – The points defining the path of the circle as a function of time. (The time is in seconds and the xy position is expressed as the delta with respect to the center of the image in arcminutes.)

• radius (the circle radius in arcseconds)

event_mask(t, ra, dec, ra0, dec0)

Return a boolean mask for a series of events.

Parameters:

• t (array_like) – The array of event times.

• ra (array_like) – The array of the event R.A. positions.

• dec (array_like) – The array of the event Dec. positions.

• ra0 (float) – The R.A. position of the target image center.

• dec0 (float) – The Dec. position of the target image center.

frame_positions(interval=100)

Return a set of frame positions, referred to the center of the target image, for given animation interval.

The positions are in the form (dx, dy, radius), where dx and dy are the distances to the image center along the two orthogonal axes, in arcmin, and radius is the circle radius in arcseconds—this has to match the signature of the xSkyAnimation.show_roi() method, which is used as the updating slot of the animation.

See https://matplotlib.org/stable/api/animation_api.html for more information about the matplotlib animation APIs.

Parameters:

interval (the animation interval in ms.)

class ixpeobssim.evt.animate.xSkyAnimation(file_path)

Small convenience class to describe an animation in sky coordinates.

add_roi(roi, **kwargs)

Add the ROI.

plot(**kwargs)

Plot the underlying image.

Note that, in addition to dispatch the plot() call to the underlying xFITSImageBase object, this is overriding the return value, returning an AxesImage object, rather than the standard figure. This is needed to properly build the animation video.

run(roi, interval=100, figsize=(12.0, 12.0))

Run the actual animation.

See https://matplotlib.org/stable/api/animation_api.html for more information about the matplotlib animation APIs.

See also https://holypython.com/how-to-save-matplotlib-animations-the-ultimate-guide/

Parameters:

interval (the animation interval in ms.)

show_frame(roi=(0.0, 0.0, 25.0), **kwargs)

Draw a circular ROI with a given position and radius.

Parameters:

roi (three-element tuple) – The position and radius of the circular ROI in the form of a three-element tuple (dx, dy, radius).

Module encapsulating image deconvolution facilities.

ixpeobssim.evt.deconvolution.calculate_inverse_kernel(intensity, true_intensity, kernel)

Calculate the gargantuan (nx x ny x nk x nk) matrix that allows to deconvolve an (I, U or Q) image, given a proxy of the true intensity.

The basic idea, here, is that if you have a sensible proxy of the true I map (e.g., a high-resolution image from Chandra) you can use it to gauge the fraction of the measured intensity ending up in each true pixel, and you can actually use that to deconvole a map (with the same binning) of either Q or U. Phrased in a slighlty different way, while the PSF kernel for the direct true -> measured convolution is identical for all pixels, the kernel for the inverse measured -> true deconvolution is in general different for all the pixels and depends on the true intensity. This is why we need to precompute a store a specific kernel for each pixels. This can be in turn use to deconvole any measured quantity with the same spatial binning (Q, U, W2).

Warning

I am not positive this is right for the border pixels, where the kernel mask ends up out of the image—need to study that.

Parameters:

• intensity (2d array) – The measured intensity map.

• true_intensity (2d array) – The proxy for the true intensity map.

• kernel (2d array) – The PSF kernel for the convolution.

ixpeobssim.evt.deconvolution.calculate_psf_kernel(pixel_size, irf_name='ixpe:obssim20240101:v13', du_id=1, max_radius=0.03, sample_size=5000000)

Calculate the digital PSF kernel for image convolution.

This is essentially creating a square grid whose pixel size is set to the value of the corresponding function argument and whose half size is (loosely) determined by the max_radius arguments. A number of points is thrown in the grid according to the PSF shape and the corresponding probability for each pixel is estimated by numerical integration.

Parameters:

• pixel_size (float) – The pixels size of the image that the kernel is operating on (in decimal degrees).

• irf_name (str) – The IRF name.

• du_id (int) – The DU id.

• max_radius (float) – The maximum radius (in decimal degrees) for defining the size of the kernel grid. This should be large enough to fully contain the PSF.

• sample_size (int) – The size of the random sample for the Monte Carlo integration of the psf profile over the square grid.

ixpeobssim.evt.deconvolution.circular_kernel(nside, normalize=False)

Build a simple circular kernel of a given side.

Parameters:

• nside (int) – The side of the kernel array, assumed to be squared.

• normalize (bool) – If True, normalize to the kernel content (i.e., use for averaging rather than summing).

ixpeobssim.evt.deconvolution.convolve_image(image, kernel)

Convolve an image with the PSF kernel.

While I am sure there exist a more clever (and faster) way to do this with scipy, I retained the explicit Python loop, here, as I don’t think this is a bottleneck in any way, and having a full implementation is useful for debugging.

Parameters:

• image (2d array) – The input image to be convolved with the PSF.

• kernel (2d array) – The PSF kernel for the convolution.

ixpeobssim.evt.deconvolution.deconvolve_image(image, kernel, K)

Deconvolve a generic image (e.g., Q or U).

Parameters:

• image (2d array) – The input image to be deconvoled.

• kernel (2d array) – The PSF kernel.

• K (4d array) – The inverse pixel-by-pixel kernel for the deconvolution.

Event display facilities.

This module provides a simple, top-level interface to track images in IXPE level-1 files.

class ixpeobssim.evt.display.Recon(absorption_point: tuple[float, float], barycenter: tuple[float, float], track_direction: float, length: float, width: float)

Container class encapsulating the event reconstruction.

draw_absorption_point()

Draw the reconstructed absorption point and track direction.

draw_barycenter()

Draw the reconstructed absorption point and track direction.

draw_track_direction(line_width=1.0, length_ratio=0.5)

Draw the track direction.

ixpeobssim.evt.display.display_event(event, grid, threshold, dbscan, file_name=None, padding=False, **kwargs)

Single-stop event display.

ixpeobssim.evt.display.load_event_list(file_path, pivot_energy=8.0, interactive=False, **kwargs)

Load the event data from the Level-2 event list for the purpose of the event display—these include, in order: mission elapsed time, energy, sky position and Stokes parameters.

This function has a few other functionalities, other than just loadinf the relevant colums from the the Level-2 files, and particularly:

• resample the input events with a given power-law spectral function;

• trimming the resampled colums to a target number of events preserving the time ordering and covering evengly the entire time span.

Parameters:

• file_path (str) – The path to the input Level-2 file.

• pivot_energy (float) – The pivot energy for the resampling of the count spectrum.

• interactive (bool) – If True, show some debug plot with the output (resampled) spectrum.

• kwargs (dict) – The keyword arguments from the xDisplayArgumentParser.

class ixpeobssim.evt.display.xDisplayArgumentParser(description)

Specialized argument parser for the event display and related facilities.

This is placed here because if needs to be used by both the single-event display and the observation carousel.

class ixpeobssim.evt.display.xDisplayCard(target_name, header)

Specialize text card to display event information.

The basic idea, here, is that one initializes the card with the EVENTS header of a Level-2 file, and then updates the information on an event-by-event basis using the set_event_data() hook.

set_event_data(met, energy, ra, dec, q, u, compact=True)

Set the event data.

update_cumulative_statistics(num_events, emin, emax)

Set the card line with the basic cumulative statistics info.

class ixpeobssim.evt.display.xHexagonCollection(x, y, radius=0.05, orientation=0.0, **kwargs)

Collection of native hexagin patches.

Parameters:

• x (array_like) – The x coordinates of the hexagon centers.

• y (array_like) – The y coordinates of the hexagon centers.

• radius (float) – The hexagon apothem.

• orientation (float) – The hexagon orientation in radians—zero means pointy topped.

• kwargs – The keyword arguments to be passed to the PatchCollection constructor.

set(*, agg_filter=<UNSET>, alpha=<UNSET>, animated=<UNSET>, antialiased=<UNSET>, array=<UNSET>, capstyle=<UNSET>, clim=<UNSET>, clip_box=<UNSET>, clip_on=<UNSET>, clip_path=<UNSET>, cmap=<UNSET>, color=<UNSET>, edgecolor=<UNSET>, facecolor=<UNSET>, gid=<UNSET>, hatch=<UNSET>, in_layout=<UNSET>, joinstyle=<UNSET>, label=<UNSET>, linestyle=<UNSET>, linewidth=<UNSET>, mouseover=<UNSET>, norm=<UNSET>, offset_transform=<UNSET>, offsets=<UNSET>, path_effects=<UNSET>, paths=<UNSET>, picker=<UNSET>, pickradius=<UNSET>, rasterized=<UNSET>, sketch_params=<UNSET>, snap=<UNSET>, transform=<UNSET>, url=<UNSET>, urls=<UNSET>, visible=<UNSET>, zorder=<UNSET>)

Set multiple properties at once.

Supported properties are

Properties:

agg_filter: a filter function, which takes a (m, n, 3) float array and a dpi value, and returns a (m, n, 3) array and two offsets from the bottom left corner of the image alpha: array-like or scalar or None animated: bool antialiased or aa or antialiaseds: bool or list of bools array: array-like or None capstyle: .CapStyle or {‘butt’, ‘projecting’, ‘round’} clim: (vmin: float, vmax: float) clip_box: ~matplotlib.transforms.BboxBase or None clip_on: bool clip_path: Patch or (Path, Transform) or None cmap: .Colormap or str or None color: :mpltype:`color` or list of RGBA tuples edgecolor or ec or edgecolors: :mpltype:`color` or list of :mpltype:`color` or ‘face’ facecolor or facecolors or fc: :mpltype:`color` or list of :mpltype:`color` figure: ~matplotlib.figure.Figure gid: str hatch: {‘/’, ‘\’, ‘|’, ‘-’, ‘+’, ‘x’, ‘o’, ‘O’, ‘.’, ‘*’} in_layout: bool joinstyle: .JoinStyle or {‘miter’, ‘round’, ‘bevel’} label: object linestyle or dashes or linestyles or ls: str or tuple or list thereof linewidth or linewidths or lw: float or list of floats mouseover: bool norm: .Normalize or str or None offset_transform or transOffset: .Transform offsets: (N, 2) or (2,) array-like path_effects: list of .AbstractPathEffect paths: unknown picker: None or bool or float or callable pickradius: float rasterized: bool sketch_params: (scale: float, length: float, randomness: float) snap: bool or None transform: ~matplotlib.transforms.Transform url: str urls: list of str or None visible: bool zorder: float

class ixpeobssim.evt.display.xHexagonalGrid(num_cols, num_rows, pitch=0.05, **kwargs)

Generic hexagonal grid.

Parameters:

• num_cols (int) – The number of columns in the grid

• num_rows (int) – The number of rows in the grid

• pitch (float) – The grid pitch in mm.

static brightness(color)

Quick and dirty proxy for the brighness of a given array of colors.

See https://stackoverflow.com/questions/9733288 and also https://stackoverflow.com/questions/30820962 for how to split in columns the array of colors.

default_roi_side(roi, min_side, pad=0.1)

Return the default physical size of the canvas necessary to fully contain a given ROI.

draw_event(event, num_clusters=1, offset=(0.0, 0.0), min_canvas_side=2.0, indices=True, padding=True, zero_sup_threshold=None, values=False, **kwargs)

Draw an actual event int the parent hexagonal grid.

This is taking over where the draw_roi() hook left, and adding the event part.

draw_roi(roi, offset=(0.0, 0.0), indices=True, padding=True, **kwargs)

Draw a specific ROI of the parent grid.

pha_to_colors(pha, zero_sup_threshold=None)

Convert the pha values to colors for display purposes.

pixel_to_world(col, row)

Transform pixel to world coordinates.

Parameters:

• col (array_like) – The input column number(s).

• row (array_like) – The input row number(s).

roi_center(roi)

Return the world coordinates of the physical center of a given ROI.

Note the small offsets that we apply serves the purpose of taking into account the fact that pixels are staggered in one direction.

static show_display(file_path=None, dpi=100, batch=False)

Convenience function to setup the matplotlib canvas for an event display.

Parameters:

file_path (str) – Optional file path to save the image immediately before the plt.show() call.

class ixpeobssim.evt.display.xL1Event(min_col: int, max_col: int, min_row: int, max_row: int, pha: array, trigger_id: int = 0, seconds: int = 0, microseconds: int = 0, timestamp: float = 0.0, livetime: int = 0, error_summary: int = 0, du_status: int = 0, recon: Recon | None = None)

A fully fledged event.

This is building up on the core logic encapsulated in the xRegionOfInterest class, and is adding all the necessary event information on top of that, i.e., the pha values, and the time and error information.

Note the dataclass field are largely mapped over the columns of the corresponding EVENTS extension in the underlying FITS files.

highest_pixel(absolute=True)

Return the coordinates (col, row) of the highest pixel.

Parameters:

absolute (bool) – If true, the absolute coordinates (i.e., those referring to the readout chip) are returned; otherwise the coordinates are intended relative to the readout window (i.e., they can be used to index the pha array).

run_clustering(engine)

Run the clustering on the track image.

class ixpeobssim.evt.display.xL1EventFile(file_path)

Simple interface to a level-1 file in FITS format.

Parameters:

• file_path (str) – The path to the input event file.

• padding (Padding instance) – The ROI padding for the event file.

bisect_met(met)

Retrieve a specific event by its mission elapsed time.

Internally this is using a binary search on the time column, and in general it can be assumed that this O(log(N)) in complexity.

value(col_name)

Return the value of a given column for a given extension for the current event.

zero_sup_threshold()

Return the zero-suppression threshold, determined from the file header.

class ixpeobssim.evt.display.xRegionOfInterest(min_col: int, max_col: int, min_row: int, max_row: int)

Class describing a region of interest (ROI).

A region of interest is the datum of the logical coorinates of its two extreme corners, in the order (min_col, max_col, min_row, max_row).

at_border()

Return True if the ROI is on the border for a given chip_size.

column_indices()

Return an array with all the valid column indices.

coordinates_in_roi(col, row)

Return a boolean mask indicaing whether elements of the (col, row) arrays lie into the ROT area.

coordinates_in_rot(col, row)

Return a boolean mask indicaing whether elements of the (col, row) arrays lie into the ROT area.

row_indices()

Return an array with all the valid row indices.

serial_readout_coordinates()

Return two one-dimensional arrays containing the column and row indices, respectively, in order of serial readout of the ROI.

Example

>>> col, row = xRegionOfInterest(0, 4, 0, 3).serial_readout_coordinates()
>>> print(col)
>>> [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]
>>> print(row)
>>> [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]

serial_readout_indices()

Return a two-dimensional array containing the readout index for each pixel of the ROI.

The ASIC serial readout starts from the top-left corner and proceeds one row at a time, i.e., for a toy 5 x 5 window starting at <0, 0> the readout indices look like +——————– | 0 1 2 3 4 | | 5 6 7 8 9 | | 10 11 12 13 14 | | 15 16 17 18 19 | | 20 21 22 23 24

It is worth noting that, provided that one loops over the row indices first and column indices last, the readout index can be determined by just accumulating a counter. This function is useful as it provides the right answer for any position in the event window, no matter what the loop order is.

class ixpeobssim.evt.display.xXpolGrid(**kwargs)

XPOL grid.

Module encapsulating the event structure and related facilities.

class ixpeobssim.evt.event.xBaseEventList(time_=None, source_id=None)

Base class for event lists.

This is the base class for objects that is created and filled at simulation time and then written to file in FITS format. It is essentially a dictionary containing all the column data for the output FITS file, providing additional facilities for addding, sorting in time and trimming event lists.

In order to be robust against name changes, subclasses should provide all the necessary interfaces to retrieve and set the column values (i.e., if the name of any of the columns changes, in principle it should not be necessary to make modifications beyond this file). Note that the colums related to times and source ID are encapsulated into this base class.

If the constructor is called with no argument an empty event list is created. Alternatively the user can initilize the event list with an array of event times (we did pick the time because typically this is the first dynamical variable being calculated). In that case the TRG_ID column is also filled right away and, if the optional source_id argument is not None, the SRC_ID as well. For completeness, this mechanism was added because we noticed that there were many calls to xEventList.fill_trigger_id() and xEventList.set_source_id() floating around in different places. (Note that, altough we need to overwrite the trigger identifier values at the end of the loop over the model components, after we have applied the deadtime, we do need to fill the corresponding column at the level of the single component because all the code downstream, e.g., the xEventList.trim() method, need all the columns to be present, and have the same length, to operate.)

apply_fiducial_area_base(detx, dety)

Trim out the events falling outside the detector active area.

Note that sub-classes are responsible for providing their own definition of the detector coordinates to apply the function.

apply_vignetting_base(ra, dec, energy, vign, ra_pnt, dec_pnt)

Base function to apply the effective area vignetting to the list.

Note that sub-classes are responsible for providing their own definition of the sky coordinates and energy to apply the function.

Note

Note that angular_separation returns values in degrees, while the vignetting function is expressed in terms of the off-axis angle expressed in arcminutes. See also https://bitbucket.org/ixpesw/ixpeobssim/issues/423

static array_is_sorted(array)

Return whether a given array is sorted.

This is an interesting topic per se, see https://stackoverflow.com/questions/47004506 I dropped this here assuming that we’re not concerned by adding a step which is O(n) from an algorithmic standpoint.

empty()

Return whether the event list is empty.

num_events()

Return the number of events in the event list.

set_source_id(src_id)

Set the source ID.

set_time_columns(time_, check_sorted=True)

Populate the time-related columns.

sort()

Sort the event list based on the event time.

static split_event_time(time_)

Split the event time into its interger and fractional part.

While in hardware we will be measuring the seconds and microseconds separately, and then build the floating-point representation of the event time by summing the two, in the simulation realm we do the opposite, i.e., we extract the event times as floating-point numbers and then reverse-engineer the integral and fractional parts.

Warning

This is a crude implementation of the concept, where we essentially take the floor of the integral and fractional parts of the event time (modulo the conversion from seconds to microseconds for the second). In the future we might make this more complicated if we need to simulate the instrument timing in a more detailed fashion.

src_id()

Return the source ID.

time()

Return the event times.

trim(mask)

Trim all the data columns according to a generic input mask.

class ixpeobssim.evt.event.xEventFile(file_path)

Read/write interface to level-2 event files.

This is the main class providing access to the information into a FITS file containing an event list.

It supports minimal output facilities, such as adding columns and writing to file.

add_column(ext_name, col_name, col_format, col_data)

Add a column to the specified extension of the event file.

Parameters:

• ext_name (str) – The name of the extension that the column should be added to.

• col_name (str) – The name of the column to be added.

• col_format (str) – The format for the new column (e.g., ‘E’).

• col_data (array_like) – The data for the new column.

add_columns(ext_name, *columns)

Add multiple columns to the specified extension of the event file.

Compared to multiple calls to add_column(), and at the cost of a slight code duplication, this avoids create a new binary table at each call.

Parameters:

• ext_name (str) – The name of the extension that the column should be added to.

• columns (fits.Column instances) – The columns to be added.

average_deadtime_per_event()

Calculate the average deadtime per event.

Warning

Note this is trickier than it seems, as the result might not be accurate if the underlying event file has already been filteres in time and/or phase.

backscal()

Return the value of the BACKSCAL header keyword, if present.

The keyword is written out by xpselect when appropriate. The function returns None when the keyword is missing in the event file.

close()

Close the underlying HDU.

copy_and_filter(mask, history=None)

Filter the event file according to a given mask of selected rows.

This is similar in spirit to the filter() method, except that the filtering is not done in place, but everything is copied beforehand, so that the original HDU list is not modified, and multiple selections can be applied in series, at the expense of a larger memory footprint.

deadtime_correction()

Return the deadtime correction.

det_position_data()

Return the detx and dety column data, in either the detector x and y measured coordinates on the detector.

ds9_region_file_mask(file_path, mc=False)

Convenience function filtering a photon list from a ds9 region file.

ds9_region_mask(*region_list, mc=False)

Check which events are inside the logical or of all the regions within a given region list and return the corresponding array mask.

du_id()

Return the detector unit number used to run the simulation.

energy_data(mc=False)

Return the ra and dec column data, in either the Monte Carlo or the measured flavor depending on the keyword arguments passed to the binning class.

The try/except clause if needed for inter-operability with the ixpesim output files, where the column naming in the MONTE_CARLO extension is not in line with ixpeobssim, see https://bitbucket.org/ixpesw/gpdsw/issues/327

file_path()

Return the path to the underlying file.

filter(mask, history=None)

Filter the event file according to a given mask of selected rows.

This filtering happens in place, so it should be used when you don’t need to apply subsequent selections on the same file multiple times.

gti_data()

Return thr GTI data.

irf_name()

Return the name of the IRF set used to run the simulation.

Note the try/except block in the function body was added for ixpeobssim to inter-operate with the ixpesim output photon lists, which might have a MONTE_CARLO extension without the IRFNAME header keyword. In this case you have to set the keyword by hand, e.g., passing the proper command-line switch to xpbin.

livetime()

Return the livetime, i.e. the duration corrected for dead time.

livetime_data()

Return the LIVETIME column.

max_good_time()

Return the largest STOP time for the GTIs in the event file.

min_good_time()

Return the smaller START time for the GTIs in the event file.

num_events()

Return the total number of events in the event file.

octi_data()

Return the OCTI data.

phase_data()

Return the PHASE column.

phi_data()

Return the PHI column.

pi_data(mc=False)

Return the PI column.

pol_deg_weighted_average(pol_deg_model, ebins=10)

Calculate the weighted average of the polarization degree as a function of energy on the given event sample for a generic input polarization model.

This is achieved by calculating the input model of the polarization degree on an event by event basis (i.e., evaluating the model itself on the proper columns of the event file), creating a 2-dimensional histogram of the polarization degree vs. energy, and calculating the average of the polarization degree in vertical slices of energy.

Warning

The polarization degree is averaged coherently, i.e., under the assumption that the polarization angle is constant over all the event sample. Therefore using this function only makes sense when either the input polarization model has a constant polarization angle, or after the photoelectron directions in the event file have been properly aligned to the given input model.

Parameters:

• pol_deg_model (callable) – A Python function or an otherwise callable object with the proper signature returning the polarization degree as a function of energy, time and sky position (in this order).

• ebins (int or arra-like (default None)) – The energy binning for the weighted average. Loosely modeled on the numpy and scipy digitization functions, if bins is and integer, a logarithmically spaced binning between 2 and 8 keV is created.

Returns:

• An xScatterPlot object encoding the average polarization degree as a

• function of energy.

primary_keywords()

Return a list of all the relevant keywords in the PRIMARY header.

This is used, e.g., to propagate all the necessary information (such as the ROI and the IRFs used) from the event files to the binned data files that are created from them.

Warning

This function needs to be refactored!

q_data()

Return the Q column.

remove_columns(ext_name, *col_names)

Remove one or more columns from the specified extension of the event file.

Parameters:

• ext_name (str) – The name of the extension that the column should be added to.

• *col_names (str) – The name(s) of the column(s) to be removed

sc_data()

Return the spacecraft data.

set_column(ext_name, col_name, col_data)

Overwrite the data for an existing column.

Parameters:

• ext_name (str) – The name of the extension that the column should be added to.

• col_name (str) – The name of the column to be added.

• col_format (str) – The format for the new column (e.g., ‘E’).

• col_data (array_like) – The data for the new column.

sky_position_data(mc=False)

Return the ra and dec column data, in either the Monte Carlo or the measured flavor depending on the keyword arguments passed to the binning class.

Note that, since flight level-2 data do not contain RA and DEC columns, this uses internally X and Y.

srcid_data()

Return the SRC_ID column.

Note the try/except block in the function body was added for ixpeobssim to inter-operate with the ixpesim output photon lists, which come have a MONTE_CARLO extension without the SCR_ID column.

start_met()

Return the start MET of the observation.

stokes_data()

Return the Q and U columns.

stop_met()

Return the stop MET of the observation.

time_data()

Return the TIME column.

timeline_data()

Return the timeline data.

total_good_time()

Return the sum of all the GTIs in the event file.

trigger_id_data()

Return the TRG_ID column.

u_data()

Return the U column.

wcs_reference()

Return the reference point of the underlying WCS object.

Note that, since we are making sure that an event file always has a valid WCS associated with it (be it the actual content of the EVENTS header or one created on the fly given RA_SRC and DEC_SRC), this should be the preferred way to gauge the center of the field—compared, e.g., to using RA_SRC, DEC_SRC or RA_PNT, DEC_PNT directly.

write(file_path, overwrite=False)

Write the underlying HDU list to file.

write_fits_selected(selection_mask, file_path, header_keywords=None, history=None, overwrite=True, filter_in_place=True)

Write to file a subselection of events.

Parameters:

• selection_mask (array) – A mask for the row selection.

• file_path (string) – The path to the output file.

xy_data()

Return the X and Y columns.

class ixpeobssim.evt.event.xEventFileFriend(file_list2, file_list1=None)

Read/write interface to put together level-1 and level-2 event files.

It supports minimal output facilities, mainly re-implementing functions from xEventFile class

energy_data(mc=False)

Wrap xEventFile.energy_data() appending values for all the LV2 files

energy_l1_data()

Rewrite xEventFile.energy_data() to get energy in keV for L1 data

l1value(val, all_events=False)

l2value(val)

sky_position_data(mc=False)

Wrap xEventFile.sky_position_data() appending values for all the LV2 files

class ixpeobssim.evt.event.xEventList(time_=None, source_id=None)

Class describing an ixpeobssim event list.

apply_charging(irf_set, **kwargs)

Apply the GEM charging to the event list.

apply_dead_time(deadtime)

Apply the dead time effect to the event list.

apply_fiducial_area()

Trim out the events falling outside the detector active area.

apply_vignetting(vign, ra_pnt, dec_pnt)

Apply the effective area vignetting to the event list.

detector_coordinates()

Return the event position in detector coordinates.

energy()

Return the reconstructed event energies.

fill_livetime(gti_list, deadtime=0.0)

Fill the LIVETIME column.

The basic rule is that the livetime for the first event is the time elapsed since the start run, and for all the other events we just calculate the time difference and subtract the deadtime per event.

Parameters:

• gti_list (xGTIList instance) – The list of good time intervals for the observation.

• deadtime (float) – The deadtime per event.

fill_livetime_unphysical(value=-1)

Fill the LIVETIME column with an unphysical value.

This is called once when the event lists for the single source components are simulated in order to have the column in the event list itself filles and of the proper length. The actual trigger IDs are then set to the actual values after the source components are merged into the final event list and the deadtime is applied.

fill_trigger_id()

Fill the TRG_ID column with sequential values starting from 1.

fill_trigger_id_unphysical(value=-1)

Fill the TRG_ID column with an unphysical value.

This is called once when the event lists for the single source components are simulated in order to have the column in the event list itself filles and of the proper length. The actual trigger IDs are then set to the actual values after the source components are merged into the final event list and the deadtime is applied.

livetime()

Return the event livetimes.

mc_energy()

Return the Monte Carlo event energies.

mc_sky_coordinates()

Return the sky positions of the events.

phi()

Return the reconstructed azimuthal angle for the events.

pi()

Return the pulse invariant for the events.

q()

Return the Q Stokes parameter.

set_detector_position_columns(detx, dety)

Set all the columns related to the event position in the GPD frame.

set_energy_columns(energy, pha, pi)

Set all the columns related to the measured energy.

set_mc_energy_columns(mc_energy, mc_pha, mc_pi)

Set all the columns related to the Monte Carlo energy.

set_mc_gain_column(gain=1.0)

Set the column corresponding to the relative GEM gain.

set_mc_sky_position_columns(mc_ra, mc_dec, mc_x, mc_y)

Set all the columns related to the Monte Carlo sky position.

set_num_clusters(value=1)

Fill the NUM_CLU column with a fixed value.

Parameters:

• value (number or array) – The values for the column to be filled with.

• in (Note the control on whether value is a number or an array is performed)

• method. (the _set_column())

set_phi_columns(phi, detphi)

Set all the columns related to the photoelectron direction.

set_rec_columns(pha, pi, energy, ra, dec, x, y, detx, dety)

Set all the columns of the EVENT extension that imply a convolution of the MC truth with the instrument response functions.

set_seed_columns(mc_energy, mc_pha, mc_pi, mc_ra, mc_dec, mc_x, mc_y, phi, detphi)

Set all the columns for the seed event list.

These is the first set of columns being simulated in xpobssim, which is why we have this convenience function to set them all at once.

Note

Historically the columns set by this function included the time, but that was abandoned for the sake of consistency when we added the possibility of setting the event times at construction time. So the basic workflow is now:

• initialize the event list using right after the event times have been extracted passing the times themselves and the source identifier to the constructor;

• call set_seed_columns();

• call set_rec_columns().

set_sky_position_columns(ra, dec, x, y)

Set all the columns related to the measured sky position.

set_weights(value=1.0)

Fill the W_MOM column.

Parameters:

• value (number or array) – The values for the column to be filled with.

• in (Note the control on whether value is a number or an array is performed)

• method. (the _set_column())

stokes_parameters()

Return the Stokes parameters.

u()

Return the U Stokes parameter.

write_fits(creator, roi_model, irf_set, **kwargs)

Write the event list and associated ancillary information to file.

Basic data format definitions for the event lists.

ixpeobssim.evt.fmt.build_standard_wcs(ra0, dec0)

Build the standard WCS object for event files programmatically.

See https://bitbucket.org/ixpesw/ixpeobssim/issues/552

ixpeobssim.evt.fmt.set_object_header_keywords(hdu, ra_obj, dec_obj, name=None)

Set the object-related header keywords for a given HDU.

ixpeobssim.evt.fmt.set_standard_xy_header_limits(hdu)

Set the TLMIN and TLMAX keywords for the X and Y columns for a given hdu.

ixpeobssim.evt.fmt.set_telescope_header_keywords(hdu, du_id)

Set the telescope-related header keywords for a given HDU.

ixpeobssim.evt.fmt.set_time_header_keywords(hdu, start_met, stop_met, duration, ontime, livetime, deadc)

Set the time-related header keywords for a given HDU.

ixpeobssim.evt.fmt.set_version_keywords(hdu, version)

Set the version-related header keywords for a given HDU.

ixpeobssim.evt.fmt.set_wcs_header_keywords(hdu)

Set the WCS-related header keywords for a given HDU.

ixpeobssim.evt.fmt.standard_radec_to_xy(ra, dec, ra0, dec0)

Convert sky coordinates to X and Y digital coordinates in the sky frame.

This is builfding on the fly a standard IXPE WCS object centered at the given sky position, and using it for the conversion.

Parameters:

• ra (array_like) – The array of input Ra coordinates.

• dec (array_like) – The array of input Dec coordinates.

• ra0 (float) – The Ra coordinate of the center of the field in the sky.

• dec0 (float) – The Dec coordinate of the center of the field in the sky.

ixpeobssim.evt.fmt.standard_xy_columns_kwargs(ra0, dec0)

Return the appropriate keywords for the X and Y columns in event files.

See https://bitbucket.org/ixpesw/ixpeobssim/issues/552

ixpeobssim.evt.fmt.standard_xy_to_radec(x, y, ra0, dec0)

Convert X and Y digital coordinates to RA and DEC.

This is builfding on the fly a standard IXPE WCS object centered at the given sky position, and using it for the conversion.

Parameters:

• x (array_like) – The array of input X coordinates.

• y (array_like) – The array of input Y coordinates.

• ra0 (float) – The Ra coordinate of the center of the field in the sky.

• dec0 (float) – The Dec coordinate of the center of the field in the sky.

class ixpeobssim.evt.fmt.xBinTableHDUEvents(ra0, dec0, data=None, keywords=None, comments=None)

Binary table description for the EVENTS extension of the observation output files.

class ixpeobssim.evt.fmt.xBinTableHDUGTI(data=None, keywords=None, comments=None)

Binary table for the good time intervals (GTI).

class ixpeobssim.evt.fmt.xBinTableHDUMonteCarlo(data=None, keywords=None, comments=None)

Binary table description for the MONTE_CARLO extension of the observation output files, including the additional Monte Carlo fields.

set_irf_name(irf_name)

Set the IRFNAME keyword.

class ixpeobssim.evt.fmt.xBinTableHDUOCTI(data=None, keywords=None, comments=None)

Binary table for the on-orbit calibration time intervals (OCTIs).

set_cal_stats(num_runs, total_time)

Set the calibration statistics keywords.

class ixpeobssim.evt.fmt.xBinTableHDURoiTable(data=None, keywords=None, comments=None)

Binary table for the good time intervals (GTI).

set_center(ra0, dec0)

Set the keywords for the ROI center.

class ixpeobssim.evt.fmt.xBinTableHDUSpacecraftData(data=None, keywords=None, comments=None)

Binary table for the spacecraft data.

set_roll_angle(roll_angle)

Set the roll angle.

class ixpeobssim.evt.fmt.xBinTableHDUTimeline(data=None, keywords=None, comments=None)

Binary table for the timeline data.

class ixpeobssim.evt.fmt.xLvl2PrimaryHDU(data=None, header=None, creator='ixpeobssim', keywords=None, comments=None)

Level 2 primary header definition.

Data structures related to the good time intervals.

class ixpeobssim.evt.gti.xGTIList(start_met, stop_met, *gtis)

Small convenience class representing a list of good time interval.

The interfaces are fairly minimal, but since we use this in quite different places, it was handy to collect the useful stuff in one place.

all_mets()

Return all the MET values corresponding to the start or stop of the GTIs in the list.

Note that Python supports the call to sum() with start set by keyword argument only since version 3.8, so we are refraining from that, here.

append_gti(start, stop)

Append a new GTI to the list.

complement()

Return the logical complement of the GTI list in the relevant time span, i.e., the list of time intervals between the start and the end of observation that are not good time intervals.

These are the ones where we do not take celestial data and we can, e.g., take calibration data.

See https://stackoverflow.com/questions/16789776/ for the use of the iterator over the list.

filter_event_times(time_)

Filter a given array of event times and return a reduced array only containing the times within the good time intervals in the list, along with the corresponding boolean mask. The latter, in turn, can be used to filter ancillary related columns, e.g., the phase in simulations of periodic sources.

start_mets()

Return the value of the start MET values for the GTI.

stop_mets()

Return the value of the stop MET values for the GTI.

total_good_time()

Return the total good time.

class ixpeobssim.evt.gti.xSimpleGTIList(start_met, stop_met)

Subclass of xGTIList with a single GTI.

class ixpeobssim.evt.gti.xUberGTIList

Subclass of xGTIList with a single GTI including all the MET from -inf to inf.

Facilities to produce photon lists at the top of the window to be fed into ixpesim.

ixpeobssim.evt.ixpesim.build_tow_response(irf_set)

Build the response at the top of the Be window that is needed to generate the photon lists.

This is assembled “by hand” using the irfgen facilities, and to do this we need to know the proper files used to build the response functions in the first place. Unfortunately, for historical reasons, these are stored in the COMMENT field of the primary header of the FITS files, and therefore we have to resort to some string parsing to make this happen. Admittedly, this is fragile, and we are guaranteed to break it if we ever change the format of the comments. (Note, however, that we are starting with sensible defaults.)

class ixpeobssim.evt.ixpesim.xBinTableHDUPhotons(data=None, keywords=None, comments=None)

Binary table description for the PHOTONS extension in ixpesim event lists.

class ixpeobssim.evt.ixpesim.xPhotonList(time_=None, source_id=None)

Class describing a photon list.

This is, in many respects, the equivalent of the evt.event.xEventList class, except that it does not represent actual events in the detector, but photons at the top of the Be window.

apply_fiducial_area()

Trim out the events falling outside the detector active area.

apply_vignetting(vign, ra_pnt, dec_pnt)

Apply the effective area vignetting to the event list.

fill(energy, ra, dec, detx, dety, pol_deg, pol_ang)

Fill all the relevant columns of the event list.

write_fits(creator, roi_model, irf_set, **kwargs)

Write the photon list and associated ancillary information to file.

Polarization analysis as described in https://arxiv.org/pdf/1409.6214.pdf

class ixpeobssim.evt.kislat2015.xModulationAnalysis(phi, weights=None)

Small class implements the polarization analysis based on the event-by-event Stokes parameters described in Kislat et al. 2015, see https://arxiv.org/pdf/1409.6214.pdf

Parameters:

phi (array_like) – The array of photoelectron azimuthal directions.

calculate_modulation(degrees=False)

Calculate the modulation.

class ixpeobssim.evt.kislat2015.xStokesAnalysis(q, u, energy, modf, aeff, livetime, weights=None, acceptcorr=True)

This class implements the polarization analysis based on the event-by-event Stokes parameters described in Kislat et al. 2015, see https://arxiv.org/pdf/1409.6214.pdf

The basic idea is that we pass to the constructor a list of photoelectron direction and energy arrays, along with the necessary response functions (and optional weights), and the proper functional combination are turned into the correponding (weighted) event-by-event Stokes parameters, that can then be easily summed in the proper energy range and turned into polarization degree and error.

Note that (modulo a factor of 2) all the book-keeping and processing is done in terms of the reconstructed Stokes parameters defined at the end of section 3 of the paper, i.e., the Stokes parameters are divided by the modulation factor on an event-by-event basis in the constructor.

Parameters:

• q (array_like) – The array of event-by-event Q Stokes parameters.

• u (array_like) – The array of event-by-event U Stokes parameters.

• energy (array_like) – The array of the event energies (in KeV)—must have the same shape as q and u.

• modf (xModulationFactor instance) – The modulation factor—this is evaluated at the event energies and used to normalize q and u.

• aeff (xEffectiveArea instance) – The effective area—this is used for the energy flux calculation, and and to calculate the acceptance correction, if necessary.

• livetime (float) – The livetime (in s) for the energy flux calculation.

• weights (array_like) – Additional (optional) multiplicative event weights—must have the same shape as q and u.

• acceptcorr (bool) – If True, the Stokes parameters are weighted by the inverse of the acceptance.

W2(mask)

Return the sum of the squares of weights.

For an un-weighted analysis, and if the acceptance correction is not applied, the weights are all equal to unity, and this sum reduces to the number of events passing the cut, which is in turn equal to the sum of the I Stokes parameter.

average_energy(emin, emax)

Return the average energy over a given energy range.

static calculate_mdp99(mu, I, W2, clip=True)

Calculate the MDP based on equation (A.8) in the paper.

Parameters:

• mu (array_like) – The effective modulation factor.

• I (array_like) – The I Stokes parameter.

• W2 (array_like) – The sum of the weights squared.

• clip (bool) – If true, the MDP is clipped within the physical bounds 0–100%.

static calculate_n_eff(counts, I, W2)

Calculate the effective number of events.

static calculate_polarization(I, Q, U, mu, W2=None, degrees=False)

Calculate the polarization degree and angle, with the associated uncertainties, for a given q and u.

This implements equations (21), (36), (22) and (37) in the paper, respectively.

Note that the Stokes parameters passed as the input arguments are assumed to be normalized to the modulation factor (for Q and U) on an event-by-event basis and summed over the proper energy range.

Great part of the logic is meant to avoid runtime zero-division errors.

static calculate_polarization_sub(I, Q, U, I_ERR, Q_ERR, U_ERR, degrees=False)

Calculate the polarization degree and angle, along with the fellow associated uncertainties.

This is using the linear error propagation and is currently only used in the polarization cube subtraction.

static calculate_stokes_errors(I, Q, U, mu, W2)

Calculation of the errors on the Stokes parameters.

effective_mu(emin, emax)

Return the effective modulation factor weighted over the input events.

mdp_map_cube(x, y, binning)

Return the necessary cube to create a MDPMAPCUBE binned object.

Note the extensions names, here, are taken from the definition of the FITS file holding the data structure, and we have to be careful in passing out the actual values in the right order.

For reference, here is a snapshot the field, in the appropriate order: [‘E_MEAN’, ‘COUNTS’, ‘MU’, ‘W2’, ‘N_EFF’, ‘FRAC_W’, ‘MDP_99’, ‘I’]

property n

Simple property function to make the accumulation of the counts easy to read.

static normalized_stokes_parameters(I, Q, U, I_ERR, Q_ERR, U_ERR)

Return the normalized Stokes parameters QN = Q / I and UN = U /I properly handling possible zero-division error issues.

polarization(emin, emax, degrees=False)

Return the average polarization degree and angle, along with the corresponding statistical uncertainties, into a given energy range.

polarization_map_cube(x, y, binning)

Return the necessary cube to create a PMAPCUBE binned object.

polarization_table(ebinning, degrees=True)

Return a table with all the relevant polarization parameters.

Note the column names, here, are taken from the definition of the FITS file holding the data structure, and we have to be careful in passing out the actual values in the right order.

Warning

It would be nice to be able to grab the column definition from the proper file (binning.fmt) but for some reason I do get a circular import error, when I try and do that. It might be a sensible thing to understand why and do the proper refactoring.

static significance(Q, U, Q_ERR, U_ERR)

Calculate the significance of the polarization detection.

Here we take advantage of the fact that Q^2 / Var(Q) + U^2 / Var(U) is a chisquare with two degrees of freedom, a.k.a. an exponential, and we’re using the corresponding cumulative function.

The significance in gaussian equivalent sigma is then calculated with the ppf of a normal distribution.

static stokes_covariance(I, QN, UN, W2)

Calculate the covariance between Q and U.

static stokes_i(phi, weights=None)

Convert the event azimuthal angle to the I Stokes parameter, see equations (9a) and (A.1a) in Kislat et al., 2015.

Parameters:

• phi (array_like) – The array of azimuthal angles.

• weights (array_like) – Optional event weights.

static stokes_q(phi, weights=None)

Convert the event azimuthal angle to the Q Stokes parameter, see equations (9b) and (A.1b) in Kislat et al., 2015.

Note that, compared to equation (9b), we have an extra factor of 2, here, that renders the calculations downstream more natural. The rest of the class is implemented consistently.

This is factored out in a staticmethod so that it can be reused consistently in other places.

Parameters:

• phi (array_like) – The array of azimuthal angles.

• weights (array_like) – Optional event weights.

static stokes_u(phi, weights)

Convert the event azimuthal angle to the U Stokes parameter, see equations (9c) and (A.1c) in Kislat et al., 2015.

Note that, compared to equation (9c), we have an extra factor of 2, here, that renders the calculations downstream more natural. The rest of the class is implemented consistently.

This is factored out in a staticmethod so that it can be reused consistently in other places.

Parameters:

• phi (array_like) – The array of azimuthal angles.

• weights (array_like) – Optional event weights.

sum_stokes_parameters(emin, emax)

Sum the Stokes parameters into a given energy range.

Sonification utilities.

class ixpeobssim.evt.sonify.ContolChangeParameter

Most common contol change codes, see https://professionalcomposers.com/midi-cc-list/

class ixpeobssim.evt.sonify.ProgramChangePrograms

Program change programs, see https://www.recordingblogs.com/wiki/midi-program-change-message

ixpeobssim.evt.sonify.midi_to_wav(midi_file_path, sound_font='')

Convert a MIDI file to audio using fluidsynth.

ixpeobssim.evt.sonify.play_midi(midi_file_path, sound_font='')

Play a MIDI file file through fluidsynth.

ixpeobssim.evt.sonify.stereo_to_mono(file_path)

Convert a stereo wave file to a mono file containing the average of the two channels.

class ixpeobssim.evt.sonify.xMidiEvent(type_, timestamp, **kwargs)

Small classe describing a midi event.

This is intended as a temporary storage for actual MIDI messages expressed in absolute time and not necessarily time-ordered.

class ixpeobssim.evt.sonify.xMidiFile(ticks_per_beat=480)

Small wrapper around the midio.MidiFile class.

add_midi_track(track_name=None, bpm=None)

Add a track to the MIDI file.

static compute_note_number(energy, mode='Ionian', key='C', pivot=4.0, dynamic_scaling=1.0)

Convert an energy (in keV) into the corresponding note number.

This is accomplished assigninig the A4 note number to the pivot energy, and scaling the energies in octave space. Since the extended IXPE bandpass over which we calculate the response functions (1–15 keV) corresponds to roughly 4 octaves, this is not a terrible match for a sonification project. The nominal IXPE energy band (2–8 keV) is only 2 octaves, and in real life we might need to tweak things a little bit, which is the purpose of the dynamic_scaling argument (see the comments in the code for a detailed explanation of the algorithm).

Parameters:

• energy (array_like) – The photon energy.

• mode (str) – The name of the musical scale to which the energy values should be snapped.

• key (str) – The key for the aforementioned musical scale.

• pivot (float) – The pivot energy (mapped to A4).

• dynamic_scaling (float) – Empirical parameter for scaling the energy in octave spaces.

static compute_pan(phi, dynamic_scaling=1.0)

Compute the pan.

This is simply mapping the photoelectron angle into the 0–127 phisical range.

Parameters:

• phi (array_like) – The photoelectron angle.

• dynamic_scaling (float) – Empirical parameter to enhance the veocity dynamics.

static compute_velocity(mc_energy, rec_energy, pivot=64.0, dynamic_scaling=1.0)

Convert the energy measurement to note velocity.

The note velocity is determined by the ratio between the reconstructed and true energy, i.e., events in the left tail of the energy dispersion are rendered with a lower velocity.

Parameters:

• mc_energy (array_like) – The true energy.

• rec_energy (array_like) – The reconstructed energy.

• pivot_velocity (float) – The pivot value for the note velocity, assigned when the reconstructed energy is equal to the true energy.

• dynamic_scaling (float) – Empirical parameter to enhance the veocity dynamics.

static control_change_message(control, value, channel=0)

Return a control change message.

fill(file_path, roi=None, **kwargs)

Fill the MIDI file with the event data in the proper format.

This is the main function where most of the action actually happens.

static pan_message(value, channel=0)

Return a pan message.

static plot_diagnostics(timestamp, note, velocity, pan, duration)

Create all the diagnostics plots.

static program_change_message(program, channel=0)

Return a program change MIDI message.

static set_tempo_meta_message(bpm, time=0)

Return a set_tempo MIDI meta-message.

static track_name_meta_message(name, time=0)

Return a track_name MIDI meta-message.

class ixpeobssim.evt.sonify.xMidiNote(note_number)

Small class encapsulating a MIDI note.

There’s an infinite number of place on the web where one can find conversion tables between MIDI notes and physical characteristics, see, e.g., https://www.inspiredacoustics.com/en/MIDI_note_numbers_and_center_frequencies

A MIDI note is univiquely identified by a note number, ranging from 0 to 127. Here we shall restrict ourselves to the piano keys, i.e., note numbers from 21 to 127 (the lower notes are too low to be useful, anyway), with the understanding that:

• note 21 is A0, at 27.500 Hz;

• note 127 is G9, at 12543.854 Hz;

• note 69 is A4, at 440.000 Hz.

A MIDI note is initialized by its MIDI note number, and encapsulates all the facilities to calculate the frequency, note name and alike.

frequency()

Return the note frequency.

name()

Return the note name and octave.

note_name()

Return the note name.

octave()

Return the note octave.

class ixpeobssim.evt.sonify.xMusicalScale(mode='Ionian', key='C')

Simple class describing a musical scale.

A scale is essentilly given by a series (in the form of a numpy array) of note number covering the piano dynamic range.

snap(values)

Snap a series of notes to the scale.

Event list filtering facilities.

class ixpeobssim.evt.subselect.xEventSelect(file_path, **kwargs)

Base class for event subselection.

get(key, default=None)

Convenience method to address the keyword aguments.

mask_selected()

Return True if selecting with event mask directly,

phase_selected()

Return True if selecting in phase, i.e., phasemin and/or phasemax are not None.

select()

Select the events and write the output file.

set(key, value)

Convenience method to set keyword arguments.

time_selected()

Return True if selecting in time, i.e., tmin and/or tmax are not None.

ixpeobssim.evt.xspec_.add_model_string(key, value)

Add a key,value pair of strings to XSPEC’s internal database. (simple wrapper around the xspec.Xset.addModelString() method.)

This database provides a way to pass string values to certain model functions (e.g., POW_EMIN and POW_EMAX) which are hardcoded to search for “key”. (See the XSPEC manual description for the “xset” command for a table showing model/key usage.)

If the key,value pair already exists, it will be replaced with the new entries.

ixpeobssim.evt.xspec_.add_model_strings(**kwargs)

Facility to add multiple model string to xspec.

ixpeobssim.evt.xspec_.calculate_confidence_intervals(model_id=1, delta_fit_statistics=1.0)

Calculate the confidence intervals for a fit.

This is running the error XSPEC command, see https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node78.html

ixpeobssim.evt.xspec_.compare_fit_data(fit_data, target_values, threshold=5.0)

Compare the best-fit parameter values with the input model.

Parameters:

• fit_data (xXspecFitData instance) – The fit output.

• target_values (array_like) – The corresponding values from the input model.

• threshold (float) – The threshold (in units of sigma) for determining the agreement between the best fit parameters and the model inputs.

Returns:

• Zero if the bes-fit parameters agree with the input, and an error code > 0

• otherwise.

ixpeobssim.evt.xspec_.current_fit_output(model_id=1)

Return the current fit output.

ixpeobssim.evt.xspec_.energy_flux(emin=2.0, emax=8.0)

Return the value of the energy flux and of the photon flux between emin and emax.

ixpeobssim.evt.xspec_.fit(stat_method='chi', max_iterations=10, error=True)

Perform a spectral fit.

ixpeobssim.evt.xspec_.fix_parameter(par_index, par_value, model_id=1)

Fix a model parameter to a given value.

ixpeobssim.evt.xspec_.load_input_files(*file_list)

Read a set of PHA1 input files and create the corresponding xspec.Spectrum objects.

ixpeobssim.evt.xspec_.load_local_models(package_name='ixpeobssim')

Load all the ixpeobssim local models.

We introduced this function, and called it automatically, version 12.0.0, based on the idea that whoever had XSPEC installed also had the XSPEC local models installed. This turned out to be a glaring usability issue, as there are many possible ways one could be unable to compile the local models, see https://bitbucket.org/ixpesw/ixpeobssim/issues/350 As of version 12.1.0 we wrap the XSPEC call in a try / except, and happily proceed even if the local models cannot be loaded.

ixpeobssim.evt.xspec_.plot(figure_name='XSPEC spectrum', logy=True)

Custom, matplotlib-based implementation of the standard XSPEC spectral-fit plot.

ixpeobssim.evt.xspec_.plot_normalized_counts(plot_data, du_id, **kwargs)

Plot the normalized counts.

Note that we use full circles for positive values and empty circles for negative ones, so that we can plot Stokes parameters (that can go negative) in logarithmic scale.

ixpeobssim.evt.xspec_.plot_residuals(plot_data)

Plot the model residuals.

ixpeobssim.evt.xspec_.reset()

Global reset.

ixpeobssim.evt.xspec_.retrieve_plot_data(bin_alg, du_id)

Retrieve the plot data (and, if available, the corresponding model) for a given spectrum.

Warning

Note that the plot data within XSPEC are intented to be normalized wrt the energy—i.e., the input from the FITS table is divided by the width of the energy bin, and the units get an additional keV^{-1}.

ixpeobssim.evt.xspec_.sample_spectral_model(expression, parameters, emin=1.0, emax=12.0, num_points=250, name=None, source_id=1)

Sample the values for a generic XSPEC model, given an expression and a set of parameters.

Most notably, this is used to feed into ixpeobssim time-independent XSPEC spectral models.

Parameters:

• expression (str) – The model expression string, using full component names.

• name (str) – The model name.

• source_id (int) – The source number.

• parameters (dict or list) – The model parameters.

• emin (double) – Energy limits, in keV.

• emax (double) – Energy limits, in keV.

• num_points (int) – The number of points to sample the spectrum.

ixpeobssim.evt.xspec_.select_energy_range(emin=None, emax=None)

Apply a global mask on the energy value for the channels to be fitted.

ixpeobssim.evt.xspec_.set_parameter(par_index, par_value, model_id=1)

Set the model parameter value corresponding to provided index.

ixpeobssim.evt.xspec_.set_parameter_range(par_index, min_value, max_value, model_id=1)

Set the range for a parameter.

ixpeobssim.evt.xspec_.setup_fit_model(expression, par_values=None)

Create a model for spectral fitting in XSPEC.

Note we abort if the model cannot be created—this can happen, e.g., if one is trying to use local models without having compiled them. This seems as the only sane option, at this point, as the result of moving on would depend on the internal XSPEC state (is there any active model?) and this would be likely more confusing to the user than stating straight away what the problem is.

class ixpeobssim.evt.xspec_.xXspecFitData(model)

Small containter class to store the output of an XSPEC spectral fit.

This includes the parameter names, best-fit values and errors, the test statistics and the number of degrees of freedom.

next()

Alternative method for Python 2 compatibility.

par_error(identifier)

Return the parameter error at a given index.

par_value(identifier)

Return the parameter value at a given index.

stat_box(**kwargs)

Return a stat box that can be overlaid onto a given plot.

Note that when the error command has been run after the fit, we do include the average of the (asymmetric) error bars in the stat box.

class ixpeobssim.evt.xspec_.xXspecPlotData(energy, data, errors, model=None)

Small container class storing the XSPEC plot data.

Plot data are typically retrieved by function calls to the global xspec.Plot object, and this contained provides a convenient interface to cache plot data for later use. The container encapsulates the x and y values, the errors on the y values and, if available, the model calculated in correspondence of the x array.

save(file_path)

To be implemented.

class ixpeobssim.evt.xspec_.xXspecSpectrumManager

Small facility to keep track of the spectra being loaded into XSPEC.

This is essentially a dict where the type of the binned spectra and the DU Ids are mapped into the data group in the XSPEC memory, so that we can easily retrieve the plot data at any time.

data_group(bin_alg, du_id, default=None)

Return the XSPEC data group mapping to a particular binning algorithm and DU ID.

register(file_path, data_group)

Register a spectrum.

Note we peek into the binned file to get an hold on the binning algorithm and the DU ID.

spectrum_types()

Return a tuple with all the binning algorithms containing at least a spectrum.

Instrument

Charging-related facilities.

ixpeobssim.instrument.charging.asymptotic_gain_delta(energy_flux, k_c=120000.0, tau_d=110000.0, delta_max=0.067)

Return the asymptotic gain decrease for the GEM in an arbitrary state.

ixpeobssim.instrument.charging.calculate_gain(time_, energy_flux, initial_gain=1.0, k_c=120000.0, tau_d=110000.0, delta_max=0.067)

Calculate the gain as a function of time for a given input energy flux.

The time grid and the corresponding energy flux are passed as numpy arrays (of the same shape). The system is integrated using finite differences.

ixpeobssim.instrument.charging.charging_tau(energy_flux, k_c=120000.0)

Return the time constant for the GEM charging in absence of discharge.

ixpeobssim.instrument.charging.create_charging_map_extension(fast_map, slow_map=None, start_date='N/A', start_time='N/A', version=1)

Create the CHRG_MAP extension for a charging map file.

ixpeobssim.instrument.charging.effective_tau(energy_flux, k_c=120000.0, tau_d=110000.0)

Return the effective time constant for the GEM in an arbitrary state.

ixpeobssim.instrument.charging.gain_profile(energy_flux, k_c=120000.0, tau_d=110000.0)

Return a lambda to generate a gain profile, given some model parameters.

to plot the actual time-profile of the gain.

ixpeobssim.instrument.charging.read_charging_map(file_path)

Open a charging map file and read the values for the slow and fast component. We get the size of the arrays from the header.

ixpeobssim.instrument.charging.read_charging_parameters(input_file_path)

Open a charging parameters calibration file, read the parameters and return them.

Parameters:

input_file_path (string) – the path to the input FITS file storing the charging parameters

class ixpeobssim.instrument.charging.xBinTableHDUCharging(data=None, keywords=None, comments=None)

Binary table description for the CHRG_MAP extension of the observation output files.

class ixpeobssim.instrument.charging.xChargingPrimaryHDU(data=None, header=None, creator='ixpeobssim', keywords=None, comments=None)

Charging map file primary header.

class ixpeobssim.instrument.charging.xEnergyFluxCube(nside, tedges, zlabel='Mission Elapsed Time [s]', wlabel='Energy flux [keV mm$^{-2}$ s$^{-1}$]')

Main data structure for the charging-induced gain correction.

calculate_gain_data(k_c_fast=120000.0, tau_d_fast=110000.0, delta_max_fast=0.067, initial_dg_fast=None, k_c_slow=0.0, tau_d_slow=0.0, delta_max_slow=0.0, initial_dg_slow=None)

Calculate the expected space-time profile of the gain for a given illumination history. We are adopting a charging model with two components: one fast and one slow. Both are regulated by the same equation, which includes a rate-dependent charging and a constant discharging, with the caracteristic time for the two processes (fast and slow) separated by roughly an order of magnitude. The system is integrated using finite differences. The returned gain profile is computed on the times defined by self.tbinning() and has the same spatial dimensions as self.content

Parameters:

• k_c_fast (float) – the charging constant for the fast process (ADC counts/mm^2)

• tau_d_fast (float) – the discharge time constant for the fast process (seconds)

• delta_max_fast (float) – maximum fractional gain drop for the fast process

• initial_dg_fast (float or array) – initial fraction of gain drop due to the fast process. If an array, the value is given for each spatial bin - thus the dimension must match the spatial dimensions of self.content. If None, we assume initial full discharge for the fast component.

• k_c_slow (float) – as k_c_fast, but for the slow process

• tau_d_slow (float) – as tau_d_fast, but for the slow process

• delta_max_slow (float) – as delta_max_fast, but for the slow process

• initial_dg_slow (float) – as initial_dg_fast, but for the slow process

fill(detx, dety, time_, energy, deadtime_correction=None, **kwargs)

Overloaded fill() method.

Here we essentially take care of the energy weighting and the normalization by the elapsed time and the pixel area. This is not completely trivial as the mean energy must be calculated in the same time slices used to bin the original histogram, to cope with the possibility that the source spectrum—or the event flux on the detector, for whatever reason—changes with time.

Warning

The deadtime correction is passed as a single number, i.e., averaged over the full data set.

gain(detx, dety, time_, **charging_params)

Return the expected corrected gain for a list of events at specific detector positions and times.

gain_map(t)

Return a two-dimensional histogram with the gain map at a given time.

static step(flux, dt, dg, k_c, delta_max, tau_d)

Compute a step of the charging model (for either the slow or fast process). For a description of the charging parameters see the calculate_gain_data() method

Parameters:

• flux (float or array) – the value of the input energy flux for this step (ADC counts/mm^2/s)

• dt (float) – the time step (in seconds)

• dg (float or array) – the current gain drop due to this process (fractional)

tbinning()

Return the histogram binning on the time axis (aka z-axis).

time_slice(tmin=None, tmax=None)

Return a time slice of the irradiation history, between generic bounds.

By default this is returning the average energy flux per unit area over the entire input sample.

Note that the times are adjusted to the closest bih edges.

xbinning()

Return the histogram binning on the detx axis.

ybinning()

Return the histogram binning on the dety axis.

ixpeobssim.instrument.du.det_name_to_du_id(det_name)

Convert the logical name of a given DU to the actual DU ID.

ixpeobssim.instrument.du.du_logical_name(du_id)

Return the logical name for a given DU (i.e., the content of the DETNAME header keyword).

ixpeobssim.instrument.du.du_physical_name(du_id)

Return the physical name for a given DU (i.e., the content of the DET_ID header keyword).

ixpeobssim.instrument.du.du_rotation_angle(du_id, roll_angle=0.0, modulo=True)

Return the rotation angle for a given DU (modulo 2 pi).

Basic GPD-related constants.

ixpeobssim.instrument.gpd.detphi_to_phi(detphi, du_id, roll_angle)

Convert the phi angle from the GPD reference frame to the focal plane (i.e., sky) reference.

Parameters:

• detphi (float or array) – The azimuthal angle in the GPD reference frale in radians

• du_id (int) – The detector unit hosting the GPD

• roll_angle (float) – The telescope roll angle in decimal degrees

ixpeobssim.instrument.gpd.fiducial_area(half_side_x=6.6, half_side_y=6.8)

Return the area of the fiducial rectangle.

This is essentially calling the base rectangle_area() function, with the proper default values for the fiducial cut.

Parameters:

• half_side_x (float) – The half side of the fiducial rectangle along the x coordinate in mm.

• half_side_y (float) – The half side of the fiducial rectangle along the y coordinate in mm.

ixpeobssim.instrument.gpd.gpd_map_binning(half_side_x, half_side_y, num_bins_x, num_bins_y=None)

Return a numpy array with an appropriate binning for a bidimensional map over the GPD area.

Note this function was refactored to adapt to a generic rectangle in response to https://github.com/lucabaldini/ixpeobssim/issues/668, so this now returns two independent arrays representing the binng on the x and y axes.

Parameters:

• half_side_x (float) – The half side of the histogram binning along the x coordinate (in mm).

• half_side_y (float) – The half side of the histogram binning along the y coordinate (in mm).

ixpeobssim.instrument.gpd.phi_to_detphi(phi, du_id, roll_angle)

Convert the azimuthal angle from the focal plane (i.e., sky) reference frame to the GPD reference.

Parameters:

• phi (float or array) – The azimuthal angle in the sky reference frame in radians

• du_id (int) – The detector unit hosting the GPD

• roll_angle (float) – The telescope roll angle in decimal degrees

ixpeobssim.instrument.gpd.rectangle_area(half_side_x, half_side_y)

Small convenience function to calulate the area of a rectangle, given the half sides along the two coordinates.

Parameters:

• half_side_x (float) – The half side of the rectangle along the x coordinate (in mm).

• half_side_y (float) – The half side of the rectangle along the y coordinate (in mm).

ixpeobssim.instrument.gpd.rotate_detxy(x, y, du_id, roll_angle=0.0, inverse=False)

Convert the (x, y) position from the focal plane reference frame to the the gpd reference frame based on the DU id and roll angle.

Warning

This would be more consistent with the rest of the code base if the direct and inverse rotations were implemented in two different functions rather than one function with an additional argument.

Parameters:

• x (float or array) – The x position or array of x positions in mm

• y (float or array) – The y position or array of y positions in mm

• du_id (int) – The detector unit hosting the GPD

• roll_angle (float) – The telescope roll angle in decimal degrees

ixpeobssim.instrument.gpd.within_fiducial_rectangle(x, y, half_side_x=6.6, half_side_y=6.8)

Return wheter an (x, y) position in mm is within a fiducial rectangle of given half-sides on the two coordinates.

Parameters:

• x (float or array) – The x position or array of x positions in mm

• y (float or array) – The y position or array of y positions in mm

• half_side_x (float) – The half side of the fiducial rectangle along the x coordinate (in mm).

• half_side_y (float) – The half side of the fiducial rectangle along the y coordinate (in mm).

ixpeobssim.instrument.gpd.within_gpd_physical_area(x, y)

Return wheter an (x, y) position in mm is within the GPD physical area.

Parameters:

• x (float or array) – The x position or array of x positions in mm

• y (float or array) – The y position or array of y positions in mm

ixpeobssim.instrument.mma.apply_dithering(time_, ra_pnt, dec_pnt, dither_params=None)

Apply dithering correction directly to pointing direction in celestial coordinates. In this way the correction gets propagated down the chain and the appropriate IRFs are used. See convolve_event_list in mma.py

Parameters:

• time (float or array) – The event time

• ra_pnt (float or array) – The pointing ra in decimal degrees

• dec_pnt (float or array) – The pointing dec in decimal degrees

• dither_params ((amplitude, pa, px, py) tuple, optional) – The parameters for the dithering of the observatory. The dithering is not applied if this is set to None.

ixpeobssim.instrument.mma.fiducial_backscal(half_side_x, half_side_y)

Calculate the BACKSCALE value for a give fiducial rectangle in detector coordinates.

This function was introduced in response to the issue https://github.com/lucabaldini/ixpeobssim/issues/668 in place of the old FIDUCIAL_BACKSCAL constant.

Parameters:

• half_side_x (float) – The half side of the fiducial rectangle along the x coordinate in mm.

• half_side_y (float) – The half side of the fiducial rectangle along the y coordinate in mm.

ixpeobssim.instrument.mma.gpd_to_sky(detx, dety, time_, ra_pnt, dec_pnt, du_id, roll_angle, dither_params=None)

Convert an array of (detx, dety) positions in detector coordinates to an array of (ra, dec) positions in the sky. This is supposed to be the inverse of sky_to_gpd(), so that the two operations, applied in series, should preserve the original data.

Parameters:

• detx (float or array) – The detx position in detector coordinates in mm

• dety (float or array) – The dety position in detector coordinates in mm

• time (float or array) – The event time

• ra_pnt (float or array) – The pointing ra in decimal degrees

• dec_pnt (float or array) – The pointing dec in decimal degrees

• du_id (int) – The detector unit id

• roll_angle (float) – The telescope roll angle in decimal degrees

• dither_params ((amplitude, pa, px, py) tuple, optional) – The parameters for the dithering of the observatory. The dithering is not applied if this is set to None.

ixpeobssim.instrument.mma.parse_dithering_kwargs(**kwargs)

Parse the keyword arguments related to the dithering.

ixpeobssim.instrument.mma.sky_to_gpd(ra, dec, time_, ra_pnt, dec_pnt, du_id, roll_angle, dither_params=None)

Convert an array of (ra, dec) positions in the sky to an array of (x, y) positions onto the gpd reference frame. This involves essentially three different steps: first we project the sky coordinates onto the focal plane reference frame, then we apply the dithering pattern (if enabled) based on the event times and finally we rotate the coordinates according to the DU id and the telescope roll angle.

Warning

Change x and y to detx and dety!

Parameters:

• ra (float or array) – The ra position in the sky in decimal degrees

• dec (float or array) – The dec position in the sky in decimal degrees

• time (float or array) – The event time

• ra_pnt (float or array) – The pointing ra in decimal degrees

• dec_pnt (float or array) – The pointing dec in decimal degrees

• du_id (int) – The detector unit id

• roll_angle (float) – The telescope roll angle in decimal degrees

• dither_params ((amplitude, pa, px, py) tuple, optional) – The parameters for the dithering of the observatory. The dithering is not applied if this is set to None.

Spacecraft-related facilities.

ixpeobssim.instrument.sc.dithering_pattern(amplitude=1.6, period_a=907.0, period_x=101.0, period_y=449.0)

Implementation of the full dithering pattern as per the report by Allyn Tennant and Kurt Dietz linked in the issue page: https://bitbucket.org/ixpesw/ixpeobssim/issues/193

Note this returns an anonymous function that can be evaluated into a generic array of time values, the basic usage being:

>>> t = numpy.linspace(0., 10000., 10001)
>>> dithering = dithering_pattern()
>>> x, y = dithering(t)

Parameters:

• amplitude (float) – The dithering amplitude in arcminutes.

• period_a (float) – The main dither period.

• period_x (float) – The x dithering period.

• period_y (float) – The y dithering period.

ixpeobssim.instrument.sc.period_to_omega(period)

Convert a characteristic period (in s) into the corresponding omega (in rad/s).

ixpeobssim.instrument.sc.pointing_direction(sc_data, met)

Return the pointing direction at a given array of MET.

ixpeobssim.instrument.sc.pointing_splines(sc_data)

Build a pair of R.A. and Dec dithering splines starting from as set of spacecraft data.

This can be used to recover the actual pointing as a function of time, given a SC_DATA binary table.

ixpeobssim.instrument.sc.pow_triangular_wave(x, amplitude, period, exponent=0.5)

Modified triangual wave, where the relative values are raised to a generic exponent, in such a way that the values of the maxima and minima are preserved.

ixpeobssim.instrument.sc.spiral_dithering_pattern(amplitude=1.6, period_theta=100.0, period_r=970.0, exponent=0.5)

Alternative, spiral-like dithering pattern.

ixpeobssim.instrument.sc.triangular_wave(x, amplitude, period)

Basic description of a (symmetric) triangular wave with a given period and amplitude.

The basic expression is taken from https://en.wikipedia.org/wiki/Triangular_distribution and, in this form, the function evaluates to 0 for x = 0.

Spacecraft trajecory related functions

class ixpeobssim.instrument.traj.xIXPETrajectory(ephemeris_file_name='de430t_2000_2040.bsp', min_sun_angle=65.0, max_sun_angle=115.0, builtin_timescale=True)

A class to characterize the spacecraft’s orbit and trajectory for ixpeobssim.

This is supposed to keep track of the necessary elements to keep track of the state of the observatory and generate the GTI.

The TLE data are read from the proper file in $IXPEOBSSIM_INSTRUMENT_DATA (which satellite precisely is controlled by the satellite_name argument).

The planetary and lunar ephemeris are taken from JPL. While JPL distributes and maintain an extensive set of files, in order to avoid large downloads at runtime, we have packaged or own trimmed versions of DE405 and DE430t within the 2000–2040 time span into the ixpeobssim/instrument/data folder. The trimming has been performed using the jplephem utility described at https://rhodesmill.org/skyfield/planets.html

The generation of the timescale information takes advantage of all the skyfield machinery. Note that, since the default skyfield method do not reasily offer any flexibility in terms of where the ancillary files are located, we have implemented our own Loader class that behaves as the the original with two notable exceptions:

• the builtin timescale ancillary files (when necessary) are loaded from IXPEOBSSIM_INSTRUMENT_DATA, rather than from the sjyfield installation;

• the updated timescale ancillary files (when necessary) are downloaded (and read from) IXPEOBSSIM_DATA.

Parameters:

• ephemeris_file_name (str) – The ephemeris file name (in the IXPEOBSSIM_INSTRUMENT_DATA folder).

• builtin_timescale (bool) – Flag passed to the skyfield timescale generation routine (if True the static files shipped with ixpeobssim are used, otehrwise they are downloaded).

static complement_epochs(epochs, start_met, stop_met)

Return the complement to a list of epochs.

earth_occultation_epochs(start_met, stop_met, target_ra, target_dec, initial_step=100.0, precision=0.001)

Calculate the epochs when a given target is occulted by the Earth.

earth_occultation_mets(start_met, stop_met, target_ra, target_dec, initial_step=100.0, precision=0.001)

Return the list of METS for which the Earth occultation status is changing.

elevation(met)

Return the altitude of the satellite (in km) at a given MET.

static equispaced_met_grid(start_met, stop_met, approximate_step)

Return an equally-spaced grid of MET values according to the input argument.

Note that the array returned ranges exactly from start_met to stop_met, and the step size is rounded to guarantee that.

geocentric(met)

Return the geocentric position at a given MET.

gti_list(start_met, stop_met, target_ra, target_dec, saa=True, occult=True, initial_step=100.0, precision=0.001)

Calculate the GTIs, i.e., the ephocs where we are not in the SAA and the target is not occulted by the Earth.

While this method seems (and is, in fact) quite convoluted, we decided to calculate separately the SAA and Earth occultation epochs, when the instrument is not taking data, and merge them at the end to extract the good time intervals. The reasoning is that, while the former epochs are quite regular over the timescale of an observation—which allows to perform a sensible binary search starting from a reasonably corse grid (say with a 100 s step)—when we put in AND the two we get all kind of resonances, to the point that the GTIs can be arbitrarily small, and we would need a much finer (and computationally intensive) initial grid in order not to introduce measurable errors.

Parameters:

• start_met (float) – The initial MET for the observation.

• stop_met (float) – The final MET for the observation.

• target_ra (float) – The right ascension of the target source in decimal degrees.

• target_dec (float) – The declination of the target source in decimal degrees.

• saa (bool) – Consider the SAA passages in the GTI calculation.

• occult (bool) – Consider the Earth occultation in the GTI calculation.

• initial_step (float) – The step (in s) for the initial equispacedtime time grid used to seed the binary search. (This can slow things down when too small.)

• precision (float) – The worst-case precision of the GTI calculation (i.e., the stop condition for the binary search)

in_saa(met)

Return whether the satellite is in the SAA at a given time.

static load_tle(tle_file_name='tle_data_science.txt', satellite_name='AGILE')

Load TLE data for a given satellite. Note the tle() method is deprecated with an entertaining comment DEPRECATED: in a misguided attempt to be overly convenient, this routine builds an unweildy dictionary of satellites with keys of two different Python types: integer keys for satellite numbers, and string keys for satellite names. It even lists satellites like ISS (ZARYA) twice, in case the user wants to look them up by a single name like ZARYA. What a mess. Users should instead call the simple tle_file() method, and themselves build any dictionaries they need.

Parameters:

• tle_file_name (str) – The name of the file (in the IXPEOBSSIM_INSTRUMENT_DATA folder) containing the TLE data.

• satellite_name (str) – The name of the satellite to be looked for in the TLE data ancillary files. In absence of sensible TLE data for IXPE, NUSTAR and AGILE are two scientific satellites in near-equatorial orbit that can be used for our purpose.

met_to_ut(met, use_iau2000b=True)

Convert MET(s) to skyfield.Time object(s), using the proper underlying timescale.

See https://github.com/skyfielders/python-skyfield/issues/373 for all the gory details of the iau2000b magic! That really saved the day.

position(met)

Return the position (longitude and latitude in decimal degrees) of the satellite at a given MET.

saa_epochs(start_met, stop_met, initial_step=100.0, precision=0.001)

Calculate the SAA epochs.

saa_mets(start_met, stop_met, initial_step=100.0, precision=0.001)

Return the list of METs for which we cross the SAA boundaries within a given time interval.

sun_constrained(met, target_ra, target_dec)

Return whether a given target at a given time cannot be observed due to the Sun constraints.

sun_constraint_epochs(start_met, stop_met, target_ra, target_dec, initial_step=1000.0, precision=0.001)

Calculate the epochs when observations are inhibited by the Sun constraints.

sun_constraint_mets(start_met, stop_met, target_ra, target_dec, initial_step=1000.0, precision=0.001)

Calculate the list of METs for which the Sun observability changes.

target_moon_angle(met, target_ra, target_dec)

Return the angular separation between a target RA and DEC and the Moon at a given MET.

target_occulted(met, target_ra, target_dec, limiting_altitude=200.0)

Return whether a given target is occulted at a given time.

target_planet_angle(met, target_ra, target_dec, planet_name)

Return the angular separation between a target RA and DEC and a planet at a given MET time.

Warning

For now the angular separation is calculated from a geocentric (i.e., not spacecraft) position, which should be sufficient for our needs.

target_sun_angle(met, target_ra, target_dec)

Return the angular separation between a target RA and DEC and the Sun at a given MET.

timeline_mets(start_met, stop_met, target_ra, target_dec, saa, occult, initial_step=100.0, precision=0.001)

Return a list of all the relevant MET values for the calculation of the observation timeline, i.e., those that pertain to the SAA and to the Earth occultation.

The calculation is started with the observaton start and stop met, and the SAA and Earth occultation MET values are addedd if necessary.

The function returns separate lists for the SAA and Earth occultation MET values, as well as a sorted logical OR of the two with no, diplicates, since all is needed downstream to calculate the timeline.

class ixpeobssim.instrument.traj.xObservationTimeline(start_met, stop_met, target_ra, target_dec, saa, occult)

Small container class encapsulating the detailed planning, in terms of GTIs, calibration intervals and SAA passages, for a specific observation (i.e., a portion of the mission with a specific, fixed pointing).

This is complete rewrite of the code that used to live into the gti_list() method of the xIXPETrajectory class, ini response to https://bitbucket.org/ixpesw/ixpeobssim/issues/417 The xIXPETrajectory class is still responsible for all of the complex logics and calculation to detect the relevant MET values for the SAA passages and the Earth occultation, and these values are turned into a list of xTimelineEpoch objects to be consumed downstream.

This class provides a coherent and unified interface to extract good time intervals and onboard calibration intervals, with optional padding on either end and with the possibility of enforcing a minimum duration (after the padding).

epoch_data()

Return the epoch data in a form that can be used to fill the TIMELINE extension in the output FITS files.

filter_epochs(min_duration=0.0, start_padding=0.0, stop_padding=0.0)

Filter the epochs by duration.

gti_list(min_duration=0.0, start_padding=0.0, stop_padding=0.0)

Return the list of the good time intervals (GTI).

The GTIs are defined as those when we are not in the SAA and we are not occulted by the Earth.

Note the return value of this function is not a plain old list, but a xGTIList object.

octi_list(min_duration=0.0, start_padding=0.0, stop_padding=0.0)

Return the list of the onboard calibration time intervals (OCTI).

The OCTIs are defined as those when we are occulted by the Earth, but outside the SAA.

sc_data(time_step, dither_params=None, saa=True, occult=True)

Return the spacecraft data to be written into the SC_DATA optional extension of the output photon list.

Spacecraft data are instantaneous data (e.g., longitude, latitude, elevation) evaluated on a regular grid with constant spacing. Rather than with the start or the stop of the observation, the grid is aligned with MET = 0, and padded as necessary on both ends to include the entire span of the observation.

The relevant data are returned into the form of a Python dictionary of the form {column_name: value_array} that can be used directly to fill the columns of the output binary table.

Arguments:

time_step: float

Interval used to save the SC information

dither_params(amplitude, pa, px, py) tuple, optional

The parameters for the dithering of the observatory. The dithering is not applied if this is set to None.

saabool

Flag indicating whether the time intervals in the SAA should be flagged—the corresponding column with be identically zero if this is false.

occultbool

Flag indicating whether the time intervals when the target is occulted should be flagged—the corresponding column with be identically zero if this is false.

class ixpeobssim.instrument.traj.xSAABoundary(file_name='saa_polygon.txt')

Small container class representing the SAA boundary.

The main purpose of this class is to encapsulate the algorithm to figure out whether any given point on the Earth surface is inside the SAA polygon or not. You can glance through the related issue at https://bitbucket.org/ixpesw/ixpeobssim/issues/232 to get an account of the underlying discussion, but essentially we opted for the SAA class to be a subclass of matplotlib.path.Path since the latter provides all the functionality that we need (includeing plotting), implemented efficiently and without bringing in yet another dependence.

contains(lon, lat)

Return a boolean mask signaling the pairs of coordinates inside the polyogon.

plot(plot_vertices=True, **kwargs)

Plot the polygon.

class ixpeobssim.instrument.traj.xSkyfieldLoader(verbose=True, expire=True)

timescale(delta_t=None, builtin=False)

Open or download the three timescale files, returning a Timescale.

This is an almost verbatim re-implementation of the timescale() method in skyfield/iokit.py, the only difference being the locations from where we retrieve the data.

class ixpeobssim.instrument.traj.xTLE

Small utility function to create TLE strings for the purpose of our simulation.

See https://en.wikipedia.org/wiki/Two-line_element_set for some background information.

Mind this is not meant to be a complete interface, which would be silly—what we are after is really the ability to generate a fictional TLE for a simple circular orbit with arbitrary altitude and inclination. We purposedly do not provide the ability to write ouput file in order not to pollute the environment—this is something that should remain internal to the simulation.

And, for completeness, here is the IXPE TLE pulled out from https://www.n2yo.com/satellite/?s=49954 shortly after launch.

1 49954U 21121A 21351.00640149 .00001120 00000-0 35770-4 0 9994 2 49954 0.2300 281.7657 0011347 134.4260 303.9164 14.90740926 1166

static lines(inclination=0.23, altitude=601.1, catalog_no=49954, launch_year=2021, launch_no=121, launch_piece='A', launch_day=351.00640149, ballistic_coefficient='.00001120', drag_term='35770-4', element_number=9994, ascending_node=281.7657, eccentricity='0011347', perigee=134.426, mean_anomaly=303.9164)

class ixpeobssim.instrument.traj.xTimelineEpoch(start_met, stop_met, in_saa, occulted)

Small convenience class to encapsulate an “epoch” within a given observation.

In this context an epoch is the datum of a start_met, a stop_met, and a series of bit flags indicating whether the observatory is, e.g., in the SAA or occulted by the Earth.

See https://bitbucket.org/ixpesw/ixpeobssim/issues/417

bitmask()

Return a bit mask corresponding to the observatory status during the epoch.

isgti()

Return True if the epoch is a good time interval, i.e., we’re not in the SAA and we’re not occulted by the Earth.

isocti()

Return True of the epoch is a good on-orbit calibration time interval, i.e., we are occulted by the Earth but not in the SAA.

shrink(start_padding, stop_padding)

Create a new epoch where the bounds are shrinked by given amounts on either side. The “shrinked” epoch is inhering the SAA and occultation flags.

Note that the original object is left untouched, i.e., the new epoch is created in place.

Response Functions

Read interface to effective area files.

class ixpeobssim.irf.arf.xEffectiveArea(file_path)

VALID_WEIGHTING_SCHEMES = (None, 'UNWEIGHTED', 'NEFF', 'SIMPLE')

Class describing the on-axis effective area.

Parameters:

file_path (str) – The path to the .arf FITS file containing the effective area table.

plot()

Plot the effective area.

weighting_scheme()

Return the weighting scheme used to calculate the effective area.

See https://bitbucket.org/ixpesw/ixpeobssim/issues/573/ for more information about why we do need to keep track of this.

The basic purpose that this hook serves is to be able to make sure that weights are used consistently across binned data products and response files, e.g., we don’t want for people to be able to create weighted polarization cubes with the effective area correction unless the arf file is created with the ‘SIMPLE’ prescription.

For backward compatibility this is returning None if the response file is missing the STOKESWT keyword in the SPECRESP extension. If this is the case, the user should not make assumptions about the actual weighting scheme used.

class ixpeobssim.irf.arf.xTowEffectiveArea(file_path)

Class describing the on-axis effective area at the top of the window.

Parameters:

file_path (str) – The path to the FITS file containing the effective area table.

Base classes for all the read interfaces to response files.

The classes in this module represent the base of the hyerarchy for all the classes descibing response files, and particularly:

• xResponseBase acts as a base class for all the response files;

• xSpecRespBase acts as a base class for the (on-axis) effective area, the modulation factor and the modulation respose function.

The latter (i.e., the classes inheriting from xSpecRespBase) all have in common being one-dimensional functions of the (true) energy, and have a SPECRESP extension containing the response in bins of energy.

The xSpecRespBase class is also equipped to allow for systematic errors on the response values, according to the proper OGIP standard.

class ixpeobssim.irf.base.xResponseBase(file_path, extension)

Base class for reading response data from FITS files.

Parameters:

• file_path (str) – The path to the input FITS file.

• extension (str) – The extension of the input FITS file (typically fits, arf or rmf)

field(col_name, ext_index=1)

Return a view over a column of the specified extension as an array.

Since most of the response files have exactly one relevant (binary table) extension, having the extension index defaulting to one is a handy shortcut to retrieve the column data in a general fashion.

Parameters:

• col_name (str) – The name of the target column in the binary table extension.

• ext_index (int) – The index of the target binary table extension.

header_comments()

Return the content of the COMMENT keyword of the primary header, split into lines.

This is essentially only used in ixpeobssim.evt.ixpesim to reverse-engineer the value of the GPD pressure for the purpose of building the telescope response at the top of the Be window.

class ixpeobssim.irf.base.xSpecRespBase(file_path, extension, k=2, pad=True)

Derived class describing a spectral response, i.e., effective area, modulation response function, or modulation factor.

has_sys_errors()

Return True if the object instance has both negative and positive systematic errors defined.

plot_base(**kwargs)

Plot the spectral response.

plot_sys_envelope()

Plot the envelope of the systematic errors.

plot_sys_errors()

Plot the fractional systematic errors.

sys_envelope()

Return the envelope of the systematic uncertainties.

Basic IRF naming conventions and CALDB-related facilities.

ixpeobssim.irf.caldb.irf_file_name(base, du_id, irf_type, intent, version, simple_weighting=False, gray_filter=False)

Return the file name for a specific response function.

Parameters:

• base (str) – Base name of the set of response functions.

• du_id (int) – Identifier of the specific DU.

• irf_type (str) – Identifier for the calibration data.

• intent (str) – The specific intent for the response file.

• version (int) – Version number.

• simple_weighting (bool) – If True, load the response file with the SIMPLE weighting scheme.

• gray_filter (bool) – If True, load the response files apppropriate for the gray filtes (where) this makes sense.

• is (The basic naming convention scheme used for the CALDB by HEASARC)

• ixpe_[instrument]_[datatype]_[date]_[ver].[ext]

• where

• (d1/d2/d3); (* [instrument] indicates the detector)

• data; (* [datatype] provides an identifier for the calibration)

• file; (* [date] indicates the first date of validity of the)

• file. (* [ver] is the CALDB version number for that)

• conventions (We approximately follow the HEASARC)

• that (except for the fact)

• date. (we don't really have a concept of validity)

• filter (Note that all the logic for handling the simple weighting and the gray)

• https (was moved here following)

ixpeobssim.irf.caldb.irf_file_path(irf_name, du_id, irf_type, caldb_path=None, check_file=True, simple_weighting=False, gray_filter=False)

Return the full file path to a particular IRF file.

ixpeobssim.irf.caldb.irf_folder_path(irf_type, caldb_path='/home/docs/checkouts/readthedocs.org/user_builds/ixpeobssim/checkouts/latest/ixpeobssim/caldb')

Return the CALDB folder for a particular IRF type.

ixpeobssim.irf.caldb.parse_irf_name(irf_name, delimiter=':')

Parse a generic IRF name and return the basic field values it encapsulates.

The irf_name is an internal designation that ixpeobssim uses to identitify a full (self-consistent) set of response functions.

Modulation-related facilities.

class ixpeobssim.irf.modf.xAzimuthalResponseGenerator

Random number generator for the azimuthal response of the polarimeter.

Here is the basic underlying math. Typically the response of a polarimeter to monochromatic, linearly polarized incident radiation is of the form:

This can be conveniently rewritten in term of the overall normalization (i.e., the total number of events) and the modulation, defined as

(the modulation can also be characterized as the fraction of modulated events in a given distribution, as can be readily demonstrated, and it is the quantity that the detector is effectively measuring). The angular response can be rewritten as

For completeness, the modulation, the modulation factor and the polarization degree for a monocromatic source are related to each other through:

(i.e., at any given energy the modulation factor is the modulation of the detector response for 100% linearly polarized incident radiation).

In terms of throwing random numbers, the phase is a trivial constant that can be added after the fact (modulo 2pi), so effectively the relevant probability density function is

The corresponding cumulative distribution function is

and it is unfortunate that it cannot be inverted, otherwise we would have no need to interpolate for generating random numbers according to this distribution.

From the prospecive of the code, this generator is a standard xUnivariateAuxGenerator where the azimuthal angle is our random variable and the modulation is our auxiliary variable. For any given value of the modulation, a vertical slice is providing the corresponding one-dimensional pdf.

The class also provide facilities to fit a histogram to recover the underlying modulation and phase.

Example

>>> import numpy
>>> from ixpeobssim.utils.matplotlib_ import plt
>>> from ixpeobssim.irf.modf import xAzimuthalResponseGenerator
>>>
>>> generator = xAzimuthalResponseGenerator()
>>> modulation = numpy.full(1000000, 0.5)
>>> phase = numpy.radians(45.)
>>> phi = generator.rvs_phi(phase, modulation)
>>> hist = plt.hist(phi, bins=numpy.linspace(0, 2*numpy.pi, 100),
                    histtype='step', label='Random angles')
>>> fit_results = generator.fit_histogram(hist)
>>> fit_results.plot()
>>> plt.show()

build_horizontal_ppf()

Overloaded method.

This is essentially done for a few reasons:

• we do have an analytical form of the cdf that we can use to reduce the numerical noise;

• this random generator is peculiar in that the pdf is linear in the auxiliary variable m;

• this is possibly the most important random engine in the framework and its accuracy must be tested carefully.

There are many parameters that can be optimized, here, in order to maximize the accuracy of the ppf representation. (This accuracy is characterized in test/test_azimuthal_response.py in terms of standard and maximum relative error). Among these are the grids sampling the modulation and quantile values, and the order of the various splines involved. We put some thoughts into getting this right, but we cannot exclude that this can be improved.

For the modulation values we just take a constant-pitch grid, while for the quantiles we calculate the cdf (for m = 1, which is where the deviations from a straight lines are largest) on a constant-pitch grid in angle.

We empirically found that k=3 is the best spline index for the ppf at any given modulation value, while we get the maximum accuracy for kx = ky = 1 in the actual bivariate spline representing the ppf.

static cdf(phi, m)

Evaluate the underlying one-dimensional cdf for a given value of the modulation, and assuming that the phase is zero.

Parameters:

• phi (float or array) – The (independent) azimuthal angle variable, in [-pi, pi].

• m (float or array) – The modulation, in [0, 1].

static pdf(phi, m)

Evaluate the underlying one-dimensional pdf for a given value of the modulation, and assuming that the phase is identically zero.

Parameters:

• phi (float or array) – The (independent) azimuthal angle variable, in [-pi, pi].

• m (float or array) – The modulation, in [0, 1].

rvs_phi(m, phase)

Generate random variates for any modulation and phase values.

This is essentially calling the underlying xUnivariateAuxGenerator.rvs() method passing the modulation array as an argument and adding the phase array (modulo 2pi) to the output.

Parameters:

• m (array) – An array of modulation values. (The function returns an equal-length array of phi values.)

• phase (float or array) – The phase of the modulation. (This can either be a vector or an array of the same length as modulation.)

class ixpeobssim.irf.modf.xModulationFactor(file_path, extension='fits', k=1)

Class describing the modulation factor.

The effective area is essentially a linear spline, with built-in facilities for evaluation and plotting.

To zero-th order, an xModulationFactor instance is an object capable of evaluating itself in a point or in a grid of points, and capable of plotting itself.

More interestingly, it can generate random phi values, given a vector of event energies and corresponding vectors (or simple floats) representing the polarization degree and angle corresponding to the energies themselves. Internally, any xModulationFactor has an xAzimuthalResponseGenerator object and when the xModulationFactor.rvs_phi() method is called, the polarization degree is multiplied by the modulation factor of the detector, evaluated at the right energy, and converted into a modulation value, after which the underlying xAzimuthalResponseGenerator.rvs_phi() is called.

plot()

Overloaded method for the xpirfview application.

rvs_phi(energy, polarization_degree, polarization_angle)

Return random variates for a given array of values of energy, polarization degree and polarization angle.

Parameters:

• energy (array) – An array of energy values. (The function returns an equal-length array of phi values.)

• polarization_degree (array or float) – The polarization degree, in [0–1]. (This can either be a vector or an array of the same length as energy.)

• polarization_angle (array or float) – The polarization angle, in radians. (This can either be a vector or an array of the same length as energy.)

weighted_average(energy)

Return the weighted average of the mudulation factor given an array of energies.

Parameters:

energy (array) – An array of energy values.

Modulation respose function.

class ixpeobssim.irf.mrf.xModulationResponse(file_path)

Class describing the modulation response, i.e., the product of the effective area times the modulation factor.

plot()

Plot the modulation response.

PSF parametrization.

ixpeobssim.irf.psf.gauss_king(r, W, sigma, N, r_c, eta)

Functional representation of the Gaussian plus King PSF profile described in Fabiani et al., 2014, equation (2):

Parameters:

• r (float or array) – The radial distance from the true source position is arcsec.

• W (float) – Normalization of the Gaussian component.

• sigma (float) – Width of the Gaussian component.

• N (float) – Normalization of the King component.

• r_c (float) – Characteristic radius of the King component.

• eta (float) – Exponent of the King component.

ixpeobssim.irf.psf.gauss_king_eef_at_infinity(W, sigma, N, r_c, eta)

Return the value of the Encircled Energy Fraction (EEF) at infinity, given the parameters of the functional representation, see equation (4) of Fabiani et al., 2014.

Parameters:

• r (float or array) – The radial distance from the true source position is arcsec.

• W (float) – Normalization of the Gaussian component.

• sigma (float) – Width of the Gaussian component.

• N (float) – Normalization of the King component.

• r_c (float) – Characteristic radius of the King component.

• eta (float) – Exponent of the King component.

class ixpeobssim.irf.psf.xPointSpreadFunction(file_path)

Class describing a (simplified, energy independent) PSF.

The effective area is essentially a linear spline, with built-in facilities for evaluation and plotting.

Parameters:

psf_file_path (str) – The path to the .psf FITS file containing the PSF parameters.

Note

The parametrization is taken from Fabiani et al., 2014, table 2.

build_eef()

Build the Encircled Energy Fraction (EEF) as a function of r.

And, while we’re at it, we also calculate and cache the HEW.

delta(size=1)

Return an array of random offset (in ra, dec or L, B) due to the PSF.

Note the output is converted in degrees.

draw_psf_circle(image, x=0.8, y=0.8, label='HPD', hpd_value=False, color='white', lw=1.5)

Add the PSF circle to the image with labels.

Parameters:

• image (xFITSImageBase instance) – The parent xFITSImageBase object

• x (float) – The position of the psf circle in relative canvas coordinates

• y (float) – The position of the psf circle in relative canvas coordinates

plot()

Overloaded plot method.

class ixpeobssim.irf.psf.xPointSpreadFunction2d(file_path)

Two-dimensional version (i.e., non azimuthally symmetric) version of the PSF.

Parameters:

file_path (str) – The file path for the PSF image

build_eef()

Calculate the (azimuthally averaged) encircled energy fraction as a function of r.

delta(size=1)

Return an array of random offset (in ra, dec) due to the PSF.

class ixpeobssim.irf.psf.xPointSpreadFunctionBase(gpd)

Base virtual class for the PSF data structures.

Added in version 28.1.0: This was added in version 27.1.0 as a base class for the old-style PSF (azimuthally symmetric and including the position resolution of the GPD) and the new-style PSF necessary to study the Stokes cross-talk (non azimuthally symmetric and limited to the X-ray optics).

This class holds a single boolean flag inficating whether the GPD position resolution is included in the PSF. Sub-classes should reimplement the delta() hook, returning random arrays of (ra, dec) offsets.

Parameters:

gpd (bool) – Boolean flag indicating whether a given PSF includes the effect of the GPD position resolution.

delta(size=1)

Return an array of random offset (ra, dec) due to the PSF.

smear(ra, dec)

Smear a pair of arrays of coordinates.

smear_single(ra, dec, size=1)

Smear a single pair of coordinates for an arbitrary number of times.

Definition of the response matrix.

class ixpeobssim.irf.rmf.xEnergyDispersion(file_path)

Class representing the energy dispersion.

Parameters:

file_path (str) – The path to the .rmf FITS file containing the energy dispersion tables.

channel_to_energy(channel)

Convenient proxy to the underlying xEnergyDispersionBounds method.

convolve_energy(mc_energy)

Convolve an array of mc_energy with the energy dispersion.

This function is less trivial than the name might suggest, and it is performing several different operations at the same time, namely:

• call xEnergyDispersionMatrix.rvs2() to smear the true energy and get the measured energy (in channel space) both before and after the digitization;

• convert the measured energy before digitization to keV;

• convert the pha to pulse invariant.

energy_to_channel(energy)

Convenient proxy to the underlying xEnergyDispersionBounds method.

pha_analysis(energy)

Perform a pulse-height analysis on a given array of energies.

static pha_to_pi(pha)

Simple placeholder method to convert a PHA value to pulse invariant.

Note

Since we’re not handling the gain corrections, yet, this is simply converting the values from integer to floating point numbers.

plot()

Plot the energy dispersion.

static scale_energy(energy)

Hook to allow to apply distorsions to the energy response of the detector.

In this default implementation this is just returning the energy, untouched, but the method can be overridden to simulate a variety of systematic effects.

class ixpeobssim.irf.rmf.xEnergyDispersionBounds(hdu)

Class encapsulating the bounds for the energy dispersion matrix, as stored in the EBOUNDS extension of a .rmf file.

This is essentially a univariate spline with the channel number on the x axis and the average energy of the corresponding interval on the y axis. At the same time we do store the bounds of the intervals themselves.

channel_to_energy(channel)

Convert a channel number to energy.

energy_to_channel(energy)

Convert a physical energy (in keV) into a channel number.

This is using a simple binary search on the underlying energy bounds.

max()

Return the maximum energy of the last channel.

min()

Return the minimum energy of the first channel.

num_channels()

Return the number of channels.

class ixpeobssim.irf.rmf.xEnergyDispersionMatrix(hdu)

Class encapsulating the energy dispersion matrix, as stored in the MATRIX extension of a .rmf file.

We don’t downsample the matrix to avoid the failure of the fit procedure in XSPEC.

Parameters:

hdu (FITS hdu) – The MATRIX hdu in the .rmf FITS file.

static digitize(value, dtype=<class 'numpy.int16'>)

Digitize an energy value.

This is essentially rounding an array to the nearest integer and cast the result to the specified type.

rvs(aux)

Overloaded method.

This is the original rvs() overloade method that was initially implemented in the energy dispersion code, and it is digitizing the measured energy before returning it, just like the detector would do in real life.

We discovered along the way that it is handy to keep track of the measured (smneared) energy before the digitization (i.e., with all the significant digits), as this allows to apply corrections (e.g., for the GEM charging) without running into rounding problems. This is why we added the rvs2() class method, which is returning both.

rvs2(aux)

Alterntative rvs() implementation where we return the measured energy both before and after the digitization.

Note that the energy before the digitization is still in the pha space, and needs to be explicitely converted to physical units (keV), if that it what one is interested. When that’s the case, the operation is deferred to the xEnergyDispersion class, who is aware of both the energy dispersion matric and the energy bounds.

PSF parametrization.

class ixpeobssim.irf.vign.xVignetting(file_path)

Class describing vignetting.

The vignetting is essentially a bivariate linear spline, with built-in facilities for evaluation and plotting.

Parameters:

file_path (str) – The path to the vign.fits file containing the vignettting table.

plot()

Plot the vignetting.

IRF Generation

Small module providing simulation tools for the alpha particles emitted by the Be window.

ixpeobssim.irfgen.alphas.gas_cell_intersection_points(ray)

Calculate the intersecting points of a given ray with the gas cell.

This is a horrible list of conditions, and should be

ixpeobssim.irfgen.alphas.generate_decays(alpha_energy=5.0, num_alphas=1000000, plot=True, energy_cut=0.1)

Generate the initial alpha particles in the Be window volume.

This is returning a list of xRay object and a numpy array of energies describing the alpha particles that emerge from the window in a direction that is intersecting the fiducial cylinder enclosing the GPD active gas volume.

This first part of the program is vectorized, as only a small fraction of the particles survive the initial cut, and the whole thing would be excruciatingly slow, if stuffed into a plain Python loop.

ixpeobssim.irfgen.alphas.generate_distributions(alpha_energy=5.0, num_alphas=10000000)

ixpeobssim.irfgen.alphas.plot_distributions(alpha_energy=5.0, num_alphas=1000000)

ixpeobssim.irfgen.alphas.plot_event(ray, enery, track_start=None, track_end=None)

Rudimentary single-event display.

ixpeobssim.irfgen.alphas.plot_raindrops(alpha_energy=5.0, num_alphas=2000, max_roi_size=5000)

ixpeobssim.irfgen.alphas.plot_rectangle(width, height, center=(0.0, 0.0), **kwargs)

Facility to plot a generic rectangle.

ixpeobssim.irfgen.alphas.simulate(alpha_energy=5.0, num_alphas=100000, plot=False, energy_cut=0.1)

Run the actual simulation.

ixpeobssim.irfgen.astar.load_alpha_stopping_power_data(identifier, density=None)

Load the stopping power data for a given element or compound.

ixpeobssim.irfgen.astar.load_dme_alpha_stopping_power_data(temperature=20.0, pressure=800.0)

Load the stopping power data for a given element or compound.

class ixpeobssim.irfgen.astar.xAlphaStoppingPowerTable(identifier, density=None)

Basic interface to the output files of the ASTAR database.

bragg_curve(energy=10.0, energy_step=0.01)

Return a spline representing the Bragg curve for the element or compound.

energy_loss(max_energy=10.0, energy_step=0.01)

process_data()

Do all the post-processing of the underlying data.

class ixpeobssim.irfgen.astar.xDMEAlphaStoppingPowerTable(temperature, pressure)

Specialized subclass for the DME compound.

Library for producing auxiliary files for the IRF generation.

ixpeobssim.irfgen.auxiliary.absorption_type_masks(mc_absz, min_dme_absz=0.832, max_dme_absz=10.83)

Return a set of three masks for the window, DME and GEM absorption.

This is based on the Monte Carlo z coordinate of the absorption point, and it is not the most elegant solution, as the correctness of the mask relies on the fact that we use the right bounds for the thickness of the absorption gap used in the simulations. (Since the latter is essentially never changed, this is less horrible than it might seems at a first glance.)

Moving forward, fixing https://bitbucket.org/ixpesw/gpdsw/issues/222 and use the exact information might be the way to go.

ixpeobssim.irfgen.auxiliary.allx_edisp_file_path(abs_label, weight_name, aux_version=3)

Return the path to the allx edisp data file path for a given absorption type and set of weights.

ixpeobssim.irfgen.auxiliary.allx_pha_model_file_path(aux_version=3)

Return the path to the file with the PHA model.

ixpeobssim.irfgen.auxiliary.allx_qeff_file_path(weight_name, aux_version=3)

Return the path to the allx lines data file path for a given absorption type and set of weights.

ixpeobssim.irfgen.auxiliary.allx_rmf_file_path(abs_label, weight_name=None, aux_version=3)

Return the path to the FITS file(s) to be used to compose the response matrix.

ixpeobssim.irfgen.auxiliary.calculate_efficiency(n, N)

Calculate the efficiency, given a number of trials and success, and attach the classical binomial error to it.

ixpeobssim.irfgen.auxiliary.calculate_pi(mc_energy, pha)

Perform the conversion from PHA to PI.

This is necessary because the PHA values in the ixpesim simulation are expressed in the terms of the native detector ADC counts, and the scale factor to the actual energy channels in the response functions is arbitrary. In real life the conversion will be based on the ground calibration with monochromatic sources, and this is meant to be a sensible proxy for that process for operating on Monte Carlo data sets.

Note this operates differently depending on whether we’re working with line data or a continoum spectrum.

Parameters:

• mc_energy (array_like) – The array of (true) energies in keV

• pha (array_like) – The array of measured PHA values.

Returns:

• rec_energy (array_like) – The array of reconstructed energies (in keV) corresponding to the input PHAs.

• pi (array_like) – The array of pulse invariants corresponding to the input PHAs.

• pha_model (xInterpolatedUnivariateSpline object) – The model to operate the conversion between PHAs and PIs.

ixpeobssim.irfgen.auxiliary.calculate_pi_broadband(mc_energy, pha, energy_binning=None)

Calculate the event pulse invariant for a broadband simulation.

This is trying to account for the bulk of the non linearities by fitting the normalized PHA in different energy slices and then fitting the scale factor as a function of the energy with an erf function. (The effect is of the order of 2% from 2 keV to 12 keV, at about 3000 ADC counts per keV nominal.)

ixpeobssim.irfgen.auxiliary.calculate_pi_line(pha, aux_version=3)

Convert the pha column into pulse invariant for a line simulation.

Note that, for line simulations, we use the PHA model derived from the allx data set.

ixpeobssim.irfgen.auxiliary.event_weights(event_data, label='alpha075')

Book-keeping function mapping text labels to weights.

ixpeobssim.irfgen.auxiliary.lines_edisp_file_path(abs_label, weight_name, aux_version=3)

Return the path to the lines edisp data file path for a given absorption type and set of weights.

ixpeobssim.irfgen.auxiliary.lines_qeff_file_path(weight_name, aux_version=3)

Return the path to the lines qeff data file path for a given absorption type and set of weights.

ixpeobssim.irfgen.auxiliary.lines_rmf_file_path(abs_label, weight_name=None, aux_version=3)

Return the path to the FITS file(s) to be used to compose the response matrix.

ixpeobssim.irfgen.auxiliary.load_allx_edisp_data(abs_label, weight_name=None, aux_version=3)

Load the (raw) response matrix data created with mkauxallx.py.

ixpeobssim.irfgen.auxiliary.load_allx_rmf_hist(abs_label, weight_name=None, aux_version=3)

Load the post-processed histogram with the rmf data.

ixpeobssim.irfgen.auxiliary.load_event_data(*file_list, absorption_masks=True, pi=True)

Load the data from the (reconstructed) event files generates with ixpesim.

This is chaining together all the data in the input files—for allx simulations you should pass all the files to a single function call, while for lines you should be open them one by one.

The event data are returned in the form of a dictionary indexed by the corresponding column names in the EVENTS and MONTE_CARLO extensions of the photon lists. Metadata are stored into a separate info dictionary.

ixpeobssim.irfgen.auxiliary.load_lines_edisp_data(abs_label, weight_name=None, aux_version=3)

Load the quantum efficiency data from file.

ixpeobssim.irfgen.auxiliary.load_lines_rmf_hist(abs_label, weight_name=None, aux_version=3)

Load the post-processed histogram with the rmf data.

ixpeobssim.irfgen.auxiliary.load_pha_model(aux_version=3)

Load a PHA to PI model from file.

ixpeobssim.irfgen.auxiliary.load_pressure_scan_table(file_path, col_name='MU2')

Load the pressure scan data from file.

This is reading the binary table stored in the FITS files and arranging them in an order that is suitable for interpolation in pressure and energy.

ixpeobssim.irfgen.auxiliary.load_qeff_table(file_path)

Load the quantum efficiency data from file.

This is reading the binary table stored in the FITS files and returning the columns in the form of a Python dictionary indexed by column name.

ixpeobssim.irfgen.auxiliary.pscan_modf_file_path(weight_name, aux_version=3)

Return the path to the modulation factor data from the analysis of the pressure scan.

ixpeobssim.irfgen.auxiliary.write_pha_model(file_path, model, info)

Write the PHA to PI conversion model table to file.

ixpeobssim.irfgen.auxiliary.write_pressure_scan_table(file_path, data, info)

Write a pressure table to file.

ixpeobssim.irfgen.auxiliary.write_qeff_table(file_path, data, info)

Write a quantum efficiency table to file.

class ixpeobssim.irfgen.auxiliary.xEdispModelGenGauss

Energy dispersion fitting model.

This is the combination of a generalized Gaussian distribution for the main peak, and a hat-shaped bridge for the tail, with the possible addition of a low-energy bump.

static value(x, norm, peak_frac, peak, sigma, beta, bump_frac)

Overloaded value() method.

class ixpeobssim.irfgen.auxiliary.xEdispModelLogNormal

Energy dispersion fitting model.

This is the combination of a log-normal distribution for the main peak, and a hat-shaped bridge for the tail, with the possible addition of a low-energy bump.

static value(x, norm, peak_frac, peak, sigma, shape, endpoint, bump_frac, bump_sigma)

Overloaded value() method.

Base module for all the facilities related to the creation of response functions.

ixpeobssim.irfgen.base.make_arf(irf_name, du_id, pressure, creator, comments, weight_name=None, aux_version=3)

Write the IXPE effective area response function.

ixpeobssim.irfgen.base.make_modf(irf_name, du_id, pressure, creator, comments, weight_name=None, aux_version=3)

Write the IXPE modulation factor response function.

ixpeobssim.irfgen.base.make_mrf(irf_name, du_id, pressure, creator, comments, weight_name=None, aux_version=3)

Write the IXPE effective area response function.

ixpeobssim.irfgen.base.make_psf(irf_name, du_id, params, creator, comments, weight_name=None, aux_version=3)

Write the IXPE PSF response function.

ixpeobssim.irfgen.base.make_rmf(irf_name, du_id, pressure, creator, comments, weight_name=None, aux_version=3)

Write the IXPE edisp response function.

The specifications are described at page ~15 of the following document: ftp://legacy.gsfc.nasa.gov/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.ps

ixpeobssim.irfgen.base.make_vign(irf_name, du_id, pressure, creator, comments, weight_name=None, aux_version=3)

Write the IXPE vignetting response function.

ixpeobssim.irfgen.du.uv_filter_transparency(energy)

Return the transparency of the UV filter evaluated on a given grid of energy points.

ixpeobssim.irfgen.du.uv_filter_transparency_spline()

Return a spline with the transparency of the UV filter as a function of the energy.

ixpeobssim.irfgen.gpd.dme_photoemission_frac_spline(pressure=800.0)

Return a spline representing the fraction of photoelectons extracted from the DME as a function of energy.

This is achieved by calculating the probability of extracting a photoelectron from the Be window and the GEM using the Geant 4 Monte Carlo simulation.

ixpeobssim.irfgen.gpd.load_ixpesim_ancillary_data(file_name='ixpesim_stdlines_20191109.txt')

Load the ancillary data calculated with the full Geant 4 Monte Carlo simulation to inform the generation of the instrument response functions.

The columns represent, in order:

• Energy [keV];

• Trigger efficiency;

• Probability of extraction of a photoelectron from the Be window;

• Probability of extraction of a photoelectron from the GEM Cu upper layer;

• Simulated value of the modulation factor.

ixpeobssim.irfgen.gpd.modulation_factor(energy, pressure, weight_name=None, aux_version=3)

Return the modulation factor calculated on a grid of points.

ixpeobssim.irfgen.gpd.photoabsorption_efficiency(energy, temperature, pressure, thickness=1.0)

Return the photoabsorption efficiency of the GPD absorption gap.

ixpeobssim.irfgen.gpd.quantum_efficiency(energy, temperature, pressure, weight_name=None, aux_version=3)

Return the overall quantum efficiency of the GPD.

Parameters:

• energy (array-like) – The energy array where the quantum efficiency is calculated

• temperature (float) – The GPD filling temperature

• pressure (float) – The GPD gas pressure

• weight_name (str) – The label for the weights to be used.

• aux_version (int) – The version of the auxiliary data products used to model the passive conversions.

ixpeobssim.irfgen.gpd.window_al_transparency(energy, thickness=5.3e-06)

Return the transparency of the 53 nm alumination of the Be window.

ixpeobssim.irfgen.gpd.window_be_transparency(energy, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0}, thickness=0.005)

Return the transparency of the 50 um Be window.

Parameters:

• energy (array-like) – The array of energy values where we calculate the transparency.

• contaminants (dict, optional) – A dictionary containing the list and concentration of contaminants in the Be window.

ixpeobssim.irfgen.gpd.window_transparency(energy, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0})

Return the overall window transparency.

class ixpeobssim.irfgen.gpd.xEdispDataInterface(weight_name=None, aux_version=3)

Basic interface to the post-processed ixpesim output that is necessary to create the response matrix.

combine(temperature=20.0, pressure=800.0, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0})

Combine the window, DME and GEM response functions in the proper proportions (depending on the target temperature and pressure, as well as the window contaminants) to create the actual response function.

class ixpeobssim.irfgen.gpd.xModfDataInterface(weight_name=None, aux_version=3)

Basic interface to the post-processed ixpesim output parametrizing the modulatuion factor as a function of the pressure.

pressure_slice(pressure)

Return a slice at the given pressure, in the form of a piecewise interpolated spline.

class ixpeobssim.irfgen.gpd.xQeffDataInterface(weight_name=None, aux_version=3)

Basic interface to the post-processed ixpesim output that is necessary to create the arf response functions.

This is essentially reading the FITS file with the quantum efficiency tabulated as a function of the energy for the different conversion types (window, DME, and GEM) and caching the necessary arrays in a form that is convenient for later use.

In addition, the class is parametrizing the extraction probabilities for the window and the GEM, so that we can readily scale all the relevant quantities at different internal pressure in a self-consistent fashion without the need to re-run the original simulations.

Warning

When we updated this to allow for weights, it became clear that the efficiency of the active gas volume could not be simply computed from the analystic formula (which was the original approach).

We therefore decided to cache the (energy-dependent) weighting efficiency along with the other quantities, based on the consideration that its dinamic range is comparatively small, and it is therefore more easily amenable to interpolation.

dme_absorption_prob(energy, temperature=20.0, pressure=800.0, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0})

Return the probabilty for a photon of being absorbed in the DME at a given energy (or array of energies).

dme_quantum_efficiency(energy, temperature=20.0, pressure=800.0, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0})

Return the quantum efficiency of the DME.

gem_absorption_prob(energy, temperature=20.0, pressure=800.0, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0})

Return the probabilty for a photon of being absorbed in the GEM at a given energy (or array of energies).

gem_extraction_prob(energy)

Return the photoelectron extraction probability for the GEM at a given energy.

Mind we clip the output of the spline interpolation to the [0, inf] interval, as the spline might accidentally go below zero at low energy due to the noise.

gem_quantum_efficiency(energy, temperature=20.0, pressure=800.0, contaminants={'Ag': 10.0, 'Al': 300.0, 'B': 3.0, 'BeO': 6000.0, 'C': 400.0, 'Ca': 100.0, 'Cd': 2.0, 'Co': 10.0, 'Cr': 100.0, 'Cu': 100.0, 'Fe': 800.0, 'Li': 3.0, 'Mg': 490.0, 'Mn': 100.0, 'Mo': 20.0, 'Ni': 100.0, 'Pb': 20.0, 'Si': 300.0, 'Th': 2.0, 'U': 140.0})

Return the quantum efficiency of the GEM.

Note that, while the GEM extraction probability is pressure-independent, this depends on the pressure due to the effect of the absorptions in the gas, changing the number of X-ray photons reaching the GEM top.