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.

Wings of the Chandra PSF

CIAO : Obtain spectrum from observed data

asphist infile=download/3662/primary/pcadf141954141N003_asol1.fits outfile=obs.asp evtfile=download/3662/primary/acisf03662N003_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/3662/primary/acisf03662N003_evt2.fits.gz asphistfile=obs.asp sourcepixelx=4174 sourcepixely=4043 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/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS][sky=rotbox(3922.6168,3556.9146,12,512,332.71498)]" 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 * 0.0064*51927.014620006084)" 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/3662/primary/pcadf141954141N003_asol1.fits

marx : Set marx parameters appropriate for observation

marx RA_Nom=254.47042153578 Dec_Nom=35.349687440319 Roll_Nom=207.24753234529 GratingType=NONE ExposureTime=51150.833339989185 DitherModel=FILE DitherFile=download/3662/primary/pcadf141954141N003_asol1.fits TStart=141954141.91872 ACIS_Exposure_Time=3.1 SourceRA=254.457583 SourceDEC=35.342433 DetectorType=ACIS-S DetOffsetX=0.0014449422646707344 DetOffsetZ=-8.777474717742791 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 = 254.47042153578
dec_pnt = 35.349687440319
roll_pnt = 207.24753234529

dither_asol_chandra{ file = "/melkor/d1/guenther/marx/test/testexp55/Wings/download/3662/primary/pcadf141954141N003_asol1.fits",
                     ra = ra_pnt, dec = dec_pnt, roll = roll_pnt }

