from utils import setup_marx, setup_saotrace
setup_marx()
trace_nest = setup_saotrace()

Point Spread Function (PSF)

The point-spread function (PSF) for Chandra describes how the light from a point source is spread over a larger area on the detector. Several effects contribute to this, e.g. scatter from not absolutely perfectly shaped mirrors, the uncertainty in the pointing, the fact that the detectors are flat, while the focal plane of the mirror is curved (specifically for large off-axis angles) and the pixalization of data on detector read-out.

The following tests compare marx simulations, SAOTrace simulations, and data to look at different aspects of the Chandra PSF.

On-axis PSF on an ACIS-BI chip

The PSF depends on many things, some of which are common to all observations like the shape of the mirror, and some are due to detector effects. For ACIS detectors, the sub-pixel event repositioning (EDSER) can improve the quality of an image by repositioning events based on the event grade. This correction depends on the type pf chip (FI or BI). This test compares the simulation of a point source on a BI ACIS-S chip to an observation. The observed object is TYC 8241 2652 1, a young star, and was observed in 1/8 sub-array mode to reduce pile-up. The pile-up fraction in the data is about 5% in the brightest pixel.

from ciao_contrib.cda.data import download_chandra_obsids
from coords.chandra import cel_to_chandra
from astropy.io import fits

# Download data
out = download_chandra_obsids([15713])

asolfile = '15713/primary/pcadf15713_000N001_asol1.fits.gz'
evt2file = '15713/primary/acisf15713N003_evt2.fits.gz'

from astropy.coordinates import SkyCoord
tyc = SkyCoord.from_name('TYC 8241 2652 1')

coos = cel_to_chandra(fits.getheader(evt2file, 1), tyc.ra.degree, tyc.dec.degree)

x = coos['x'][0]
y = coos['y'][0]
r = 4
region = f'circle({x},{y},{r})'
psffrac = 1.0

bg_x = 4044
bg_y = 4123

from utils import write_spectrum_from_fluxcorrection

write_spectrum_from_fluxcorrection(asolfile, evt2file, x, y, region, psffrac=psffrac)

WARNING: UnitsWarning: 'channel' did not parse as fits unit: At col 0, Unit 'channel' not supported by the FITS standard.  If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

from utils import marxpars_from_asol

pars = marxpars_from_asol(asolfile, evt2file)
pars['SpectrumType'] = 'FILE'
pars['SpectrumFile'] = 'input_spec_marx.tbl'
pars['OutputDir'] = 'marx_only'
pars['SourceFlux'] = -1

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)
marxcall

'marx RA_Nom=182.25561765981 Dec_Nom=-51.34945070544 Roll_Nom=33.971360241142 GratingType=NONE ExposureTime=10190.8068100214 DitherModel=FILE DitherFile=15713/primary/pcadf15713_000N001_asol1.fits TStart=508966864.84768 ACIS_Exposure_Time=0.7 SourceRA=182.259167 SourceDEC=-51.344667 DetectorType=ACIS-S DetOffsetX=0.0014449422646707344 DetOffsetZ=-0.007542945902798692 SpectrumType=FILE SpectrumFile=input_spec_marx.tbl OutputDir=marx_only SourceFlux=-1'

'''Currently, MARX can only be run from the command line. So we build
a command line call from the parameters dictionary and run it via
subprocess.call().'''
import subprocess
out = subprocess.run(marxcall, shell=True, capture_output=True)

asolh = fits.getheader(asolfile, 1)
evth = fits.getheader(evt2file, 1)

def write_saotrace_lua(asolfile, asolh, evth, spectrumfile='input_spec_saotrace.rdb'):
    """Write an SAOTrace Lua source parameter file."""
    lua_text = f'''
    ra_pnt = {asolh['RA_NOM']}
    dec_pnt = {asolh['DEC_NOM']}
    roll_pnt = {asolh['ROLL_NOM']}

    dither_asol_chandra{{ file = "{asolfile}",
                        ra = ra_pnt, dec = dec_pnt, roll = roll_pnt }}

    point{{ position = {{ ra = {evth['RA_TARG']},
            dec = {evth['DEC_TARG']},
            ra_aimpt = ra_pnt,
            dec_aimpt = dec_pnt,
        }},
        spectrum = {{ {{ file = "{spectrumfile}",
                        units = "photons/s/cm2",
                        scale = 1,
                        format = "rdb",
                        emin = "ENERG_LO",
                        emax = "ENERG_HI",
                        flux = "flux"}} }}
        }}
    '''
    with open('saotrace_source.lua', 'w') as f:
        f.write(lua_text)

write_saotrace_lua(asolfile, asolh, evth)

'''
CXO time numbers are large and round-off error can appear
which causes SAOTrace to fail.
Therefore, shorten all times by about 0.01 sec to make sure.
'''
limit = asolh['TSTOP'] - asolh['TSTART']
limit = limit - 0.02

out = subprocess.run(trace_nest +
    f'tag=saotrace srcpars=saotrace_source.lua tstart={asolh["TSTART"] + 0.01} limit={limit} limit_type=sec', shell=True)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

