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 on an ACIS-BI chip

CIAO : Obtain spectrum from observed data

asphist infile=download/primary/pcadf508966864N002_asol1.fits outfile=obs.asp evtfile=download/primary/acisf15713N002_evt2.fits clobber=yes
# It would be more accuarte to use mkwarf.
# However, that's slower and for our purposes an approximate solution is good enough.
mkarf detsubsys=ACIS-7 grating=NONE outfile=obs.arf obsfile=download/primary/acisf15713N002_evt2.fits asphistfile=obs.asp sourcepixelx=4096 sourcepixely=4074 engrid="0.3:8.0:0.1" maskfile=NONE pbkfile=NONE dafile=NONE verbose=1 mode=h clobber=yes

# We have to use the same energy binning in energy space here that we used for the arf!
# So, first convert the PI to energy (to a precision that's good enough).
dmtcalc "download/primary/acisf15713N002_evt2.fits[EVENTS][sky=circle(4096, 4074, 4)]" evt2_with_energy.fits  expr="energy=(float)pi*0.0149" clobber=yes
dmextract "evt2_with_energy.fits[bin energy=.3:7.999:0.1]" obs.spec clobber=yes opt=generic
dmcopy "obs.arf[cols energ_lo,energ_hi,specresp]" input_spec.fits clobber=yes
dmpaste input_spec.fits "obs.spec[cols counts]" input_spec_1.fits clobber=yes
dmtcalc input_spec_1.fits input_spec_2.fits  expr="flux=counts/((float)specresp * 1*11543.55063)" clobber=yes

# devide by binwidth to turn the flux into a flux DENSITY for marx
dmtcalc input_spec_2.fits input_spec_3.fits  expr="fluxdens=flux/0.1" clobber=yes
dmcopy "input_spec_3.fits[cols energ_hi,fluxdens]" "input_spec_marx.tbl[opt kernel=text/simple]" clobber=yes
dmcopy "input_spec_3.fits[cols energ_lo,energ_hi,flux]" "input_spec_saotrace.tbl[opt kernel=text/simple]" clobber=yes

# We now use a combination of some of the most obscure UNIX tools to
# bring the SAOTRACE input spectrum into the right format
# see http://cxc.harvard.edu/cal/Hrma/RDB/FileFormat.html

# Remove the leading "#" from the line with the column names
awk '{ sub(/\# ENERG_LO/,"ENERG_LO"); print }' < input_spec_saotrace.tbl > input_spec_saotrace.temp
# Then, replace spaces with tabs
tr ' ' \\t < input_spec_saotrace.temp > input_spec_saotrace.rdb

marx : Set marx parameters appropriate for observation

marx RA_Nom=182.259098914 Dec_Nom=-51.3416620121 Roll_Nom=33.9720282593 GratingType=NONE ExposureTime=10190.80681 DitherModel=FILE DitherFile=download/primary/pcadf508966864N002_asol1.fits TStart=508966864.848 ACIS_Exposure_Time=0.7 SourceRA=182.259167 SourceDEC=-51.344667 DetectorType=ACIS-S DetOffsetX=0.00144494226467 DetOffsetZ=-0.0075429459028 SpectrumType=FILE SpectrumFile=input_spec_marx.tbl OutputDir=marx_only SourceFlux=-1

marx2fits : Use the EDSER subpixel algorithm

marx2fits --pixadj=EDSER marx_only marx_only.fits

Lua input for SAOTrace : SAOTrace input matching observation

ra_pnt = 182.259098914
dec_pnt = -51.3416620121
roll_pnt = 33.9720282593

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

point{ position = { ra = 182.259167,
          dec = -51.344667,
          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=508966864.858  limit=10190.78681 limit_type=sec

marx : Run marx with SAOTrace ray file as input

marx RA_Nom=182.259098914 Dec_Nom=-51.3416620121 Roll_Nom=33.9720282593 GratingType=NONE ExposureTime=10190.80681 DitherModel=FILE DitherFile=download/primary/pcadf508966864N002_asol1.fits TStart=508966864.848 ACIS_Exposure_Time=0.7 SourceRA=182.259167 SourceDEC=-51.344667 DetectorType=ACIS-S DetOffsetX=0.00144494226467 DetOffsetZ=-0.0075429459028 SpectrumType=FILE SpectrumFile=input_spec_marx.tbl OutputDir=marx_saotrace SourceType=SAOSAC SAOSACFile=saotrace.fits SourceFlux=-1

marx2fits : Same setting as above for comparison

marx2fits --pixadj=EDSER marx_saotrace marx_saotrace.fits

marx2fits : USE RANDOMIZE as an alternative to EDSER

marx2fits --pixadj=RANDOMIZE marx_only marx_only_rand.fits

marx2fits : Again, same setting as above

marx2fits --pixadj=RANDOMIZE marx_saotrace marx_saotrace_rand.fits

CIAO : Reprocess the observation with RANDOMIZE

acis_process_events infile=download/primary/acisf15713N002_evt2.fits outfile=obs_rand.fits acaofffile=download/primary/pcadf508966864N002_asol1.fits doevtgrade=no calculate_pi=no pix_adj=RANDOMIZE clobber=yes

Python : Compare radial event distributions

# 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

'''Compare radial event distributions'''
import matplotlib.pyplot as plt

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

evt2file = self.get_data_file('evt2')

fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

for ax, files, title in zip([ax1, ax2],
                            [[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)

    if title == 'EDSER':
        self.save_test_result('marx90', err_marx)
        self.save_test_result('saotrace90', err_sao)

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