point{ position = { ra = 254.457583,
          dec = 35.342433,
          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=141954141.92872  limit=51150.81333998919 limit_type=sec

marx : Run marx with SAOTrace ray file as input

marx RA_Nom=254.47042153578 Dec_Nom=35.349687440319 Roll_Nom=207.24753234529 GratingType=NONE ExposureTime=51150.833339989185 DitherModel=FILE DitherFile=download/3662/primary/pcadf141954141N003_asol1.fits TStart=141954141.91872 ACIS_Exposure_Time=3.1 SourceRA=254.457583 SourceDEC=35.342433 DetectorType=ACIS-S DetOffsetX=0.0014449422646707344 DetOffsetZ=-8.777474717742791 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

CIAO : Make image and run cell detect to find exclusion regions for extraction.

dmcopy "download/3662/primary/acisf03662N003_evt2.fits.gz[bin x=3700:4300:1,y=3700:4300:1][opt type=i4]" im.fits option=image clobber=yes
mkpsfmap im.fits psf.map 1.4 ecf=0.5
celldetect im.fits im_src.fits psffile=psf.map clobber=yes
asphist download/3662/primary/pcadf141954141N003_asol1.fits asphist_7.fits evtfile="download/3662/primary/acisf03662N003_evt2.fits.gz[ccd_id=7]" clobber=yes

CIAO : Create exposure maps

Also, make exposure map for S-3 chip to correct for dither etc. Mirror effective area is not included in that calculation, because all the photons we care for come from the same source, not from different positions on the sky, thus the mirror effective area is the constant for the source we look at.

The instrument map depends strongly on the energy. We will thus pass in the source spectrum. This gives us one map appropriately scaled for the whole spectrum. However, in later analysis, we want to split the data in narrower bands. When looking at that accurately, we need to have one map for each band.

Also, the observed data contains bad pixels, while the marx simulated data does not. This means that we need to make a second set of maps with no bad pixels in them.

mkinstmap instmap_0.4.fits monoenergy=0.45 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.4.fits instmap_0.4.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.4_marx.fits monoenergy=0.45 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.4_marx.fits instmap_0.4_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.5.fits monoenergy=0.55 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.5.fits instmap_0.5.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.5_marx.fits monoenergy=0.55 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.5_marx.fits instmap_0.5_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.6.fits monoenergy=0.6499999999999999 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.6.fits instmap_0.6.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.6_marx.fits monoenergy=0.6499999999999999 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.6_marx.fits instmap_0.6_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.7.fits monoenergy=0.8 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.7.fits instmap_0.7.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.7_marx.fits monoenergy=0.8 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.7_marx.fits instmap_0.7_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.9.fits monoenergy=0.95 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.9.fits instmap_0.9.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_0.9_marx.fits monoenergy=0.95 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_0.9_marx.fits instmap_0.9_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.0.fits monoenergy=1.1 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.0.fits instmap_1.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.0_marx.fits monoenergy=1.1 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.0_marx.fits instmap_1.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.2.fits monoenergy=1.2999999999999998 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.2.fits instmap_1.2.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.2_marx.fits monoenergy=1.2999999999999998 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.2_marx.fits instmap_1.2_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.4.fits monoenergy=1.5 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.4.fits instmap_1.4.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.4_marx.fits monoenergy=1.5 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.4_marx.fits instmap_1.4_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.6.fits monoenergy=1.7000000000000002 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.6.fits instmap_1.6.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.6_marx.fits monoenergy=1.7000000000000002 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.6_marx.fits instmap_1.6_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.8.fits monoenergy=1.9 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.8.fits instmap_1.8.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_1.8_marx.fits monoenergy=1.9 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_1.8_marx.fits instmap_1.8_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.0.fits monoenergy=2.1 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.0.fits instmap_2.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.0_marx.fits monoenergy=2.1 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.0_marx.fits instmap_2.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.2.fits monoenergy=2.3 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.2.fits instmap_2.2.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.2_marx.fits monoenergy=2.3 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.2_marx.fits instmap_2.2_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.4.fits monoenergy=2.5 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.4.fits instmap_2.4.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.4_marx.fits monoenergy=2.5 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.4_marx.fits instmap_2.4_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.6.fits monoenergy=2.7 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.6.fits instmap_2.6.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.6_marx.fits monoenergy=2.7 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.6_marx.fits instmap_2.6_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.8.fits monoenergy=2.9 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.8.fits instmap_2.8.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_2.8_marx.fits monoenergy=2.9 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_2.8_marx.fits instmap_2.8_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_3.0.fits monoenergy=3.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_3.0.fits instmap_3.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_3.0_marx.fits monoenergy=3.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_3.0_marx.fits instmap_3.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_3.5.fits monoenergy=3.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_3.5.fits instmap_3.5.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_3.5_marx.fits monoenergy=3.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_3.5_marx.fits instmap_3.5_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_4.0.fits monoenergy=4.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_4.0.fits instmap_4.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_4.0_marx.fits monoenergy=4.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_4.0_marx.fits instmap_4.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_4.5.fits monoenergy=4.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_4.5.fits instmap_4.5.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_4.5_marx.fits monoenergy=4.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_4.5_marx.fits instmap_4.5_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_5.0.fits monoenergy=5.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_5.0.fits instmap_5.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_5.0_marx.fits monoenergy=5.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_5.0_marx.fits instmap_5.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_5.5.fits monoenergy=5.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_5.5.fits instmap_5.5.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_5.5_marx.fits monoenergy=5.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_5.5_marx.fits instmap_5.5_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_6.0.fits monoenergy=6.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_6.0.fits instmap_6.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_6.0_marx.fits monoenergy=6.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_6.0_marx.fits instmap_6.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_6.5.fits monoenergy=6.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_6.5.fits instmap_6.5.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_6.5_marx.fits monoenergy=6.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_6.5_marx.fits instmap_6.5_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_7.0.fits monoenergy=7.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_7.0.fits instmap_7.0.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_7.0_marx.fits monoenergy=7.25 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_7.0_marx.fits instmap_7.0_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_7.5.fits monoenergy=7.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_7.5.fits instmap_7.5.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no
mkinstmap instmap_7.5_marx.fits monoenergy=7.75 pixelgrid="1:1024:#1024,1:1024:#1024" obsfile="download/3662/primary/acisf03662N003_evt2.fits.gz[EVENTS]" dafile=NONE detsubsys="ACIS-S3;BPMASK=0x00000" mirror="HRMA;AREA=1" grating=NONE maskfile=NONE spectrumfil=NONE clobber=yes
mkexpmap asphist_7.fits expmap_7.5_marx.fits instmap_7.5_marx.fits xygrid="2975:4400:#512,3250:4650:#512" clobber=yes useavgaspect=no

Python : Clean bogus sources from src file and write regions

celldetect always detects plenty of sources on the read-out streak and around the piled-up regions. For processing speed and clarity we want to remove those bogus detections from the source list.

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

'''Clean bogus sources from src file and write regions

celldetect always detects plenty of sources on the read-out streak and
around the piled-up regions.
For processing speed and clarity we want to remove those bogus detections
from the source list.
'''
import numpy as np
from astropy.table import Table
src = Table.read(os.path.join(self.basepath, 'im_src.fits'), hdu=1)
# Get sources within 50 pix of Her X-1. We don't look there anyway because of pile-up
src['d_herx1'] = np.sqrt((src['X'] - self.source['x'])**2 +
                         (src['Y'] - self.source['y'])**2)
d_src = src['d_herx1'] < 50.
# Get sources on the read-out streak
# 2 points on streak
p1 = np.array([3772., 3264.])
p2 = np.array([4251., 4193.])
# n is normal to read-out streak
n = np.array([1., - (p2[0] - p1[0]) / (p2[1] - p1[1])])
n = n / np.linalg.norm(n)
x = np.vstack([src['X'], src['Y']])
d_line = np.abs(np.dot(n, (p1 - x.T).T)) < 2.
src = src[~d_src & ~d_line]

# box for readout streak
streak = 'rotbox(4010.0285,3725.3416,18.240189,1085.0354,332.74569)'
# Assemble regions for extraction:
n_rings = 50
r_rings = np.logspace(np.log10(7.), np.log10(500.), n_rings + 1) / 0.492
with open(os.path.join(self.basepath, 'annuli.stk'), 'w') as f:
    for i in range(n_rings):
        reg = 'annulus({0}, {1}, {2}, {3})'.format(self.source['x'], self.source['y'],
                                                   r_rings[i], r_rings[i + 1])
        reg += '-' + streak
        # find sources that overlap this ring
        ind1 = src['d_herx1'] + src['R'][:, 0] > r_rings[i]
        ind2 = src['d_herx1'] - src['R'][:, 0] < r_rings[i + 1]
        src_intersect = src[ind1 & ind2]
        for s in src_intersect:
            reg += '-ellipse({0}, {1}, {2}, {3}, {4})'.format(s['X'], s['Y'],
                                                                s['R'][0], s['R'][1],
                                                                s['ROTANG'])
        f.write(reg + '\n')

CIAO : Extract data in energy bins

dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=400.0:500.0][bin sky=@annuli.stk]" profile_0.4_obs.fits exp=expmap_0.4.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=500.0:600.0][bin sky=@annuli.stk]" profile_0.5_obs.fits exp=expmap_0.5.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=600.0:700.0][bin sky=@annuli.stk]" profile_0.6_obs.fits exp=expmap_0.6.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=700.0:900.0][bin sky=@annuli.stk]" profile_0.7_obs.fits exp=expmap_0.7.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=900.0:1000.0][bin sky=@annuli.stk]" profile_0.9_obs.fits exp=expmap_0.9.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=1000.0:1200.0][bin sky=@annuli.stk]" profile_1.0_obs.fits exp=expmap_1.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=1200.0:1400.0][bin sky=@annuli.stk]" profile_1.2_obs.fits exp=expmap_1.2.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=1400.0:1600.0][bin sky=@annuli.stk]" profile_1.4_obs.fits exp=expmap_1.4.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=1600.0:1800.0][bin sky=@annuli.stk]" profile_1.6_obs.fits exp=expmap_1.6.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=1800.0:2000.0][bin sky=@annuli.stk]" profile_1.8_obs.fits exp=expmap_1.8.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=2000.0:2200.0][bin sky=@annuli.stk]" profile_2.0_obs.fits exp=expmap_2.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=2200.0:2400.0][bin sky=@annuli.stk]" profile_2.2_obs.fits exp=expmap_2.2.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=2400.0:2600.0][bin sky=@annuli.stk]" profile_2.4_obs.fits exp=expmap_2.4.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=2600.0:2800.0][bin sky=@annuli.stk]" profile_2.6_obs.fits exp=expmap_2.6.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=2800.0:3000.0][bin sky=@annuli.stk]" profile_2.8_obs.fits exp=expmap_2.8.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=3000.0:3500.0][bin sky=@annuli.stk]" profile_3.0_obs.fits exp=expmap_3.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=3500.0:4000.0][bin sky=@annuli.stk]" profile_3.5_obs.fits exp=expmap_3.5.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=4000.0:4500.0][bin sky=@annuli.stk]" profile_4.0_obs.fits exp=expmap_4.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=4500.0:5000.0][bin sky=@annuli.stk]" profile_4.5_obs.fits exp=expmap_4.5.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=5000.0:5500.0][bin sky=@annuli.stk]" profile_5.0_obs.fits exp=expmap_5.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=5500.0:6000.0][bin sky=@annuli.stk]" profile_5.5_obs.fits exp=expmap_5.5.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=6000.0:6500.0][bin sky=@annuli.stk]" profile_6.0_obs.fits exp=expmap_6.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=6500.0:7000.0][bin sky=@annuli.stk]" profile_6.5_obs.fits exp=expmap_6.5.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=7000.0:7500.0][bin sky=@annuli.stk]" profile_7.0_obs.fits exp=expmap_7.0.fits opt=generic2 clobber=yes
dmextract "download/3662/primary/acisf03662N003_evt2.fits.gz[energy=7500.0:8000.0][bin sky=@annuli.stk]" profile_7.5_obs.fits exp=expmap_7.5.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=400.0:500.0][bin sky=@annuli.stk]" profile_0.4_marx.fits exp=expmap_0.4_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=500.0:600.0][bin sky=@annuli.stk]" profile_0.5_marx.fits exp=expmap_0.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=600.0:700.0][bin sky=@annuli.stk]" profile_0.6_marx.fits exp=expmap_0.6_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=700.0:900.0][bin sky=@annuli.stk]" profile_0.7_marx.fits exp=expmap_0.7_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=900.0:1000.0][bin sky=@annuli.stk]" profile_0.9_marx.fits exp=expmap_0.9_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=1000.0:1200.0][bin sky=@annuli.stk]" profile_1.0_marx.fits exp=expmap_1.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=1200.0:1400.0][bin sky=@annuli.stk]" profile_1.2_marx.fits exp=expmap_1.2_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=1400.0:1600.0][bin sky=@annuli.stk]" profile_1.4_marx.fits exp=expmap_1.4_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=1600.0:1800.0][bin sky=@annuli.stk]" profile_1.6_marx.fits exp=expmap_1.6_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=1800.0:2000.0][bin sky=@annuli.stk]" profile_1.8_marx.fits exp=expmap_1.8_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=2000.0:2200.0][bin sky=@annuli.stk]" profile_2.0_marx.fits exp=expmap_2.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=2200.0:2400.0][bin sky=@annuli.stk]" profile_2.2_marx.fits exp=expmap_2.2_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=2400.0:2600.0][bin sky=@annuli.stk]" profile_2.4_marx.fits exp=expmap_2.4_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=2600.0:2800.0][bin sky=@annuli.stk]" profile_2.6_marx.fits exp=expmap_2.6_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=2800.0:3000.0][bin sky=@annuli.stk]" profile_2.8_marx.fits exp=expmap_2.8_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=3000.0:3500.0][bin sky=@annuli.stk]" profile_3.0_marx.fits exp=expmap_3.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=3500.0:4000.0][bin sky=@annuli.stk]" profile_3.5_marx.fits exp=expmap_3.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=4000.0:4500.0][bin sky=@annuli.stk]" profile_4.0_marx.fits exp=expmap_4.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=4500.0:5000.0][bin sky=@annuli.stk]" profile_4.5_marx.fits exp=expmap_4.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=5000.0:5500.0][bin sky=@annuli.stk]" profile_5.0_marx.fits exp=expmap_5.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=5500.0:6000.0][bin sky=@annuli.stk]" profile_5.5_marx.fits exp=expmap_5.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=6000.0:6500.0][bin sky=@annuli.stk]" profile_6.0_marx.fits exp=expmap_6.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=6500.0:7000.0][bin sky=@annuli.stk]" profile_6.5_marx.fits exp=expmap_6.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=7000.0:7500.0][bin sky=@annuli.stk]" profile_7.0_marx.fits exp=expmap_7.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_only.fits[energy=7500.0:8000.0][bin sky=@annuli.stk]" profile_7.5_marx.fits exp=expmap_7.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=400.0:500.0][bin sky=@annuli.stk]" profile_0.4_saotrace.fits exp=expmap_0.4_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=500.0:600.0][bin sky=@annuli.stk]" profile_0.5_saotrace.fits exp=expmap_0.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=600.0:700.0][bin sky=@annuli.stk]" profile_0.6_saotrace.fits exp=expmap_0.6_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=700.0:900.0][bin sky=@annuli.stk]" profile_0.7_saotrace.fits exp=expmap_0.7_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=900.0:1000.0][bin sky=@annuli.stk]" profile_0.9_saotrace.fits exp=expmap_0.9_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=1000.0:1200.0][bin sky=@annuli.stk]" profile_1.0_saotrace.fits exp=expmap_1.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=1200.0:1400.0][bin sky=@annuli.stk]" profile_1.2_saotrace.fits exp=expmap_1.2_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=1400.0:1600.0][bin sky=@annuli.stk]" profile_1.4_saotrace.fits exp=expmap_1.4_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=1600.0:1800.0][bin sky=@annuli.stk]" profile_1.6_saotrace.fits exp=expmap_1.6_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=1800.0:2000.0][bin sky=@annuli.stk]" profile_1.8_saotrace.fits exp=expmap_1.8_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=2000.0:2200.0][bin sky=@annuli.stk]" profile_2.0_saotrace.fits exp=expmap_2.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=2200.0:2400.0][bin sky=@annuli.stk]" profile_2.2_saotrace.fits exp=expmap_2.2_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=2400.0:2600.0][bin sky=@annuli.stk]" profile_2.4_saotrace.fits exp=expmap_2.4_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=2600.0:2800.0][bin sky=@annuli.stk]" profile_2.6_saotrace.fits exp=expmap_2.6_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=2800.0:3000.0][bin sky=@annuli.stk]" profile_2.8_saotrace.fits exp=expmap_2.8_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=3000.0:3500.0][bin sky=@annuli.stk]" profile_3.0_saotrace.fits exp=expmap_3.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=3500.0:4000.0][bin sky=@annuli.stk]" profile_3.5_saotrace.fits exp=expmap_3.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=4000.0:4500.0][bin sky=@annuli.stk]" profile_4.0_saotrace.fits exp=expmap_4.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=4500.0:5000.0][bin sky=@annuli.stk]" profile_4.5_saotrace.fits exp=expmap_4.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=5000.0:5500.0][bin sky=@annuli.stk]" profile_5.0_saotrace.fits exp=expmap_5.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=5500.0:6000.0][bin sky=@annuli.stk]" profile_5.5_saotrace.fits exp=expmap_5.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=6000.0:6500.0][bin sky=@annuli.stk]" profile_6.0_saotrace.fits exp=expmap_6.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=6500.0:7000.0][bin sky=@annuli.stk]" profile_6.5_saotrace.fits exp=expmap_6.5_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=7000.0:7500.0][bin sky=@annuli.stk]" profile_7.0_saotrace.fits exp=expmap_7.0_marx.fits opt=generic2 clobber=yes
dmextract "marx_saotrace.fits[energy=7500.0:8000.0][bin sky=@annuli.stk]" profile_7.5_saotrace.fits exp=expmap_7.5_marx.fits opt=generic2 clobber=yes

