Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 8 Aug 2013

URL: http://cxc.harvard.edu/ciao/threads/wresp_multiple_sources/

Extract Spectrum and Response Files for Multiple Sources

CIAO 4.7 Science Threads


Overview

Synopsis:

A set of sources may be treated as multiple point sources that are coadded together for spectral analysis or as a single, multi-source region. The appropriate Response Matrix Files (RMFs) and Ancillary Response Files (ARFs) are created for each spectrum. In the latter case, the responses are weighted because the RMF and ARF vary with detector location: the RMF is defined as constant across FEF tiles, whose size depends on both the position on and focal-plane temperature of the ACIS chip, while the ARF varies with position in the detector plane.

Purpose:

Create a spectrum for a set of sources in two ways:

The specextract script automates the analysis for each case.

Related Links:

Last Update: 8 Aug 2013 - Updated for the contributed scripts 4.5.4 release: the "ancillary" files - e.g. aspect solution and bad-pixel files - can now be picked up from information in the event file (in many cases).


Contents


Get Started

Download the sample data: 3207 (ACIS-S, 3C 294)

unix% chandra_repro 3207 outdir=
unix% cd 3207

A region file created by wavdetect is used as well; you should therefore have completed the Running wavdetect thread (or copy the region file listed below). Here we show a brief run through of creating such a source list; note that the parameter choices made here are for example purposes only and should not be assumed to be scientifically valid:

unix% punlearn fluximage
unix% fluximage "repro/acisf03207_repro_evt2.fits[ccd_id=7]" fimg/ bin=1
Running fluximage
Version: 06 March 2013

Using CSC ACIS broad science energy band.
Aspect solution repro/pcadf131209709N003_asol1.fits found.
Bad-pixel file repro/acisf03207_repro_bpix1.fits found.
Mask file repro/acisf03207_000N003_msk1.fits found.
Pbk file repro/acisf131209384N003_pbk0.fits found.

The output images will have 1312 by 1312 pixels, pixel size of 0.492 arcsec,
    and cover x=3361.5:4673.5:1,y=3153.5:4465.5:1.

Running tasks in parallel with 4 processors.
Creating aspect histogram for obsid 3207
Creating instrument map for obsid 3207
Creating exposure map for obsid 3207
Thresholding data for obsid 3207
Exposure-correcting image for obsid 3207

The following files were created:

 The clipped counts image is:
     fimg/broad_thresh.img

 The clipped exposure map is:
     fimg/broad_thresh.expmap

 The exposure-corrected image is:
     fimg/broad_flux.img

unix% cd fimg
unix% punlearn mkpsfmap
unix% mkpsfmap broad_thresh.img broad_thresh.psfmap 1.49 ecf=0.393
unix% punlearn wavdetect
unix% pset wavdetect infile=broad_thresh.img
unix% pset wavdetect outfile=wav.fits
unix% pset wavdetect scellfile=wav.scell
unix% pset wavdetect imagefile=wav.img
unix% pset wavdetect defnbkgfile=wav.nbkg
unix% pset wavdetect scales="2 4 8 16 32 64"
unix% pset wavdetect psffile=broad_thresh.psfmap
unix% pset wavdetect expfile=broad_thresh.expmap
unix% pset wavdetect regfile=wav.reg
unix% wavdetect mode=h

The source list has been edited to remove the three sources associated with the cluster (Figure 1) and to require a minimum of 100 net counts in the source (this latter constraint is just to reduce the number of sources to a manageable level for this example).

unix% dmlist wav.fits counts
44
unix% dmlist "wav.fits[exclude pos=circle(4070,3945,20)]" counts
41
unix% dmcopy "wav.fits[exclude pos=circle(4070,3945,20)]" no-clus.fits
unix% dmstat "no-clus.fits[cols net_counts]"
NET_COUNTS[count]
    min:	5.5470795631 	      @:	40
    max:	497.15258789 	      @:	28
   mean:	57.384727618
  sigma:	85.158878925
    sum:	2352.7738323
   good:	41
   null:	0

unix% dmlist "no-clus.fits[net_counts=100:]" counts
8
unix% dmcopy "no-clus.fits[net_counts=100:]" ../src.fits
[Thumbnail image: The regions are spread over ACIS-S3.]

[Version: full-size]

[Print media version: The regions are spread over ACIS-S3.]

Figure 1: The source regions

unix% ds9 -region width 2 broad_thresh.img -smooth -scale log -cmap b \
     -zoom 0.5 -pan to 4020 3990 physical -region ../src.fits &

