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.

Two thermal components on ACIS-I

Sherpa : Generate input spectrum

# set source properties
set_source(xsvapec.a1 + xsvapec.a2)
a1.Ne.frozen = False
a2.Ne = a1.Ne
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004
# get source
my_src = get_source()

# set energy grid
bin_width = 0.01
energies = np.arange(0.15, 12., bin_width)

# evaluate source on energy grid
flux = my_src(energies)

save_arrays("marx_input.tbl", [energies[1:], flux[:-1] / bin_width], ["keV","photons/s/cm**2/keV"], ascii=True, clobber=True)

marx : Point source

marx SourceFlux=-1 SpectrumType=FILE SpectrumFile=marx_input.tbl ExposureTime=30000 TStart=2015.5 OutputDir=marx GratingType=NONE DetectorType=ACIS-I DitherModel=INTERNAL RA_Nom=30 Dec_Nom=40 SourceRA=30 SourceDEC=40

marx2fits : turn into fits file

marx2fits --pixadj=NONE marx marx_evt2.fits

marxasp : Make asol file for CIAO

marxasp MarxDir=marx OutputFile=marx_asol1.fits

CIAO : Extract spectra

Use special settings to account for steps that marx does not include: ACIS QE maps and current (non FEF) rmfs.

phagrid="pi=1:1024:1"

evtfile="marx_evt2.fits"
asolfile="marx_asol1.fits"
phafile="marx_pha.fits"
rmffile="marx_rmf.fits"
arffile="marx_arf.fits"
asphistfile="marx_asp.fits"

asphist infile="$asolfile" outfile="$asphistfile" evtfile="$evtfile" clobber=yes

dmextract infile="$evtfile[sky=circle(4096.5,4096.5,20)][bin $phagrid]" outfile="$phafile" clobber=yes

dmcoords $evtfile asol=$asolfile option=sky x=4096.5 y=4096.5
ccdid=$(pget dmcoords chip_id)
chipx=$(pget dmcoords chipx)
chipy=$(pget dmcoords chipy)

detname="ACIS-$ccdid;UNIFORM;bpmask=0"
grating="NONE"
# For ACIS-I, use engrid="0.3:11.0:0.003". This reflects a limitation of mkrmf.
engrid="0.3:12.0:0.003"

mkarf mirror="hrma" detsubsys="$detname" grating="$grating"         outfile="$arffile" obsfile="$evtfile" engrid="$engrid" asphistfile="$asphistfile"         sourcepixelx=4096.5 sourcepixely=4096.5 maskfile=NONE pbkfile=NONE dafile=NONE clobber=yes verbose=0

fef="$CALDB/data/chandra/acis/fef_pha/acisD2000-01-29fef_phaN0005.fits"

cxfilter="chipx_hi>=$chipx,chipx_lo<=$chipx"
cyfilter="chipy_hi>=$chipy,chipy_lo<=$chipy"
mkrmf infile="$fef[ccd_id=$ccdid,$cxfilter,$cyfilter]" outfile="$rmffile" axis1="energy=$engrid" axis2="$phagrid" thresh=1e-8 clobber=yes verbose=0

CIAO : Use default CIAO tools to extract a spectrum

specextract "marx_evt2.fits[sky=circle(4096.5,4096.5,20)]" spec asp=marx_asp.fits weight=no clobber=yes mskfile=NONE badpixfile=NONE

Sherpa : Compare spectra to input model

load_pha(1, 'marx_pha.fits')
load_rmf(1, 'marx_rmf.fits')
load_arf(1, 'marx_arf.fits')
load_data(2, 'spec_grp.pi')
group_counts(1, 25)
group_counts(2, 25)
ignore(None, 0.3)
ignore(7, None)
# set source properties
set_source(xsvapec.a1 + xsvapec.a2)
a1.Ne.frozen = False
a2.Ne = a1.Ne
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004
a1.kT = 0.7
a1.Ne = 2.0
a1.norm = 0.0002
a2.kT = 2.0
a2.norm = 0.0004

set_source(2, a1 + a2)

# Check how good the input parameters work
statres = get_stat_info()
out = {'chi2marx': statres[0].rstat,
        'chi2default': statres[1].rstat}

# Make the plots
plot_data(1, color='g')
plot_model(1, overplot=True, color='g')
fit(1)
conf(1)
c1 = get_conf_results()
plot_model(1, overplot=True)
plot_data(2, overplot=True, yerrorbars=False, marker='^', color='orange')
fit(2)
conf(2)
c2 = get_conf_results()
plot_model(2, overplot=True, color='orange')
xscale("log")
yscale("log")
savefig('/melkor/d1/guenther/marx/doc/source/tests/figures/SpectrumAPECACISI_spec.png', bbox_inches='tight')

plot_arf(1, color='g')
plot_arf(2, overplot=True, color='orange')
savefig('/melkor/d1/guenther/marx/doc/source/tests/figures/SpectrumAPECACISI_arf.png', bbox_inches='tight')

# compile conf outputs for saving
out['fitmarx'] = {'parnames': c1.parnames,
                   'parmins': c1.parmins,
                   'parvals': c1.parvals,
                   'parmaxes': c1.parmaxes}
out['fitdefault'] = {'parnames': c2.parnames,
                      'parmins': c2.parmins,
                      'parvals': c2.parvals,
                      'parmaxes': c2.parmaxes}

# Calculate maximum relative difference in arf
a1 = get_arf(1)
a2 = get_arf(2)
# bring on same energy grid
engrid = np.arange(0.3, 10., 0.05)
a1 = np.interp(engrid, 0.5 * (a1.energ_lo + a1.energ_hi), a1.get_y())
a2 = np.interp(engrid, 0.5 * (a2.energ_lo + a2.energ_hi) , a2.get_y())

out['arfdiff'] = np.mean(np.abs(a2/a1 -1))

# write numbers to file for later
import json
with open('sherpaout.json', 'w') as f:
    json.dump(out, f)

Python : Compare input with fit model

The input model has several parameters. The calculation here summarizes all that into a single number: The mean deviation between input and fit result as measured in terms of Gaussian sigmas. That is a massive simplification since, the distribution is typically not Gaussian shaped and parameters might not be independent (e.g. there is ambiguity between the absorbing column density and the amount of cool plasma).

# 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.

'''Compare input with fit model

The input model has several parameters. The calculation here summarizes
all that into a single number: The mean deviation between input and fit
result as measured in terms of Gaussian sigmas.
That is a massive simplification since, the distribution is typically
not Gaussian shaped and parameters might not be independent (e.g. there
is ambiguity between the absorbing column density and the amount of
cool plasma).
'''
with open(os.path.join(self.basepath, 'sherpaout.json')) as f:
    sherpaout = json.load(f)

self.save_test_result('chi2marx', sherpaout['chi2marx'])
self.save_test_result('chi2default', sherpaout['chi2default'])
self.save_test_result('arfdiff', sherpaout['arfdiff'])