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-FI chip¶
CIAO : Obtain spectrum from observed data¶
asphist infile=download/4469/primary/pcadf188536756N003_asol1.fits outfile=obs.asp evtfile=download/4469/primary/acisf04469N003_evt2.fits.gz 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-3 grating=NONE outfile=obs.arf obsfile=download/4469/primary/acisf04469N003_evt2.fits.gz asphistfile=obs.asp sourcepixelx=4140 sourcepixely=4102 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/4469/primary/acisf04469N003_evt2.fits.gz[EVENTS][sky=circle(4140, 4102, 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*100755.4543800056)" 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
shell : Unzip fits file.¶
MARX cannot read zipped fits files, so we need to unzip the .fits.gz asol files that we downloaded from the archive. On the other hand, CIAO tools work on both zipped or unzipped files, so there is no need to unzip all of them, just the files that MARX reads as input.
gunzip -f download/4469/primary/pcadf188536756N003_asol1.fits
marx : Set marx parameters appropriate for observation¶
marx RA_Nom=109.68346446157 Dec_Nom=-24.95518362722 Roll_Nom=22.312764819377 GratingType=NONE ExposureTime=99150.81681001186 DitherModel=FILE DitherFile=download/4469/primary/pcadf188536756N003_asol1.fits TStart=188536756.92716 ACIS_Exposure_Time=3.2 SourceRA=109.678333 SourceDEC=-24.955139 DetectorType=ACIS-I DetOffsetX=0.0014398546217030406 DetOffsetZ=0.0050286306017994775 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 = 109.68346446157
dec_pnt = -24.95518362722
roll_pnt = 22.312764819377
dither_asol_chandra{ file = "/melkor/d1/guenther/marx/test/testexp55/ACISIPSF/download/4469/primary/pcadf188536756N003_asol1.fits",
ra = ra_pnt, dec = dec_pnt, roll = roll_pnt }
point{ position = { ra = 109.678333,
dec = -24.955139,
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=188536756.93716 limit=99150.79681001186 limit_type=sec
marx : Run marx with SAOTrace ray file as input¶
marx RA_Nom=109.68346446157 Dec_Nom=-24.95518362722 Roll_Nom=22.312764819377 GratingType=NONE ExposureTime=99150.81681001186 DitherModel=FILE DitherFile=download/4469/primary/pcadf188536756N003_asol1.fits TStart=188536756.92716 ACIS_Exposure_Time=3.2 SourceRA=109.678333 SourceDEC=-24.955139 DetectorType=ACIS-I DetOffsetX=0.0014398546217030406 DetOffsetZ=0.0050286306017994775 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/4469/primary/acisf04469N003_evt2.fits.gz outfile=obs_rand.fits acaofffile=download/4469/primary/pcadf188536756N003_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 radial_distribution(x, y):
'''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 - astropy.stats.sigma_clipped_stats(x, sigma=1.5, maxiters=20)[0]
centy = y - astropy.stats.sigma_clipped_stats(y, sigma=1.5, maxiters=20)[0]
xy = np.vstack([centx, centy])
return np.linalg.norm(xy, axis=0)
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(list(self.figures.keys())[0]),
bbox_inches='tight')