We can convert the FITS format file to ASCII format using dmmakereg:

unix% cd ..
unix% punlearn dmmakereg
unix% dmmakereg region="region(src.fits)" out=src.reg kernel=ascii
unix% cat src.reg
# Region file format: DS9 version 3.0
global color=blue font="helvetica 10 normal" select=1 edit=1 move=1 delete=1 include=1 fixed=0
physical;Ellipse(3939.65,3732.18,5.28798,3.19063,15.4731) #
physical;Ellipse(3806.98,3754.33,7.08125,4.18562,26.0989) #
physical;Ellipse(4312.56,3827.05,3.9379,2.64072,50.7592) #
physical;Ellipse(4048.22,3970.57,2.44139,1.96386,20.4309) #
physical;Ellipse(4179.67,4078.33,2.33642,2.05695,9.03762) #
physical;Ellipse(4439.89,4091.26,6.46662,3.97055,92.5572) #
physical;Ellipse(4190.34,4140.68,2.40601,2.13286,46.3592) #
physical;Ellipse(3687.97,4356.76,7.13883,4.4737,149.623) #

Option A: Coadding Multiple Spectra

A pointlike spectrum, an aperture-corrected ARF, and an RMF are extracted for each source in the list. The spectra and responses are then coadded for the spectral analysis. For a step-by-step explanation of the analysis, read the Extract Spectrum and Response Files for a Pointlike Source thread and the Coadding Spectra and Responses thread.

For this approach we need a stack which applies each source region from src.reg to the event file, which can be obtained using manual editing or a command like:

unix% grep physical src.reg | sed -e "s:^physical;:repro/acisf03207_repro_evt2.fits[sky=:" -e "s/Ellipse/ellipse/" -e "s/ #/]/" > multi.lis
unix% cat multi.lis
repro/acisf03207_repro_evt2.fits[sky=ellipse(3939.65,3732.18,5.28798,3.19063,15.4731)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(3806.98,3754.33,7.08125,4.18562,26.0989)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(4312.56,3827.05,3.9379,2.64072,50.7592)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(4048.22,3970.57,2.44139,1.96386,20.4309)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(4179.67,4078.33,2.33642,2.05695,9.03762)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(4439.89,4091.26,6.46662,3.97055,92.5572)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(4190.34,4140.68,2.40601,2.13286,46.3592)]
repro/acisf03207_repro_evt2.fits[sky=ellipse(3687.97,4356.76,7.13883,4.4737,149.623)]

We set an output root, multi/3c294, and allow specextract to number the spectra and responses (multi/3c294_src1, multi/3c294_src2, etc):

unix% mkdir multi
unix% cd multi
unix% punlearn specextract
unix% pset specextract infile=@../multi.lis
unix% pset specextract outroot=3c294
unix% pset specextract weight=no
unix% pset specextract correctpsf=yes
unix% pset specextract combine=yes

The weight parameter is set to "no" to tell the script to make a pointlike ARF and the correctpsf parameter is set to "yes" so that arfcorr will be run to correct the ARF. The final pset - combine=yes - means the specextract will call combine_spectra to coadd the spectra and responses.

Note that the script picks up the location of the ancillary files - aspect solution, bad-pixel, mask file and parameter block - from keywords in the event file, even though they are not in the current working directory:

unix% specextract mode=h

Running: specextract
  Version:    12 July 2013
Using event file ../repro/acisf03207_repro_evt2.fits[sky=ellipse(3939.65,3732.18,5.28798,3.19063,15.4731)]
Aspect solution(s) ../repro/pcadf131209709N003_asol1.fits found.
Bad-pixel file ../repro/acisf03207_repro_bpix1.fits found.
ACIS parameter block file ../repro/acisf131209384N003_pbk0.fits found.
Mask file ../repro/acisf03207_000N003_msk1.fits found.
Using event file ../repro/acisf03207_repro_evt2.fits[sky=ellipse(3806.98,3754.33,7.08125,4.18562,26.0989)]
Aspect solution(s) ../repro/pcadf131209709N003_asol1.fits found.
Bad-pixel file ../repro/acisf03207_repro_bpix1.fits found.
ACIS parameter block file ../repro/acisf131209384N003_pbk0.fits found.
Mask file ../repro/acisf03207_000N003_msk1.fits found.
...

Setting bad pixel file for item 1 of 8 in input list


Converting source region to physical coordinates for item 1 of 8 in input list.


Extracting src spectra for item 1 of 8 in input list


Creating src ARF for item 1 of 8 in input list