pars = marxpars_from_asol(asolfile, evt2file)
pars['OutputDir'] = 'marx_saotrace'
pars['SourceType'] = 'SAOSAC'
pars['SAOSACFile'] = 'saotrace.fits'
pars['SourceFlux'] = -1

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)

out = subprocess.run(marxcall, shell=True, capture_output=True)

out = subprocess.run('marx2fits --pixadj=EDSER marx_only marx_only.fits', shell=True, 
                     capture_output=True)
out2 = subprocess.run('marx2fits --pixadj=EDSER marx_saotrace marx_saotrace.fits', shell=True, 
                     capture_output=True)
out3 = subprocess.run('marx2fits --pixadj=RANDOMIZE marx_only marx_only_rand.fits', shell=True, 
                     capture_output=True)
out4 = subprocess.run('marx2fits --pixadj=RANDOMIZE marx_saotrace marx_saotrace_rand.fits', shell=True, 
                     capture_output=True)

'''Reprocess the observation with RANDOMIZE'''
import ciao_contrib.runtool as rt
rt.acis_process_events(infile=evt2file, outfile='obs_rand.fits', acaofffile=asolfile,
                        doevtgrade=False, calculate_pi=False, apply_cti=False, pix_adj='RANDOMIZE', clobber=True)

# acis_process_events (CIAO 4.18.0): WARNING: The tgain adjustment will not be (re)computed because the doevtgrade is set to no
# acis_process_events (CIAO 4.18.0): WARNING: The values of PHA, ENERGY, PI, FLTGRADE and GRADE may be inaccurate because doevtgrade=no.
# acis_process_events (CIAO 4.18.0): WARNING: The values of ENERGY and PI may be inaccurate because calculate_pi=no.

import numpy as np
from astropy.stats import sigma_clipped_stats
from astropy.table import Table, Column
import matplotlib.pyplot as plt

# on import, this registers the plotting scale "power"
import plot_utils