Python : Analyze and plot results - Plot 1

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

'''Analyze and plot results - Plot 1'''
import numpy as np
from matplotlib import pyplot as plt
from astropy.table import Table
from sherpa.data import Data1D
from sherpa.models.basic import Const1D, PowLaw1D
from sherpa.stats import Chi2
from sherpa.optmethods import LevMar
from sherpa.fit import Fit

# set up Sherpa fitting
const1d = Const1D()
const1d.c0.min = 0
pow1d = PowLaw1D()
stat = Chi2()
opt = LevMar()

fig = plt.figure(figsize=(10, 8))
energy_index = [1, 8, 18, 24]
fullname = {'obs': 'Observation', 'marx': 'MARX',
            'saotrace': 'SAOTrace + MARX'}
for i in range(len(energy_index)):
    e = self.energy_bins[energy_index[i]]
    e_upper = self.energy_bins[energy_index[i] + 1]
    ax = fig.add_subplot(2, 2, i + 1)
    for j, name in enumerate(['obs', 'marx', 'saotrace']):
        tab = Table.read(f'profile_{e:3.1f}_{name}.fits', hdu=1)
        r = np.mean(tab['R'], axis=1)
        r_arcsec = r * 0.492
        indfit = r_arcsec > 15.
        plots = ax.errorbar(r_arcsec, tab['SUR_FLUX'],
                            yerr=tab['SUR_FLUX_ERR'],
                            label='__no_legend__')

        if name == 'obs':
            mymod = const1d + pow1d
            const1d.c0.val = tab['SUR_FLUX'][-1]
        else:
            mymod = pow1d
        data = Data1D('', r[indfit], tab['SUR_FLUX'][indfit],
                      staterror=tab['SUR_FLUX_ERR'][indfit])
        gfit = Fit(data, mymod, stat=stat, method=opt)
        gfit.fit()

        ax.plot(r_arcsec, mymod(r), label=fullname[name],
                color=plots[0].get_color())
        ax.set_title('{0:3.1f} - {1:3.1f} keV'.format(e, e_upper))
        ax.loglog()
        ax.set_xlim(6, 550)
        if i in [0, 2]:
            ax.set_ylabel('flux per area')
        if i in [2, 3]:
            ax.set_xlabel('radius [arcsec]')

