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'))
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/13182/primary/pcadf408913236N002_asol1.fits
marx : Set marx parameters appropriate for observation¶
marx RA_Nom=332.17106315813 Dec_Nom=45.73817165926 Roll_Nom=301.41391297181 GratingType=NONE ExposureTime=18161.46342998743 DitherModel=FILE DitherFile=download/13182/primary/pcadf408913236N002_asol1.fits TStart=408913236.56726 SourceRA=332.17 SourceDEC=45.742306 DetectorType=HRC-I DetOffsetX=0.0014262321171452097 DetOffsetZ=-0.0025143152978017724 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.17106315813
dec_pnt = 45.73817165926
roll_pnt = 301.41391297181
dither_asol_chandra{ file = "/melkor/d1/guenther/marx/test/testexp55/HRCIPSF/download/13182/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.57726 limit=18161.44342998743 limit_type=sec
marx : Run marx with SAOTrace ray file as input¶
marx RA_Nom=332.17106315813 Dec_Nom=45.73817165926 Roll_Nom=301.41391297181 GratingType=NONE ExposureTime=18161.46342998743 DitherModel=FILE DitherFile=download/13182/primary/pcadf408913236N002_asol1.fits TStart=408913236.56726 SourceRA=332.17 SourceDEC=45.742306 DetectorType=HRC-I DetOffsetX=0.0014262321171452097 DetOffsetZ=-0.0025143152978017724 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 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
'''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(list(self.figures.keys())[0]),
bbox_inches='tight')
self.save_test_result('marx90', err_marx)
self.save_test_result('saotrace90', err_sao)