def radial_distribution(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    '''Calculate the radial distance from the mean position for events.

    Parameters
    ----------
    x, y : np.array
        x and y coordinates of events

    Returns
    -------
    r : np.array
        radial distance from mean position of the events.
    '''
    centx = x - sigma_clipped_stats(x, sigma=1.5, maxiters=20)[0]
    centy = y - sigma_clipped_stats(y, sigma=1.5, maxiters=20)[0]

    xy = np.vstack([centx, centy])
    return np.linalg.norm(xy, axis=0)


def plot_ecf(ax:plt.Axes, 
             files: list[str], 
             filterargs: dict, 
             bgfilterargs: dict) -> tuple[float, float]:
    '''
    Parameters
    ----------
    ax : matplotlib axes
        Axes where the plot will be placed
    file : list
        List of three strings with filename for observations,
        Marx simulation, and SAOTRace + marx simulation
    filterargs : dict
        See ``filter_events`` for details
    bgfilterargs : dict
        Select a background region (same size, no automatic scaling)

    Returns
    -------
    err_marx : float
        Flux error if the MARX simulated PSF is assumed to be right.
        This is calculated as follows: We find the radius that encircles 90% of
        all counts in the MARX simulation. Then, we extract the observed counts
        using that radius. If the simulation is correct, that radius should
        contain 90% of the observed counts, too. If, e.g. the simulated radius
        is too small, we may extract only 80 % of the counts. The ratio between
        the two would be 8/9=0.88, meaning that all fluxes extracted using
        this radius are 12 % too small.
    err_sao: float
        Same for the simulation that used SAOTrace + MARX.
    '''
    # In this energy range we have the most counts and we limited
    # the simulations to the same range.
    obs = filter_events(Table.read(files[0], hdu=1), **filterargs)
    bkg = filter_events(Table.read(files[0], hdu=1), **bgfilterargs)
    r_obs = radial_distribution(obs['x'], obs['y'])
    r_obs.sort()
    ax.plot(r_obs, np.linspace(0, 1 + len(bkg) / len(r_obs), len(r_obs)),
            label='Observation')

    simmarx = filter_events(Table.read(files[1], hdu=1), **filterargs)
    r_marx = radial_distribution(simmarx['X'], simmarx['Y'])
    r_marx.sort()
    ax.plot(r_marx, np.linspace(0, 1, len(r_marx)), label='MARX')

    simsao = filter_events(Table.read(files[2], hdu=1), **filterargs)
    r_sao = radial_distribution(simsao['X'], simsao['Y'])
    r_sao.sort()
    ax.plot(r_sao, np.linspace(0, 1, len(r_sao)), label='SAOTrace + MARX')

    ax.set_xscale('power', power=0.5)
    ax.set_xlabel('radius [pixel]')
    ax.set_ylabel('enclosed count fraction')
    ax.legend(loc='lower right')
    ax.set_xticks([.1, .4, .7, 1, 2, 3, 4, 5])
    ax.set_xlim([0.1, 5])
    ax.set_ylim([0, 1.])

    psf90_marx = np.percentile(r_marx, 90)
    ecf_at_that_rad = np.argmin(np.abs(r_obs - psf90_marx)) / len(r_obs)
    err_marx = ecf_at_that_rad / 0.9 - 1

    psf90_sao = np.percentile(r_sao, 90)
    ecf_at_that_rad = np.argmin(np.abs(r_obs - psf90_sao)) / len(r_obs)
    err_sao = ecf_at_that_rad / 0.9 - 1

    return err_marx, err_sao


def colname_case(table: Table, name: str) -> Column:
    '''Extract a column from a Table case insensitive to the column name.

    If several columns in the table match ``name`` in a case-insensitive way,
    this just returns the first match.

    Parameters
    ----------
    table : `astropy.table.Table`
        Table with columns
    name : string
        column name, case insensitive.

    Returns
    -------
    col : `astropy.table.Column`
        Table column
    '''
    for n in table.colnames:
        if name.lower() == n.lower():
            return table[n]
    return table[name]


def filter_events(evt: Table, 
                  circle: tuple | None = None, 
                  energy: tuple | None = None, 
                  time: tuple | None = None) -> Table:
    '''A simple event filter function.

    The parameters offer several filters, it set to ``None`` a filter is not
    applied.
    This simple function does not duplicate the entire CIAO syntax, only
    a few simple filters useful for plotting in pure-python functions.

    Parameters
    ----------
    evt : `astropy.table.Table`
        Event table
    circle : tuple or ``None``
        Tuple of the form (x, y, r) where (x, y) is the center of
        a circle and r the radius. All number are in detector coordinates.
    energy : tuple or ``None``
        Tuple of the form (lower, upper) with all energies in eV.
    time : tuple or ``None``
        Tuple of the form (star, end) with all times in seconds.

    Returns
    -------
    evt_filt : ~astropy.table.Table`
        Filtered event list
    '''
    if circle is not None:
        d = np.sqrt((colname_case(evt, 'x') - circle[0])**2 + (colname_case(evt, 'y') - circle[1])**2)
        indcirc = d < circle[2]
    else:
        indcirc = np.ones(len(evt), dtype=bool)

    if energy is not None:
        en = colname_case(evt, 'energy')
        inden = (en > energy[0]) & (en < energy[1])
    else:
        inden = np.ones_like(indcirc)

    if time is not None:
        t = colname_case(evt, 'time')
        indt = (t > time[0]) & (t < time[1])
    else:
        indt = np.ones_like(indcirc)

    return evt[indcirc & inden & indt]

filterargs = {'energy': [300, 3000], 'circle': (x, y, 5 * r)}
bgfilterargs = {'energy': [300, 3000], 'circle': (bg_x, bg_y, 5 * r)}

fig, axes = plt.subplots(figsize=(6, 3), ncols=2)

for ax, files, title in zip(axes,
                            [[evt2file, 'marx_only.fits', 'marx_saotrace.fits'],
                                ['obs_rand.fits', 'marx_only_rand.fits', 'marx_saotrace_rand.fits']],
                            ['EDSER', 'RANDOMIZE']):
    err_marx, err_sao = plot_ecf(ax, files, filterargs, bgfilterargs)

    ax.set_title(title)
    ax.grid()

WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

Enclosed count fraction for observation and simulations. The simulations are run with a count number similar to the (fairly short) observation and thus there is some wiggling due to Poisson noise. On the left, the results using the default EDSER algorithms are shown, on the right, the sub-pixel positions are randomized within the pixel independent of the flight grade, which smooths out the inner PSF and hides the problem that marx does not make good predictions for the flight grades.

For BI chips (here ACIS-S3) the marx simulated PSF is too narrow. This effect is most pronounced at small radii around 1 pixel. At larger radii the marx simulation comes closer to the observed distribution. Running the simulation with a combination of SAOTrace and marx gives PSF distributions that are closer to the observed numbers. However, in the range around 1 pix, the simulations are still too wide. Depending on the way sub-pixel information is handled, there is a notable difference in the size of the effect. Using energy-dependent sub-pixel event repositioning (EDSER) requires not only a good mirror model, but also a realistic treatment of the flight grades assigned on the detector. This is where the difference between simulation and observations is largest.

The figures show that there are at least two factors contributing to the difference: The mirror model and problems in the simulation of the EDSER algorithm.

The put those numbers into perspective: If I used the marx simulation to determine the size of the extraction region that encloses 80% of all source counts, I would find a radius close to 0.9 pixel. However, in reality, such a region will only contain about 65% of all flux, causing me to underestimate the total X-ray flux in this source by 15%. (The exact number depends on the spectrum of the source in question.)

On-axis PSF on an ACIS-FI chip

Same as above for a FI chip. The target of this observations is tau Canis Majores.

# Download data
out = download_chandra_obsids([4469])

asolfile = '4469/primary/pcadf04469_000N001_asol1.fits.gz'
evt2file = '4469/primary/acisf04469N004_evt2.fits.gz'

sky_coos = SkyCoord.from_name('tau Canis Majoris')
coos = cel_to_chandra(fits.getheader(evt2file, 1), sky_coos.ra.degree, sky_coos.dec.degree)

x = coos['x'][0]
y = coos['y'][0]
r = 4
region = f'circle({x},{y},{r})'
psffrac = 1.0

bg_x = 4182
bg_y = 4127

write_spectrum_from_fluxcorrection(asolfile, evt2file, x, y, region, psffrac=psffrac)

WARNING: UnitsWarning: 'channel' did not parse as fits unit: At col 0, Unit 'channel' not supported by the FITS standard.  If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

pars = marxpars_from_asol(asolfile, evt2file)
pars['SpectrumType'] = 'FILE'
pars['SpectrumFile'] = 'input_spec_marx.tbl'
pars['OutputDir'] = 'marx_only'
pars['SourceFlux'] = -1

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)