Creating src RMF for item 1 of 8 in input list


Using mkacisrmf...

Grouping src spectrum for item 1 of 8 in input list

Updating header of 3c294_src1.pi with RESPFILE and ANCRFILE keywords.
Updating header of 3c294_src1_grp.pi with RESPFILE and ANCRFILE keywords.

Setting bad pixel file for item 2 of 8 in input list

...

Grouping src spectrum for item 8 of 8 in input list

Updating header of 3c294_src8.pi with RESPFILE and ANCRFILE keywords.
Updating header of 3c294_src8_grp.pi with RESPFILE and ANCRFILE keywords.

The contents of the parameter file may be checked with plist specextract.

The important output files for each source are:

Filename Description
3c294_src<n>.pi Ungrouped spectrum
3c294_src<n>_grp.pi Grouped spectrum
3c294_src<n>.corr.arf ARF after being corrected by arfcorr
3c294_src<n>.arf ARF generated by mkarf
3c294_src<n>.rmf RMF generated by mkacisrmf

and the coadded spectrum and response files to use in the fitting are:

3c294_combined_src.arf
3c294_combined_src.pi
3c294_combined_src.rmf

Note that combine_spectra ignored any grouping flags in the input spectra, so the output spectrum is ungrouped.

Proceed to the Fitting section for instructions on how to analyse the data.


Option B: Extracting a Weighted Spectrum

A weighted spectrum and response files are extracted from the set of sources, essentially treating them as an extended source. For a step-by-step explanation of the analysis, read the Extract Spectrum and Response Files for an Extended Source thread.

unix% punlearn specextract
unix% pset specextract infile="repro/acisf03207_repro_evt2.fits[sky=region(src.fits)]"
unix% pset specextract outroot=field/combined

Note that unlike the point-source approach, we use the default settings for the weight, correctpsf, and combine parameters, namely yes, no, and no respectively.

unix% specextract mode=h

Running: specextract
  Version:    12 July 2013
Using event file repro/acisf03207_repro_evt2.fits[sky=region(src.fits)]
Aspect solution(s) repro/pcadf131209709N003_asol1.fits found.
Bad-pixel file repro/acisf03207_repro_bpix1.fits found.
ACIS parameter block file repro/acisf131209384N003_pbk0.fits found.
Mask file repro/acisf03207_000N003_msk1.fits found.

Setting bad pixel file for item 1 of 1 in input list


Extracting src spectra for item 1 of 1 in input list


Creating src ARF for item 1 of 1 in input list


Creating src RMF for item 1 of 1 in input list


Using mkacisrmf...

Grouping src spectrum for item 1 of 1 in input list

Updating header of field/combined.pi with RESPFILE and ANCRFILE keywords.
Updating header of field/combined_grp.pi with RESPFILE and ANCRFILE keywords.

The contents of the parameter file may be checked with plist specextract.

The output files to use in fitting are:

Filename Description
field/combined.pi Ungrouped spectrum
field/combined_grp.pi Grouped spectrum
field/combined.warf ARF
field/combined.wrmf RMF

Proceed to the Fitting section for instructions on how to analyse the data.


Fitting

Either source spectrum can now be analyzed - using the related ARF and RMF - in Sherpa, as described in the Introduction to Fitting PHA Spectra thread.

In Figure 2 we plot the two data sets using Sherpa:

unix% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.4 Sherpa version 1 Friday, December 2, 2011

sherpa-1> load_pha("multi/3c294_combined_src.pi")
read ARF file multi/3c294_combined_src.arf
read RMF file multi/3c294_combined_src.rmf
sherpa-2> load_pha(2, "field/combined.pi")
read ARF file field/combined.warf
read RMF file field/combined.wrmf
sherpa-3> notice(0,5, 7)
sherpa-4> group_counts(1, 20)
sherpa-5> group_counts(2, 20)
sherpa-6> plot_data(1)
sherpa-7> set_curve(["*.color", "red"])
sherpa-8> plot_data(2, overplot=True)
sherpa-8> log_scale()
[Thumbnail image: The two spectra look to have the same shape but a different scaling factor.]

[Version: full-size]

[Print media version: The two spectra look to have the same shape but a different scaling factor.]

Figure 2: Both spectra plotted in Sherpa

The red curve shows multi/3c294_combined_src.pi and the black curve field/combined.pi, both having been grouped to 20 counts per bin.


Improving the weights