ax.legend(fontsize='small', loc='lower left')
fig.savefig(self.figpath(list(self.figures.keys())[0]),
            bbox_inches='tight')

Python : Analyze and plot results - Plot 2

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

'''Analyze and plot results - Plot 2'''
import numpy as np
from matplotlib import pyplot as plt
from astropy.table import Table
from astropy.modeling import models
from sherpa.data import Data1D
from sherpa.models.basic import Const1D, PowLaw1D
from sherpa.stats import Chi2
from sherpa.optmethods import LevMar
from sherpa.fit import Fit

# set up Sherpa fitting
const1d = Const1D()
const1d.c0.min = 0
pow1d = PowLaw1D()
stat = Chi2()
opt = LevMar()

slopes = np.zeros((len(self.energy_bins) - 1, 3))
for i in range(len(self.energy_bins) - 1):
    e = self.energy_bins[i]
    for j, name in enumerate(['obs', 'marx', 'saotrace']):
        tab = Table.read(f'profile_{e:3.1f}_{name}.fits', hdu=1)
        r = np.mean(tab['R'], axis=1)
        r_arcsec = r * 0.492
        indfit = r_arcsec > 15.
        if name == 'obs':
            mymod = const1d + pow1d
            const1d.c0.val = tab['SUR_FLUX'][-1]
        else:
            mymod = pow1d
        data = Data1D('', r[indfit], tab['SUR_FLUX'][indfit],
                      staterror=tab['SUR_FLUX_ERR'][indfit])
        gfit = Fit(data, mymod, stat=stat, method=opt)
        gfit.fit()

        if name == 'obs':
            slopes[i, j] = pow1d.gamma.val
        else:
            slopes[i, j] = pow1d.gamma.val
        print(gfit)

fig = plt.figure()
ax = fig.add_subplot(111)
e_mid = 0.5 * (self.energy_bins[:-1] + self.energy_bins[1:])
for i, name in enumerate(['Observation', 'MARX', 'SAOTrace + MARX']):
    ax.plot(e_mid, slopes[:, i], label=name)
    ax.legend()
    ax.set_xlabel('energy [keV]')
    ax.set_ylabel('slope of powerlaw')
fig.savefig(self.figpath(list(self.figures.keys())[1]),
            bbox_inches='tight')