out = subprocess.run(marxcall, shell=True, capture_output=True)

asolh = fits.getheader(asolfile, 1)
evth = fits.getheader(evt2file, 1)
write_saotrace_lua(asolfile, asolh, evth)

'''
CXO time numbers are large and round-off error can appear
which make SAOTrace fail.
Therefore, shorten all times by about 0.01 sec to make sure.
'''
limit = asolh['TSTOP'] - asolh['TSTART']
limit = limit - 0.02

out = subprocess.run(trace_nest +
    f'tag=saotrace srcpars=saotrace_source.lua tstart={asolh["TSTART"] + 0.01} limit={limit} limit_type=sec', shell=True)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

pars['OutputDir'] = 'marx_saotrace'
pars['SourceType'] = 'SAOSAC'
pars['SAOSACFile'] = 'saotrace.fits'

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)

out = subprocess.run(marxcall, shell=True, capture_output=True)

out = subprocess.run('marx2fits --pixadj=EDSER marx_only marx_only.fits', shell=True, 
                     capture_output=True)
out2 = subprocess.run('marx2fits --pixadj=EDSER marx_saotrace marx_saotrace.fits', shell=True, 
                     capture_output=True)
out3 = subprocess.run('marx2fits --pixadj=RANDOMIZE marx_only marx_only_rand.fits', shell=True, 
                     capture_output=True)
out4 = subprocess.run('marx2fits --pixadj=RANDOMIZE marx_saotrace marx_saotrace_rand.fits', shell=True, 
                     capture_output=True)

rt.acis_process_events(infile=evt2file, outfile='obs_rand.fits', acaofffile=asolfile,
                        doevtgrade=False, calculate_pi=False, apply_cti=False, pix_adj='RANDOMIZE', clobber=True)

# acis_process_events (CIAO 4.18.0): WARNING: The tgain adjustment will not be (re)computed because the doevtgrade is set to no
# acis_process_events (CIAO 4.18.0): WARNING: The values of PHA, ENERGY, PI, FLTGRADE and GRADE may be inaccurate because doevtgrade=no.
# acis_process_events (CIAO 4.18.0): WARNING: The values of ENERGY and PI may be inaccurate because calculate_pi=no.

filterargs = {'energy': [300, 3000], 'circle': (x, y, 5 * r)}
bgfilterargs = {'energy': [300, 3000], 'circle': (bg_x, bg_y, 5 * r)}


fig, axes = plt.subplots(figsize=(8, 4), ncols=2)

for ax, files, title in zip(axes,
                            [[evt2file, 'marx_only.fits', 'marx_saotrace.fits'],
                                ['obs_rand.fits', 'marx_only_rand.fits', 'marx_saotrace_rand.fits']],
                            ['EDSER', 'RANDOMIZE']):
    err_marx, err_sao = plot_ecf(ax, files, filterargs, bgfilterargs)

    ax.set_title(title)
    ax.grid()

WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'Chan' did not parse as fits unit: At col 0, Unit 'Chan' not supported by the FITS standard. Did you mean chan? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

Enclosed count fraction for observation and simulations in the same format as above, but for an ACIS-I observation.

As for the BI chip, marx simulations indicate a PSF that is narrower than the observed distribution. Using SAOTrace as a mirror model, gives better results but still differs significantly from the observed distribution.

On-axis PSF for an HRC-I observation

Same as above for an HRC-I observations of AR Lac.

Since the HRC has no intrinsic energy resolution, the source spectrum can not be taken from the observed data. Instead, we fitted a model to a grating spectrum of the same source and saved that model spectrum to a file. AR Lac is known to be time variable and to show stellar flares, so this spectrum is only an approximation, but it is the best we can do here.

# Download data
out = download_chandra_obsids([13182])

asolfile = '13182/primary/pcadf13182_000N001_asol1.fits.gz'
evt2file = '13182/primary/hrcf13182N004_evt2.fits.gz'

sky_coos = SkyCoord.from_name('AR Lac')
coos = cel_to_chandra(fits.getheader(evt2file, 1), sky_coos.ra.degree, sky_coos.dec.degree)

x, y = coos['x'][0], coos['y'][0]
r = 15
region = f'circle({x},{y},{r})'
psffrac = 1.0

bg_x, bg_y = 16043, 16700

pars = marxpars_from_asol(asolfile, evt2file)
pars['SpectrumType'] = 'FILE'
pars['SpectrumFile'] = 'data/ARLac_input_spec_marx.tbl'
pars['OutputDir'] = 'marx_only'
pars['SourceFlux'] = -1

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)
out = subprocess.run(marxcall, shell=True, capture_output=True)

