Using srcflux plugins

Reprocess data to apply latest calibrations

We begin by applying the latest calibrations by running chandra_repro

unix% chandra_repro 18915,20906 outdir= clob+

Here chandra_repro is given a stack of directories to process. With outdir=, a blank string, a repro directory will be created under each individual input directory.

Create initial images

The next step is to create images, exposure maps, and PSF maps for each individual observations; these are needed to detect sources using wavdetect. For simplicity we are going to restrict data to just the aimpoint chip, ACIS-7.

unix% /bin/ls */repro/*evt2.fits > evt.lis
unix% merge_obs @evt.lis"[ccd_id=7]" out=merge/out bin=1 band=broad psfecf=0.5 mode=h clob+

The output from merge_obs includes the images, exposure maps, PSF maps, and the event files reprojected to a common tangent point for each individual observation as well as the combined products.

For more details see Using merge_obs to combine observations and create exposure-corrected images.

Detect sources in each individual observation

Next we will detect sources in each individaul observation using wavdetect

unix% for obsid in 18915 20906
do
  punlearn wavdetect wrecon wtransform
  pset wavdetect \
     infile=merge/out_${obsid}_broad_thresh.img \
     expfile=merge/out_${obsid}_broad_thresh.expmap \
     psffile=merge/out_${obsid}_broad_thresh.psfmap \
     outfile=wav_${obsid}.src \
     scales="1.4 2 4" \
     scell=wav.cell defn=wav.nbg imagef=wav.recon clob+ 
  wavdetect mode=h
done

For more details see the Running wavdetect thread.

Compute fine astrometric offsets

Users analyzing a single observations can skip to the next section.

We strongly recommend that users inspect the output from wavdetect and possibly filter their source lists before using them for any analysis, including cross-matching for astrometric correction. We did not find any spurious sources that would affect the fine astrometric correction in this example.

The next step is to cross-match the sources. For this example thread we simply choose to match the observations to OBS_ID 20906. This was an arbitrary choice. Users may want to cross-match with an external catalog (such as Gaia) or develop their own heuristics for picking a reference observation (eg longest observation, most sources, etc.).

unix% ref=20906
unix% for obsid in 18915 20906
do
  punlearn wcs_match
  wcs_match infile=wav_${obsid}.src ref=wav_${ref}.src out=xform_${obsid}.fits method=trans radius=4 clob+ \
          wcs=merge/out_${obsid}_reproj_evt.fits
  dmlist xform_${obsid}.fits"[cols t1,t2]" data,clean
done

The output looks like

#  t1                   t2
    -0.16420759778384    -0.59824529646612
#  t1                   t2
                    0                    0

t1 and t2 represent the X and Y offsets, in physical pixels (0.492") between the observations. We can see that there is over a half pixel shift between the two observations which can be significant especially for sources near the optical axis.

For more details and options see the Correcting Absolute Astrometry thread and the Calculate source count rates and fluxes for combined datasets thread.

Apply astrometric corrections

The next step is to apply the fine astrometric correction to each observation's event file and aspect solution files. We also copy the necessary auxiliary files: the bad pixel file and detector mask files into the same directory so that the analysis scripts can locate them automatically.

unix% mkdir -p fine_astro

unix% for obsid in 18915 20906
do
  cp merge/out_${obsid}_reproj_evt.fits fine_astro/out_${obsid}_reproj_fa_evt.fits 
  wcs_update fine_astro/out_${obsid}_reproj_fa_evt.fits trans=xform_${obsid}.fits outfile= \
    wcs=merge/out_${obsid}_reproj_evt.fits
  wcs_update ${obsid}/repro/pcadf${obsid}_00?N00?_asol1.fits \
    outfile=fine_astro/out_${obsid}_fa_asol.fits \
    trans=xform_${obsid}.fits wcs=merge/out_${obsid}_reproj_evt.fits clob+
  dmhedit fine_astro/out_${obsid}_reproj_fa_evt.fits file= op=add key=ASOLFILE value=out_${obsid}_fa_asol.fits

  cp ${obsid}/repro/*bpix1.fits fine_astro/
  cp ${obsid}/repro/*msk1.fits fine_astro/
done

The event file is copied because wcs_update updates it in-place (only header keywords are modified), whereas a new aspect solution file is created since each row is updated.

Again for more details see the Correcting Absolute Astrometry thread.

Create final merged products and detect sources

The final data preparation step then is to re-merge the astrometric corrected data and run wavdetect on the co-aligned image.

unix% /bin/ls fine_astro/*evt.fits > fa_evt2.lis
unix% merge_obs @fa_evt2.lis merge_fa/out bin=1 band=broad psfecf=0.5 mode=h clob+
unix% punlearn wavdetect wrecon wtransform
unix% pset wavdetect \
  infile=merge_fa/out_broad_thresh.img \
  expfile=merge_fa/out_broad_thresh.expmap \
  psffile=merge_fa/out_broad_thresh.psfmap \
  outfile=wav_merge.src scell=wav.cell defn=wav.nbkg imagef=wav.recon \
  scales="1.4 2 4" clob+
unix% wavdetect mode=h

unix% dmlist wav_merge.src counts
6

The combined image with the source detections are shown in Figure 1.

Figure 1: Broad band counts image for combined dataset

[Version: full-size]

Figure 1: Broad band counts image for combined dataset

The broad band image for the merged dataset. Both observations were performed using a 256 row (one-quarter chip) subarray. The wavdetect identified source regions are the overlaid green ellipses.

Introduction to Plugins

Plugins allow users to write code to perform their own custom analysis and have srcflux automatically run that code and collect the outputs. Plugins were added in CIAO 4.15. srcflux supports two plugins: one to be run for each individual observation, and one to be run on the merged results.

There are a few simple requirements when writing a plugin.

• The plugin is written in Python.
• For the per-observation plugin, the Python script must define a routine named srcflux_obsid_plugin. The routine is run separately for each observation, for each source, and for each energy band. The routine's parameters are:
  • infile : the name of the input event file.
  • outroot : the name of the output root for this source number.
  • band : the name of the energy band.
  • elo : the low energy bound for this energy band, in keV.
  • ehi : the high energy bound for this energy band, in keV.
  • src_num : the number of the source (1 to N) being processed.

  For example:

  def srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
      return []
  

• For the merged results plugin, the Python script must define a routine named srcflux_merge_plugin. The parameters are the same as srcflux_obsid_plugin, except that the infile is a list of all the input event file names.

  def srcflux_merge_plugin(infile, outroot, band, elo, ehi, src_num):
      return []
  

• The plugin must always return a list of values to be added to the output .flux file. If there are no values then it still must return and empty list, ie []. The list must be returned even if there is an Exception raised.
• Each item in the returned list must have the following attributes: name, units, value, and description. This may be a NamedTuple, a dataclass object, or any other object type will work. The output column name must not already exist in the .flux file and must adhere to FITS standards (letters, numbers, hyphens and underscores only).
• The plugin can create files, but it must clean up any temporary files itself.

Sample Plugins

We have provided several sample plugins in the $ASCDS_CONTRIB/data/sample_srcflux_plugins.py file.

unix% cat $ASCDS_CONTRIB/data/sample_srcflux_plugins.py
"Sample plugin file for srcflux"


"""
These are the files that will be available for the plugin to use.
For a single input event file, the ${outroot} will be the same as 
the outroot parameter.  For multiple input event files, 
the obi number is appended to the outroot parameter.

    The main output files w/ all the src properties:
        ${outroot}_${band}.flux

    Files used for spectral fits:
        ${outroot}.arf
        ${outroot}.pi
        ${outroot}.rmf
        ${outroot}_bkg.pi
        ${outroot}_bkg.arf   (only if bkgresp=yes)
        ${outroot}_bkg.rmf   (only if bkgresp=yes)
        ${outroot}_grp.pi
        ${outroot}_nopsf.arf

    Region files:
        ${outroot}_bkgreg.fits
        ${outroot}_srcreg.fits

    Light curves
        ${outroot}_${band}.gllc
        ${outroot}_${band}.lc

    Postage stamp images
        ${outroot}_${band}_flux.img
        ${outroot}_${band}_thresh.expmap
        ${outroot}_${band}_thresh.img

    PSF (only if psfmethod=marx)
        ${outroot}_${band}_projrays.fits
        ${outroot}_${band}.psf

    aprates probability density function
        ${outroot}_${band}_rates.prob

For multiple input event files, these files will be available to the 
plugin. The outroot is the same as the outroot parameter.

    The main output files w/ all the src properties:
        ${outroot}_${band}.flux

    Files used for spectral fits:
        ${outroot}_src.arf   \
        ${outroot}_src.pi     - Note: merge file names have _src 
        ${outroot}_src.rmf   /
        ${outroot}_bkg.pi
        ${outroot}_bkg.arf   (only if bkgresp=yes)
        ${outroot}_bkg.rmf   (only if bkgresp=yes)

"""
...

In this file users will find plugins to

• def srcflux_merge_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin creates merged images.
• def spectral_fit_sample_srcflux_merge_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin performs a spectral fit using the combined data products.
• def hardness_ratio_sample_srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin computes hardness ratios for individual observations.
• def radial_profile_example_srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin create a radial profile.
• def arestore_sample_srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin performs a blind deconvolution.
• def srcextent_sample_srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin tries to estimate the extent of the source compared to the PSF.
• def spectral_fit_sample_srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin fits the spectra for each individual observation.
• def aplimits_srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
  This plugin computes upper limits and false source probabilities.

Our plugin: Spectral fitting an absorbed APEC model

Based on the above examples we have written the following plugin. Users may download the file here: my_plugin.py.

unix% cat -n my_plugin.py
     1	
     2	'An example plugin for the scrflux script'
     3	
     4	import os
     5	import numpy as np
     6	import matplotlib
     7	import matplotlib.pylab as plt
     8	from sherpa.astro.ui import *
     9	
    10	from collections import namedtuple
    11	ReturnValue = namedtuple('ReturnValue', 'name value units description')
    12	
    13	matplotlib.use("agg")
    14	
    15	def srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num):
    16	    """
    17	    Sample srcflux plugin: fitting spectrum for each individual obsid
    18	
    19	    This sample plugin uses sherpa to fit a spectral model, and
    20	    return an estimate of the flux w/ errors calculated with
    21	    the sample_flux routine.
    22	
    23	    We wrap the actual fitting code below in a try/except block
    24	    here so that we can always return nan's in the case of any error.
    25	    """
    26	
    27	    try:
    28	        return doit_specfit(infile, outroot, band, elo, ehi, src_num)
    29	    except Exception as wrong:
    30	        print(str(wrong))
    31	        print(f"Problem fitting {outroot} spectrum. Skipping it.")
    32	
    33	        return [ReturnValue("fitted_Nh", np.nan, "cm**22", "Fitted Absorption value"),
    34	                ReturnValue("kT", np.nan, "keV", "Plasma Temperature"),
    35	                ReturnValue("normalization", np.nan, "", "Spectrum Normalization"),
    36	                ReturnValue("reduced_statistic", np.nan, "", "Reduced Fit Statistic"),
    37	                ReturnValue("fit_statistic", np.nan, "", "Fit Statistic"),
    38	                ReturnValue("dof", np.nan, "", "Degrees of Freedom"),
    39	                ReturnValue("sample_flux", np.nan, "", "2-10 keV Sample Flux"),
    40	                ReturnValue("sample_flux_lo", np.nan, "", "2-10 Sample Flux Uncertainty Low"),
    41	                ReturnValue("sample_flux_hi", np.nan, "", "2-10 Sample Flux Uncertainty Low"),
    42	                ]
    43	
    44	
    45	def srcflux_merge_plugin(infile, outroot, band, elo, ehi, src_num):
    46	    """
    47	    Sample srcflux plugin: fitting spectrum for the combine spectra.
    48	
    49	    In this example it is the same plugin as the per-observation plugin,
    50	    but it does not have to be.
    51	    """
    52	
    53	    return srcflux_obsid_plugin(infile, outroot, band, elo, ehi, src_num)
    54	
    55	
    56	def doit_specfit(infile, outroot, band, elo, ehi, src_num):
    57	    'Perform the spectral fit: an absorbed APEEC model'
    58	    
    59	    from sherpa.utils.logging import SherpaVerbosity
    60	
    61	    # Find the spectrum file (per-obi vs. merged file names)
    62	    if os.path.exists(outroot+".pi"):
    63	        pi_file = outroot+".pi"
    64	    elif os.path.exists(outroot+"_src.pi"):
    65	        pi_file = outroot+"_src.pi"
    66	    else:
    67	        raise IOError(f"Unable to locate source spectrum for this source number: {src_num}")
    68	
    69	    # Suppress some of the sherpa chatter
    70	    with SherpaVerbosity('WARN'):
    71	        load_data(pi_file)
    72	        group_counts(10)
    73	        set_source(xsphabs.abs1*xsapec.thrm)
    74	        set_method("neldermead")
    75	        subtract()
    76	        counts = calc_data_sum(lo=0.5, hi=7.0)
    77	        if counts < 100:
    78	            raise ValueError(f"{pi_file} only has {counts} counts. Need 100 to do spectral fit")
    79	        notice(0.5, 8.0)
    80	        fit()
    81	        fit_info = get_fit_results()
    82	        fflux, cflux, vals = sample_flux(lo=2.0, hi=10.0)
    83	        f0, fhi, flo = cflux
    84	        plot_fit_delchi()
    85	        plt.savefig(outroot+".png")
    86	
    87	    return [ReturnValue("fitted_Nh", abs1.nH.val, "cm**-22", "Fitted Absorption value"),
    88	            ReturnValue("kT", thrm.kT.val, "keV", "Plasma Temperature"),
    89	            ReturnValue("normalization", thrm.norm.val, "", "Spectrum Normalization"),
    90	            ReturnValue("reduced_statistic", fit_info.rstat, "", "Reduced Fit Statistic"),
    91	            ReturnValue("fit_statistic", fit_info.statval, "", "Fit Statistic"),
    92	            ReturnValue("dof", fit_info.dof, "", "Degrees of Freedom"),
    93	            ReturnValue("sample_flux", f0, "", "2-10 keV Sample Flux"),
    94	            ReturnValue("sample_flux_lo", flo, "", "2-10 keV Sample Flux Uncertainty Low"),
    95	            ReturnValue("sample_flux_hi", fhi, "", "2-10 Sample Flux Uncertainty Low"),
    96	            ]

Below we go into a little more detail about the code.

• Lines 10 and 11: srcflux requires that the plugin return an object that has 4 attributes: name, value, units, and description. We have chosen to use a namedtuple.
• Line 13: Tell matplotlib to use the non-interactive agg backend for plotting.
• Line 15: this is the function definition line for the per-observation plugin with each of the required parameters.
• Lines 33 through 41: srcflux requires that the plugin return a list of the return values even when an Exception occurs. To make this easier, the routine that does the real work, doit_specfit, is wrapped in a try/except block so that any Exceptions raised during the fit are handled and the script returns the list of parameters with all the values set to NaN.
• Line 45: this is the function definition line for the plugin to run on the merged dataset.
• Line 53: We have written our code to handle both the per-observation and merged datasets so the plugin for the merged data just calls the plugin for the per-observation data.
• Line 56: The routine that does the actual spectral fitting on either the per-obsid or the merged products.
• Lines 62-67: Identifies the spectrum file to load based on the outroot; outroot is different for the per-observation and the merged datasets.
• Lines 59 and 70: this allow us to suppress some of the Sherpa chatter and is completely optional.
• Lines 71 through 85: this is the actual Sherpa fitting code. It performs some setup by loading the file, grouping the data, setting the source model, and subtracting the background. It also checks to see if there are fewer than 100 counts in which case it skips the fit. Then the fit is performed, the results are used to compute an estimate of the flux in the 2-10 keV range, and we create a plot of the fit with residuals. This could be expanded to compute confidence intervals on the fitted parameters, to do something iterative like freeze some parameter to do an initial fit and then thaw them to perform a final fit, use set_xsxset to set any of the internal XSPEC parameter, etc.
• Lines 87 to 96: This returns a list of the parameter values that we have chose to return. The list of the parameters returned in Lines 30-39 when an Exception occurs must match the same list of parameters being returned here. Users may chose to return an scalar values; this could be confidence intervals, execution time, etc.