In these examples, we used the detected counts to weight the responses, when what what we really should be using is the incident flux - i.e. the source spectrum before it passes through the telescope/detector system. If we assume an incident spectrum, we can correct the detected counts to give the incident flux, and use that for the weights. Since we do not know the source spectrum, we can use an iterative scheme - where the best-fit spectrum of the source from the previous iteration is used to correct the WMAP for the current iteration, until the results converge. An alternative scheme is to continue to use the detected counts - i.e. make no correction for the telescope/detector response - but only use a restricted energy range where the response + model does not vary much; an example being the soft energy band for clusters of galaxies.

The spectrumfile parameter of mkwarf accepts an ASCII file containing a model of the source spectrum, which may be created by running the Calculating Spectral Weights thread. Note that the energy range of the spectral model should match that used to create the WMAP. At this time specextract does not accept a spectrumfile.

Note that, unlike mkinstmap, the input spectrum does not need to be normalized to unit flux.


Analysis Caveats

The algorithm used to calculate the weighted responses requires that there are no spatial variations in the source spectrum.

Users should be cautious about analyzing the data for sources near the edges of the ACIS CCDs.

  1. For X-rays passing through the mirrors, the very bottom of each CCD is obscured by the frame store. As a result, some of the events in rows with CHIPY <= 8 are not detected. (The set of rows affected varies from CCD to CCD.) Since the CIAO tools do not compensate for this effect, the ARFs and exposure maps for sources in these regions may be inaccurate.

  2. For sources within about thirty-two pixels of any edge of a CCD, the source may be dithered off the CCD during part of an observation. The aspect histogram, which is used to create ARFs and exposure maps, is designed to compensate for this effect.

  3. An ARF calculated at the edge of a chip will not be accurate since the response tools for spectral extraction (specifically the ARF) assume that 100% of the PSF is enclosed - i.e. on the chip - all the time, which may not be the case. The amount of error introduced depends on how close the source is to the edge, the morphology of the source, and the characteristics of the PSF, which depends on the source spectrum.

  4. A contaminant has accumulated on the optical-blocking filters of the ACIS detectors, as described in the ACIS QE Contamination why topic. Since there is a gradient in the temperature across the filters (the edges are colder), there is a gradient in the amount of material on the filters. (The contaminant is thicker at the edges.) Within about 100 pixels of the outer edges of the ACIS-I and ACIS-S arrays, the gradient is relatively steep. Therefore, the effective low-energy (' 1 keV) detection efficiency may vary within the dither pattern in this region. The ARF and instrument map tools are designed to read a calibration file which describes this spatial dependence.



Parameters for /home/username/cxcds_param/specextract.par


        infile = @../multi.lis    Source event file(s)
       outroot = 3c294            Output directory path + root name for output files
      (bkgfile = )                Background event file(s)
          (asp = )                Source aspect solution or histogram file(s)
      (pbkfile = )                Parameter block file (input to mkwarf)
      (mskfile = )                Maskfile (input to mkwarf)
      (rmffile = CALDB)           rmffile input for CALDB
   (badpixfile = )                Bad pixel file for the observation
       (dafile = CALDB)           Dead area file (input to mkwarf)
      (bkgresp = yes)             Create background ARF and RMF?
       (weight = no)              Should response files be weighted?
   (correctpsf = yes)             Apply point source aperture correction to ARF?
      (combine = yes)             Combine ungrouped output spectra and responses?
    (grouptype = NUM_CTS)         Spectrum grouping type (same as grouptype in dmgroup)
      (binspec = 15)              Spectrum grouping specification (NONE,1:1024:10,etc)
(bkg_grouptype = NONE)            Background spectrum grouping type (NONE, BIN, SNR, NUM_BINS, NUM_CTS, or ADAPTIVE)
  (bkg_binspec = )                Background spectrum grouping specification (NONE,10,etc)
        (ptype = PI)              PI or PHA
       (energy = 0.3:11.0:0.01)   Energy grid
      (channel = 1:1024:1)        RMF binning attributes
  (energy_wmap = 300:2000)        Energy range for (dmextract) WMAP input to mkacisrmf
      (binwmap = tdet=8)          Binning factor for (dmextract) WMAP input to mkacisrmf
   (binarfwmap = 1)               Binning factor for (sky2tdet) WMAP input to mkwarf
      (clobber = no)              OK to overwrite existing output file?
      (verbose = 1)               Debug Level(0-5)
         (mode = ql)
    


Parameters for /home/username/cxcds_param/specextract.par


        infile = repro/acisf03207_repro_evt2.fits[sky=region(src.fits)] Source event file(s)
       outroot = field/combined   Output directory path + root name for output files
      (bkgfile = )                Background event file(s)
          (asp = )                Source aspect solution or histogram file(s)
      (pbkfile = )                Parameter block file (input to mkwarf)
      (mskfile = )                Maskfile (input to mkwarf)
      (rmffile = CALDB)           rmffile input for CALDB
   (badpixfile = )                Bad pixel file for the observation
       (dafile = CALDB)           Dead area file (input to mkwarf)
      (bkgresp = yes)             Create background ARF and RMF?
       (weight = yes)             Should response files be weighted?
   (correctpsf = no)              Apply point source aperture correction to ARF?
      (combine = no)              Combine ungrouped output spectra and responses?
    (grouptype = NUM_CTS)         Spectrum grouping type (same as grouptype in dmgroup)
      (binspec = 15)              Spectrum grouping specification (NONE,1:1024:10,etc)
(bkg_grouptype = NONE)            Background spectrum grouping type (NONE, BIN, SNR, NUM_BINS, NUM_CTS, or ADAPTIVE)
  (bkg_binspec = )                Background spectrum grouping specification (NONE,10,etc)
        (ptype = PI)              PI or PHA
       (energy = 0.3:11.0:0.01)   Energy grid
      (channel = 1:1024:1)        RMF binning attributes
  (energy_wmap = 300:2000)        Energy range for (dmextract) WMAP input to mkacisrmf
      (binwmap = tdet=8)          Binning factor for (dmextract) WMAP input to mkacisrmf
   (binarfwmap = 1)               Binning factor for (sky2tdet) WMAP input to mkwarf
      (clobber = no)              OK to overwrite existing output file?
      (verbose = 1)               Debug Level(0-5)
         (mode = ql)
    

History

04 Jan 2005 updated for CIAO 3.2: added information about mkacisrmf, see Creating the RMF: mkrmf vs mkacisrmf section (deprecated)
20 Dec 2005 updated for CIAO 3.3: default value of dmextract error and bkgerror parameters is "gaussian"; updated screen output accordingly
01 Feb 2006 added information about specextract thread
01 Dec 2006 updated for CIAO 3.4: CIAO and ChIPS versions; set mkrmf verbosity > 0 for screen output; parameter file updates for mkwarf
02 Feb 2007 updated for CALDB 3.3.0.1 patch
06 Mar 2007 added ACIS dead area correction section and example of setting the pbkfile and dafile parameters
24 Jan 2008 updated for CIAO 4.0: rewritten with ObsID 3207; turn off the ACIS dead area correction in the mkwarf step (application of the dead area correction is on by default); show_wgt.sl not included in this release; links point to Sherpa Beta website ; removed outdated calibration updates
04 Feb 2009 updated for CIAO 4.1: "ARDLIB warning ... Assuming the first "interesting" extension." no longer printed; updated path for CALDB 4; input data must have a CTI_APP keyword
19 Feb 2009 added a section on Fitting
12 Jan 2010 updated for CIAO 4.2: calibration update - the ACIS QE contamination model has been upgraded to vN0005.
09 Mar 2010 The ACIS detector is calibrated over the range 0.224004 - 12 keV; choosing values outside this range may result in errors from mkwarf.
15 Dec 2010 updated for CIAO 4.3: new ACIS contamination calibration file
01 Mar 2011 CALDB 4.4.2 release: fix to the header of the ACIS QE contamination file. Prior to this release, CIAO would fail when trying to look up the contamination model correction for chips ACIS-8 (S4) and ACIS-9 (S5).
26 Apr 2011 install version 2 of the tools package for CIAO 4.3 to fix the mkrmf bug.
20 Jul 2011 required software updates are listed in Synopsis
05 Mar 2012 revised for CIAO 4.4 (previous title was “Weighting ARFs and RMFs: multiple sources”): thread rewritten to show two cases - stack of sources and single, multi-source region; analysis is done with the specextract script.
03 Dec 2012 Review for CIAO 4.5; updated file versions
08 Aug 2013 Updated for the contributed scripts 4.5.4 release: the "ancillary" files - e.g. aspect solution and bad-pixel files - can now be picked up from information in the event file (in many cases).


Last modified: 8 Aug 2013
Smithsonian Institute Smithsonian Institute

The Chandra X-Ray Center (CXC) is operated for NASA by the Smithsonian Astrophysical Observatory. 60 Garden Street, Cambridge, MA 02138 USA.   Email:   cxcweb@head.cfa.harvard.edu Smithsonian Institution, Copyright © 1998-2014. All rights reserved.