marxcall

'marx RA_Nom=332.16369387119 Dec_Nom=45.739883051807 Roll_Nom=301.39740584162 GratingType=NONE ExposureTime=18161.46342998743 DitherModel=FILE DitherFile=13182/primary/pcadf13182_000N001_asol1.fits TStart=408913236.56726 SourceRA=332.17 SourceDEC=45.742306 DetectorType=HRC-I DetOffsetX=0.0014262321171452097 DetOffsetZ=-0.0025143152978017724 SpectrumType=FILE SpectrumFile=data/ARLac_input_spec_marx.tbl OutputDir=marx_only SourceFlux=-1'

asolh = fits.getheader(asolfile, 1)
evth = fits.getheader(evt2file, 1)
write_saotrace_lua(asolfile, asolh, evth, spectrumfile='data/ARLac_input_spec_saotrace.rdb')

'''
CXO time numbers are large and round-off error can appear
which make SAOTrace fail.
Therefore, shorten all times by about 0.01 sec to make sure.
'''
limit = asolh['TSTOP'] - asolh['TSTART']
limit = limit - 0.02

out = subprocess.run(trace_nest +
    f'tag=saotrace srcpars=saotrace_source.lua tstart={asolh["TSTART"] + 0.01} limit={limit} limit_type=sec', shell=True)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

pars['OutputDir'] = 'marx_saotrace'
pars['SourceType'] = 'SAOSAC'
pars['SAOSACFile'] = 'saotrace.fits'

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)
out = subprocess.run(marxcall, shell=True, capture_output=True)

# HRC does not record grade information, so EDSER is not applicable
out3 = subprocess.run('marx2fits --pixadj=RANDOMIZE marx_only marx_only.fits', shell=True, 
                     capture_output=True)
out4 = subprocess.run('marx2fits --pixadj=RANDOMIZE marx_saotrace marx_saotrace.fits', shell=True, 
                     capture_output=True)

filterargs = {'circle': (x, y, 5 * r)}
bgfilterargs =  {'circle': (bg_x, bg_y, 5 * r)}
files = [evt2file, 'marx_only.fits', 'marx_saotrace.fits']

fig, ax = plt.subplots()

err_marx, err_sao = plot_ecf(ax, files, filterargs, bgfilterargs)
ax.set_xlim([0, 20.])
ax.grid()

Enclosed count fraction for observation and simulations.

The PSF predicted by SAOTrac + marx tracks the general shape pf the observed PSF well, though differences are apparent in the core and in the wings. Like in ACIS, the PSF simulated with marx alone is narrower than the observed PSF.

Off-axis PSF

For an off-axis source the differences between detector pixel size and different event repositioning algorithms is not important any longer. Instead, the size of the PSF is dominated by optics errors because the detector plane deviates from the curved focal plane and because a Wolter type I optic has fundamental limitations for off-axis sources.

The PSF is much larger no longer circular. As this test shows, the shadows of the support struts become visible.

out = download_chandra_obsids([1068])

asolfile = '1068/primary/pcadf01068_000N001_asol1.fits.gz'
evt2file = '1068/primary/acisf01068N004_evt2.fits.gz'

x, y = 5246, 6890
r = 200
region = f'circle({x},{y},{r})'
psffrac = 1.0

write_spectrum_from_fluxcorrection(asolfile, evt2file, x, y, region, psffrac=psffrac)

WARNING: UnitsWarning: 'channel' did not parse as fits unit: At col 0, Unit 'channel' not supported by the FITS standard.  If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

pars = marxpars_from_asol(asolfile, evt2file)
pars['SpectrumType'] = 'FILE'
pars['SpectrumFile'] = 'input_spec_marx.tbl'
pars['OutputDir'] = 'marx_only'
pars['SourceFlux'] = -1

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)
out = subprocess.run(marxcall, shell=True, capture_output=True)

asolh = fits.getheader(asolfile, 1)
evth = fits.getheader(evt2file, 1)
write_saotrace_lua(asolfile, asolh, evth)

'''
CXO time numbers are large and round-off error can appear
which make SAOTrace fail.
Therefore, shorten all times by about 0.01 sec to make sure.
'''
limit = asolh['TSTOP'] - asolh['TSTART']
limit = limit - 0.02

out = subprocess.run(trace_nest +
    f'tag=saotrace srcpars=saotrace_source.lua tstart={asolh["TSTART"] + 0.01} limit={limit} limit_type=sec', shell=True)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

pars['OutputDir'] = 'marx_saotrace'
pars['SourceType'] = 'SAOSAC'
pars['SAOSACFile'] = 'saotrace.fits'

marxcall = ['marx'] + ['{0}={1}'.format(k, v) for k, v in pars.items()]
marxcall = ' '.join(marxcall)

out = subprocess.run(marxcall, shell=True, capture_output=True)

out = subprocess.run('marx2fits --pixadj=EDSER marx_only marx_only.fits', shell=True, 
                     capture_output=True)
out2 = subprocess.run('marx2fits --pixadj=EDSER marx_saotrace marx_saotrace.fits', shell=True, 
                     capture_output=True)

