Note

The following shows the code used to run this marx test. You can inspect it and adapt it to your needs, but you cannot copy and paste it directly because it depends on local $PATH and other environment variables. For example, we use a python function to manage the directory structure for all the images generated by all the tests instead of giving the file name directly to save images.

On-axis PSF for an HRC-I observation

Python : Copy input spectra

Since the HRC has no intrisinc 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.

Spectrum files are part of the marx-test distribution. Copy them from the source code into the right directory.

# Note that this code might not run if you directly copy and paste it:
# - Not all import statements are shown here
# - `self` is a reference to a test instance, which allows access to
#   parameters such as the directory where the test is run etc.

'''Copy input spectra

Since the HRC has no intrisinc 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.

Spectrum files are part of the marx-test distribution.
Copy them from the source code into the right directory.
'''
import shutil
shutil.copy(os.path.join(self.pkg_data, 'ARLac_input_spec_marx.tbl'),
            os.path.join(self.basepath, 'input_spec_marx.tbl'))
shutil.copy(os.path.join(self.pkg_data, 'ARLac_input_spec_saotrace.rdb'),
            os.path.join(self.basepath, 'input_spec_saotrace.rdb'))

marx : Set marx parameters appropriate for observation

marx RA_Nom=332.171063158 Dec_Nom=45.7381716593 Roll_Nom=301.413912972 GratingType=NONE ExposureTime=18161.46343 DitherModel=FILE DitherFile=download/primary/pcadf408913236N002_asol1.fits TStart=408913236.567 SourceRA=332.17 SourceDEC=45.742306 DetectorType=HRC-I DetOffsetX=0.00142623211715 DetOffsetZ=-0.0025143152978 SpectrumType=FILE SpectrumFile=input_spec_marx.tbl OutputDir=marx_only SourceFlux=-1

marx2fits : No EDSER is available for HRC data

marx2fits --pixadj=NONE marx_only marx_only.fits

Lua input for SAOTrace : SAOTrace input matching observation

ra_pnt = 332.171063158
dec_pnt = 45.7381716593
roll_pnt = 301.413912972

dither_asol_chandra{ file = "/melkor/d1/guenther/marx/test/testexp4/HRCIPSF/download/primary/pcadf408913236N002_asol1.fits",
                     ra = ra_pnt, dec = dec_pnt, roll = roll_pnt }

point{ position = { ra = 332.17,
          dec = 45.742306,
          ra_aimpt = ra_pnt,
          dec_aimpt = dec_pnt,
       },
       spectrum = { { file = "input_spec_saotrace.rdb",
                      units = "photons/s/cm2",
                      scale = 1,
                      format = "rdb",
                      emin = "ENERG_LO",
                      emax = "ENERG_HI",
                      flux = "FLUX"} }
    }

SAOTrace :

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.

trace-nest tag=saotrace srcpars=saotrace_source.lua tstart=408913236.577  limit=18161.44343 limit_type=sec

marx : Run marx with SAOTrace ray file as input

marx RA_Nom=332.171063158 Dec_Nom=45.7381716593 Roll_Nom=301.413912972 GratingType=NONE ExposureTime=18161.46343 DitherModel=FILE DitherFile=download/primary/pcadf408913236N002_asol1.fits TStart=408913236.567 SourceRA=332.17 SourceDEC=45.742306 DetectorType=HRC-I DetOffsetX=0.00142623211715 DetOffsetZ=-0.0025143152978 SpectrumType=FILE SpectrumFile=input_spec_marx.tbl OutputDir=marx_saotrace SourceType=SAOSAC SAOSACFile=saotrace.fits SourceFlux=-1

marx2fits : Same settings as the marx2fits run above

marx2fits --pixadj=NONE marx_saotrace marx_saotrace.fits

Python : Extract radial count distribution

# Note that this code might not run if you directly copy and paste it:
# - Not all import statements are shown here
# - `self` is a reference to a test instance, which allows access to
#   parameters such as the directory where the test is run etc.

def filter_events(evt, circle=None, energy=None, time=None):
    '''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:`dm` 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((cc(evt, 'x') - circle[0])**2 + (cc(evt, 'y') - circle[1])**2)
        indcirc = d < circle[2]
    else:
        indcirc = np.ones(len(evt), dtype=bool)

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

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

    return evt[indcirc & inden & indt]

def plot_ecf(ax, files, filterargs, bgfilterargs):
    '''
    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]), **filterargs)
    bkg = filter_events(Table.read(files[0]), **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]), **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]), **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

'''Extract radial count distribution'''
import matplotlib.pyplot as plt
filterargs = {'circle': (self.source['x'], self.source['y'],
                         5 * self.source['r'])
              }

bgfilterargs =  {'circle': (self.source['bg_x'], self.source['bg_y'],
                         5 * self.source['r'])
                }

evt2file = self.get_data_file('evt2')
files = [evt2file, 'marx_only.fits', 'marx_saotrace.fits']

fig = plt.figure()
ax = fig.add_subplot(111)

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

fig.savefig(self.figpath(self.figures.keys()[0]))

self.save_test_result('marx90', err_marx)
self.save_test_result('saotrace90', err_sao)