Color-color diagram

In X-ray astronomy we usually define "color" as the hardness ratio.

The hardness ratio is the ratio of the difference in the number of counts in two independent energy bands divided by the total number of counts. The total number of counts may be the sum of the two energy bands, or may be the total over a wider energy range that full encloses the bands being used.

Thus defined, the hardness ratio will be between -1 and 1. Values closer to -1 are "soft", low energy band dominates. Values closer to 1 are "hard", high energy band dominates.

A color-color plot then shows the hardness ratio in 1 pair of energy bands along the X-axis, and a different pair of energy bands along the Y-axis.

This collection of classes simulates a spectrum with the user specified spectral model and instrument response. Two of the model parameters are varied over an input grid and the color in the different energy bands is computed. This grid is then plotted in the color-color diagram.

The idea is that for a source, if you compute the colors and assume a spectral model shape, you can get an estimate of the model parameters by looking at the color-color diagram.

These routines are dependent on Sherpa for the models. Plotting is done using matplotlib.

Examples

Command line tool

The plot shown above was created using the following command line

$ color_color acis.arf clr.fits mode=h showplot=yes outplot=clr.png \
    model="xsphabs.abs1*xspowerlaw.pwrlaw" \
    param1=pwrlaw.PhoIndex \
    grid1=1,2,3,4 \
    param2=abs1.nH \
    grid2=0.01,0.1,0.2,0.5,1,10 \
    soft=csc medium=csc hard=csc \
    clobber=yes

The values are also stored in the output file clr.fits

$ dmlist clr.fits cols
 
--------------------------------------------------------------------------------
Columns for Table Block TABLE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   PhoIndex                          Real8          -Inf:+Inf            
   2   nH                                Real8          -Inf:+Inf            
   3   S_COUNTS                          Real8          -Inf:+Inf            
   4   M_COUNTS                          Real8          -Inf:+Inf            
   5   H_COUNTS                          Real8          -Inf:+Inf            
   6   HARD_HM                           Real8          -Inf:+Inf            
   7   HARD_MS                           Real8          -Inf:+Inf            

Python module

from color_color import *

import sherpa.astro.ui as ui

#
# Define our energy bands
#
soft = EnergyBand( 0.5, 1.2, 'S')
medium = EnergyBand(1.2, 2.0, 'M')
hard = EnergyBand(2.0, 7.0, 'H')

#
# Define model
#
mymodel = ui.xswabs.abs1 * ui.xspowerlaw.pwrlaw
arffile = "acissD2006-10-26pimmsN0009.fits"

#
# First model parameter axis
#
pho_grid = [ 1., 2., 3., 4. ]
photon_index = ModelParameter( pwrlaw.PhoIndex, pho_grid)

#
# Second model parameter axis
# 
sg = [ 1.e20, 1.e21, 2.e21, 5.e21, 1.e22, 1e23] 
nh_grid = [x/1e22 for x in sg ]
absorption = ModelParameter( abs1.nH, nh_grid)

#
# Get to work.  
#
hm_ms = ColorColor( mymodel, arffile )
hm_ms_results = hm_ms( photon_index, absorption, soft, medium, hard)

#
# Setup plot styles
#
photon_index.set_curve_style(marker="", color="black" )
absorption.set_curve_style(marker="", linestyle="--", color="forestgreen") 
absorption.set_label_style(color="forestgreen")

hm_ms_results.plot("color_color.png")

vals = hm_ms_results.get_results()

The output looks like

Now suppose there was a source with (H-M)=-0.25 and (M-S)=-0.2. If the assumption about the spectral model is correct, then it would be expected to have a photon index between 2 and 3 and absorption between 0.2 and 0.5 e10^22.

Total energy band

If the total energy band is None then the routines will use the sum of the individual bands independently for each axis. In this case the color values will fill the entire -1:1 parameter space.

If the total energy band is specified, then it may be the case that only part of the parameter space is filled.

Background

Background is not treated specifically with these routines, but you can include the background term in the model expression, something like

mymodel = ui.xswabs.abs1 * ui.xspowerlaw.pwrlaw + ui.const1d.bkg
bkg.c0 = 0.01

to see how it affects things. The background model component can be as complicated as needed.