from astropy.visualization import LogStretch
from astropy.visualization.mpl_normalize import ImageNormalize
import pycrates


fig, axes = plt.subplots(1, 3, figsize=(10,6), subplot_kw={'aspect':'equal'})
for ax, src, title in zip(axes.flatten(), 
                   [evt2file, "marx_only.fits", "marx_saotrace.fits"],
                   ['Observation', 'MARX only', 'SAOTrace + MARX']):
    evt = pycrates.read_file(src)
    out = ax.hist2d(evt.get_column('X').values, evt.get_column('Y').values, 
                    bins=[400, 300], range=[[x-200,x+200], [y-150,y+150]],
                    norm=ImageNormalize(stretch=LogStretch(), vmin=1, vmax=20))
    out = ax.set_title(title)
    ax.set_axis_off()

The structure of the PSFs is very similar, emphasizing how good both mirror models are. On closer inspection, there is a small shadow just above and to the right of the point where the support strut shadows meet. This feature is a little smaller in marx than in SAOTrace or the real data due to the simplification that the marx mirror model makes.

On-axis PSF at different energies

The gold standard to test the fidelity of marx and SAOTrace obviously is to compare simulations to observations. However, it is also instructive to look at a few idealized cases with no observational counterpart so we can simulate high fidelity PSFs with a large number of counts without worrying about background or pile-up. The Chandra Proposers Observatory Guide. contains a long section on the PSF and encircled energy based on SAOTrace simulations. Some of those simulations are repeated here to compare them to pure marx simulations.

energy = [0.25, 0.5, 1, 2, 3, 4, 6, 8]

# Default marx settings for all simulations
default_marx ='GratingType=NONE DetectorType=HRC-I ' + \
              'ExposureTime=10000 ' + \
              'SourceRA=0. SourceDEC=0. ' + \
              'RA_Nom=0. Dec_Nom=0. Roll_Nom=0. '

import subprocess
for e in energy:
    out = subprocess.run(f'marx MinEnergy={e} MaxEnergy={e} DitherModel=INTERNAL OutputDir="marx_{e}" ' +
                    # Higher energies have larger PSF and thus need more photons to plot smooth
                    # contours. This scaling is a compromise between runtime and looking good.
                         f'SourceFlux={0.01*np.sqrt(e)} ' + 
                     default_marx, shell=True, capture_output=True)
    out2 = subprocess.run(f'marx2fits --pixadj=NONE marx_{e} marx{e}.fits', 
                          shell=True, capture_output=True)

'''Generate asol file

Since all simulations use the same pointing and exposure time, it is
enough to run `marxasp` once.
'''
out = subprocess.run(f'marxasp MarxDir="marx_{energy[0]}" ',
                     shell=True, capture_output=True)

asolh = fits.getheader('sim_asol.fits', 1)

limit = asolh['TSTOP'] - asolh['TSTART']
limit = limit - 0.02


for e in energy:
    lua_text = f'''
point{{ position = {{ ra = 0.,
        dec = 0.,
        ra_aimpt = 0.,
        dec_aimpt = 0.,
        }},
        spectrum = {{ {{ {e}, {0.01 * np.sqrt(e)} }} }} 
    }}
roll(0)
dither_asol_chandra{{ file = "sim_asol.fits", ra = 0., dec = 0., roll = 0. }}
'''
    with open('saotrace_source.lua', 'w') as f:
        f.write(lua_text)
    subprocess.run(trace_nest +
                   f'tag=sao{e} srcpars=saotrace_source.lua ' +\
                   f'tstart={asolh["TSTART"] + 0.01} limit={limit} limit_type=sec', shell=True)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

for e in energy:
    out = subprocess.run(f'marx SourceType=SAOSAC SAOSACFile="sao{e}.fits" DitherModel=FILE DitherFile="sim_asol.fits" OutputDir="saomarx_{e}" ' +
                     default_marx, shell=True, capture_output=True)
    out2 = subprocess.run(f'marx2fits --pixadj=NONE saomarx_{e} saomarx_{e}.fits', 
                          shell=True, capture_output=True)

from mpl_toolkits.axes_grid1 import AxesGrid
from matplotlib.ticker import StrMethodFormatter, FuncFormatter

fig = plt.figure(figsize=(10, 4))

grid = AxesGrid(fig, 111,  # similar to subplot(142)
                nrows_ncols=(2, len(energy)),
                axes_pad=0.0,
                share_all=True,
                label_mode="L",
                cbar_location="top",
                cbar_mode="single"
                )

