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/15713/primary/pcadf508966864N002_asol1.fits outfile=obs.asp evtfile=download/15713/primary/acisf15713N002_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-7 grating=NONE outfile=obs.arf obsfile=download/15713/primary/acisf15713N002_evt2.fits.gz 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/15713/primary/acisf15713N002_evt2.fits.gz[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.550629973412)" 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/15713/primary/pcadf508966864N002_asol1.fits
marx : Set marx parameters appropriate for observation¶
marx RA_Nom=182.25909891372 Dec_Nom=-51.341662012085 Roll_Nom=33.972028259258 GratingType=NONE ExposureTime=10190.806809961796 DitherModel=FILE DitherFile=download/15713/primary/pcadf508966864N002_asol1.fits TStart=508966864.84771 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
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.25909891372
dec_pnt = -51.341662012085
roll_pnt = 33.972028259258
dither_asol_chandra{ file = "/melkor/d1/guenther/marx/test/testexp55/ACISSPSF/download/15713/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.85771 limit=10190.786809961795 limit_type=sec
marx : Run marx with SAOTrace ray file as input¶
marx RA_Nom=182.25909891372 Dec_Nom=-51.341662012085 Roll_Nom=33.972028259258 GratingType=NONE ExposureTime=10190.806809961796 DitherModel=FILE DitherFile=download/15713/primary/pcadf508966864N002_asol1.fits TStart=508966864.84771 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_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/15713/primary/acisf15713N002_evt2.fits.gz outfile=obs_rand.fits acaofffile=download/15713/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 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')