for i, e in enumerate(energy):

    for prog in ['marx', 'saomarx_']:
        f = f'{prog}{e}.fits'
        offset = 0 if prog == 'marx' else len(energy)
        
        with fits.open(f) as hdus:
            header = hdus[1].header
        tab = Table.read(f, hdu=1)
        # The old question: Is an index the center of a pixel or the corner?
        # Differs by 0.5...
        d_ra = (tab['X'] - header['TCRPX9'] - 0.5) * header['TCDLT9'] * 3600.
        # And do I count from the left or right (RA is reversed on the sky)?
        d_dec = (tab['Y'] - header['TCRPX10'] + 0.5) * header['TCDLT10'] * 3600.

        im = grid[i + offset].hist2d(d_ra, d_dec,
                                        range=[[-1, 1], [-1, 1]],
                                        bins=[20, 20],
                                        cmap=plt.get_cmap('OrRd'))
        levels = np.max(im[0]) * np.array([0.3, 0.6, 0.9])
        grid[i + offset].contour(0.5 * (im[1][1:] + im[1][:-1]),
                                    0.5 * (im[2][1:] + im[2][:-1]),
                                    # for some reason, the data is
                                    # ordered for other orientation
                                    im[0].T, levels=levels, colors='k') #,
                                    #origin=im[3].origin)
        grid[i + offset].grid()
        if prog == 'marx':
            grid[i].text(-.5, .5, '{0} keV'.format(e))
        else:
            grid[i].xaxis.set_major_formatter(StrMethodFormatter('{x:2.1g}'))

grid[0].yaxis.set_major_formatter(StrMethodFormatter('{x:2.1g}'))
grid[len(energy)].yaxis.set_major_formatter(StrMethodFormatter('{x:2.1g}'))
cbar = grid.cbar_axes[0].colorbar(im[3],
             format=FuncFormatter(lambda x, pos: ''))

# This affects all axes as share_all = True.
grid.axes_llc.set_xticks([-.5, 0, .5])
grid.axes_llc.set_yticks([-.5, 0, .5])

grid[0].set_ylabel('Marx')
grid[offset].set_ylabel('Marx +\nSAOTrace')
grid[int(1.5 * offset)].set_xlabel('arcsec (measured from nominal source position)')

fig.subplots_adjust(top=1, bottom=0, left=0.1, right=0.99)

PSF images for different discrete energies for pure marx and SAOTrace + marx simulations. The color scale is linear, but the absolute scaling is different for different images, because the effective area and thus the number of detected photons is lower at higher energies. Contour lines mark flux levels at 30%, 60%, and 90% of the peak flux level. At higher energies, the PSF becomes asymmetric, because mirror pair 6, which is most important at high energies, is slightly tilted with respect to the nominal aimpoint. At the same time, the scatter, and thus the size of the PSF, increase at higher energies.

Because the simulations are done with Monte-Carlo methods and the number of photons simulated is finite (but large compared to typical observations), there is some noise in the images.

fig = plt.figure(figsize=(7, 4))
axecf = fig.add_subplot(111)
color = plt.cm.jet_r(np.linspace(0, 1, len(energy)))
for i, e in enumerate(energy):

    for prog in ['marx', 'saomarx_']:
        f = f'{prog}{e}.fits'
        with fits.open(f) as hdus:
            header = hdus[1].header

        tab = Table.read(f, hdu=1)
        # The old question: Is an index the center of a pixel of the corner?
        # Differs by 0.5...
        d_ra = (tab['X'] - header['TCRPX9'] - 0.5) * header['TCDLT9'] * 3600.
        # Count from the left or right (RA is reversed on the sky)?
        d_dec = (tab['Y'] - header['TCRPX10'] + 0.5) * header['TCDLT10'] * 3600.

        # The simulation is set up to make this simple: RA=DEC=0
        # so cos(DEC) = 1 and we can approximate with Euklidian distance
        r = np.linalg.norm(np.vstack([d_ra, d_dec]), axis=0)
        val, edges = np.histogram(r, range=[0, 5], bins=50)
        bin_mid_marx = 0.5 * (edges[:-1] + edges[1:])
        ecf_marx = 1.0 * val.cumsum() / val.sum()
        if prog == 'marx':
            line, = axecf.plot(bin_mid_marx, ecf_marx, color=color[i],
                                label='{0} keV'.format(e))
        else:
            axecf.plot(bin_mid_marx, ecf_marx, color=line.get_color(), ls=':')
axecf.set_xscale('power', power=0.5)
axecf.legend(loc='lower right')
axecf.set_ylabel('encircled count fraction')
axecf.set_xlabel('radius [arcsec]')
axecf.set_xticks([0, .1, .2, .4, .6, .8, 1, 2, 3, 4, 5])
axecf.grid()

A different way to present the width of the PSF is the encircled count fraction. For each photon, the radial distance is calculated from the nominal source position. Because the center is offset for hard photons, the PSF appears wider at those energies. solid line: marx only, dotted lines: SAOTrace + marx.

marx and SAOTrace simulations predict a very similar PSF shape, but for most energies the SAOTrace model predicts a slightly broader PSF.

Simulated off-axis PSF

Compare marx and SAOTrace + marx simulations for different off-axis angles. For simplicity, this is done here for a single energy. We pick 4 keV because it is roughly the center of the Chandra band.

off_axis_angle = [0, .5, 1, 2, 3, 5, 7, 10, 15]  # in arcmin

# Default marx settings for all simulations\
default_marx ='GratingType=NONE DetectorType=HRC-I ' + \
              'ExposureTime=10000 SourceDEC=0. ' + \
              'RA_Nom=0. Dec_Nom=0. Roll_Nom=0. ' + \
              'MinEnergy=4 MaxEnergy=4 SourceFlux=0.01 '

import subprocess
for a in off_axis_angle:
    out = subprocess.run(f'marx DitherModel=INTERNAL OutputDir="marx_{a}" ' +
                         f'SourceRA={a / 60.} ' +
                     default_marx, shell=True, capture_output=True)
    out2 = subprocess.run(f'marx2fits --pixadj=NONE marx_{a} marx{a}.fits', 
                          shell=True, capture_output=True)

# We can reuse the asol file from above since we set up the same pointing, so the
# variable asolh, limit are still valid.
for a in off_axis_angle:
    lua_text = f'''
point{{ position = {{ ra = {a / 60.},
        dec = 0.,
        ra_aimpt = 0.,
        dec_aimpt = 0.,
        }},
        spectrum = {{ {{ 4, 0.01 }} }} 
    }}
roll(0)
dither_asol_chandra{{ file = "sim_asol.fits", ra = 0., dec = 0., roll = 0. }}
'''
    with open('saotrace_source.lua', 'w') as f:
        f.write(lua_text)
    subprocess.run(trace_nest +
                   f'tag=sao{a} srcpars=saotrace_source.lua ' +\
                   f'tstart={asolh["TSTART"] + 0.01} limit={limit} limit_type=sec', shell=True)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

Error in parameter file: line 1: field Name
# 331: multilayer_reflect: /opt/saotrace/2.0.6.1/bin/../xbin/multilayer_reflect: error opening parameter file
# 331: multilayer_reflect: Stack:
# 334: /opt/saotrace/2.0.6.1/bin/../xbin/scatter: pipe error
# 334: /opt/saotrace/2.0.6.1/bin/../xbin/scatter: Stack:
# 337: aperture: unable to attach input stream `stdin' to bpipe: pipe error
# 337: aperture: stack: init_bpipe<-process_all_assemblies<-<TOP>
# 340: projray: couldn't open input file `stdin': pipe error
# 340: projray: stack: process<-<TOP>
# 343: quef: stdin: pipe error
# 343: quef: stack: init_bpipe<-quantum_efficiency<-<TOP>
# 345: tot_wt: couldn't open input file `stdin': pipe error
# 345: tot_wt: stack: tot_wt<-<TOP>
# 349: bpmanip: binary pipe input init failure: pipe error
# 349: bpmanip: Stack:
# 349: bpmanip:  [ 1] BinaryPipe::BinaryPipe
# 349: bpmanip:  [ 0] <TOP>

# 304 trace-shell: error in raytrace pipe

# 1 trace-nest: error in raytrace pipe
WARNING: image platform (linux/amd64) does not match the expected platform (linux/arm64)

for a in off_axis_angle:
    out = subprocess.run(f'marx SourceType=SAOSAC SAOSACFile="sao{a}.fits" DitherModel=FILE DitherFile="sim_asol.fits" OutputDir="saomarx_{a}" ' + 
        # The following line is needed because of https://github.com/Chandra-MARX/marx/issues/83
                         'SourceRA={0} '.format(a / 60.) +
                         default_marx, shell=True, capture_output=True)
    out2 = subprocess.run(f'marx2fits --pixadj=NONE saomarx_{a} saomarx_{a}.fits', 
                          shell=True, capture_output=True)

ecf = [0.5, 0.7, 0.9]
fig, axecf = plt.subplots()
out = np.zeros((len(off_axis_angle), 2, len(ecf)))

for i, a in enumerate(off_axis_angle):
    for prog in ['marx', 'saomarx_']:
        f = f'{prog}{a}.fits'
        with fits.open(f) as hdus:
            header = hdus[1].header
        tab = Table.read(f, hdu=1)
        # The old question: Is an index the center of a pixel of the corner?
        # Differs by 0.5...
        d_ra = (tab['X'] - header['TCRPX9'] - 0.5) * header['TCDLT9'] * 3600.
        # Count from the left or right (RA is reversed on the sky)?
        d_dec = (tab['Y'] - header['TCRPX10'] + 0.5) * header['TCDLT10'] * 3600.

        # The simulation is set up to make this simple: RA=DEC=0
        # so cos(DEC) = 1 and we can approximate with Euklidian distance
        r = radial_distribution(d_ra, d_dec)
        out[i, 0 if prog == 'marx' else 1, :] = np.percentile(r, np.array(ecf) * 100.)

for i, e in enumerate(ecf):
    line, = axecf.plot(off_axis_angle, out[:, 0, i],
                        label='{0} %'.format(e*100))
    axecf.plot(off_axis_angle, out[:, 1, i], color=line.get_color(),
                ls=':')
axecf.legend(loc='lower right')
axecf.set_ylabel('radius [arcsec]')
axecf.set_xlabel('off-axis angle [arcmin]')
axecf.grid()

Radius of encircled count fraction for different off-axis angles. The values shown include the effect of the HRC-I positional uncertainty and the uncertainty in the aspect solution. solid line: marx only, dotted lines: SAOTrace + marx

Both simulations show essentially identical off-axis behavior in this test.