Extract Spectrum and Response Files for Multiple Sources
CIAO 4.12 Science Threads
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.
Create a spectrum for a set of sources in two ways:
- a number of point sources present in the field, coadded
- treating the field of sources as a single, extended object
The specextract script automates the analysis for each case when combining data from separate observations.
- Analysis Guide: Extended Sources
Last Update: 8 Apr 2019 - Updated sherpa commands for matplotlib plotting.
- Get Started
- Option A: Coadding Multiple Spectra
- Option B: Extracting a Weighted Spectrum
- Improving the weights
- Analysis Caveats
- Parameter files:
Download the sample data: 3207 (ACIS-S, 3C 294)
unix% download_chandra_obsid 3207
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 psfecf=0.393 Running fluximage Version: 28 November 2018 Using CSC ACIS broad science energy band. Aspect solution pcadf131209709N003_asol1.fits found. Bad-pixel file acisf03207_repro_bpix1.fits found. Mask file acisf03207_000N003_msk1.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 Creating PSF map for obsid 3207 The following files were created: The clipped counts image is: ./broad_thresh.img The clipped exposure map is: ./broad_thresh.expmap The PSF map is: ./broad_thresh.psfmap The exposure-corrected image is: ./broad_flux.imgunix% cd fimg 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
Figure 1: The source regions
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 4.1 global color=blue dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1 physical Ellipse(3939.65,3732.18,5.28798,3.19063,15.4731) # Ellipse(3806.98,3754.33,7.08125,4.18562,26.0989) # Ellipse(4312.56,3827.05,3.9379,2.64072,50.7592) # Ellipse(4048.22,3970.57,2.44139,1.96386,20.4309) # Ellipse(4179.67,4078.33,2.33642,2.05695,9.03762) # Ellipse(4439.89,4091.26,6.46662,3.97055,92.5572) # Ellipse(4190.34,4140.68,2.40601,2.13286,46.3592) # 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 Ellipse src.reg | sed 's:Ellipse(.*):repro/acisf03207_repro_evt2.fits[sky=&]:;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=no
The final pset - combine=no - means the specextract will not combine_spectra to coadd the spectra and responses. Earlier versions of this thread ran with combine=yes however, that is only appropraite when combining data from different observations. Since this thread is combining data from the same observation, it must be combined seprately.
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: 30 November 2016 Using event file acisf03207_repro_evt2.fits[sky=ellipse(3939.65,3732.18,5.30179,3.19864,15.4853)] Aspect solution file pcadf131209709N003_asol1.fits found. Bad-pixel file acisf03207_repro_bpix1.fits found. Mask file acisf03207_000N003_msk1.fits found. Using event file acisf03207_repro_evt2.fits[sky=ellipse(3806.98,3754.33,7.08125,4.18562,26.0989)] Aspect solution file pcadf131209709N003_asol1.fits found. ... Setting bad pixel file [1 of 8] Converting source region to physical coordinates [1 of 8] Extracting src spectra [1 of 8] Creating src ARF [1 of 8] Calculating aperture correction for src ARF [1 of 8] Creating src RMF [1 of 8] Using mkacisrmf... Grouping src spectrum [1 of 8] Updating header of 3c294_src1.pi with RESPFILE and ANCRFILE keywords. Updating header of 3c294_src1_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:
|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|
These spectra can now be combined together using the combine_spectra script. First we make a stack of the ungrouped source spectra:
unix% /bin/ls 3c294_src*.pi | grep -v _grp > src.lis unix% cat src.lis 3c294_src1.pi 3c294_src2.pi 3c294_src3.pi 3c294_src4.pi 3c294_src5.pi 3c294_src6.pi 3c294_src7.pi 3c294_src8.pi
and then we combine them. The default is for combine_spectra to sum the exposure times from the spectra. However, since all these spectra are from the same observation, the exposure time in the combined output is unchanged. This is done by setting method=avg as shown here:
unix% pset combine_spectra email@example.com unix% pset combine_spectra out=3c294_combined unix% pset combine_spectra method=avg unix% combine_spectra mode=h clobber=yes Prepared to combine 8 spectra source PHA: 3c294_src1.pi ARF: 3c294_src1.corr.arf RMF: 3c294_src1.rmf source PHA: 3c294_src2.pi ARF: 3c294_src2.corr.arf RMF: 3c294_src2.rmf source PHA: 3c294_src3.pi ARF: 3c294_src3.corr.arf RMF: 3c294_src3.rmf source PHA: 3c294_src4.pi ARF: 3c294_src4.corr.arf RMF: 3c294_src4.rmf source PHA: 3c294_src5.pi ARF: 3c294_src5.corr.arf RMF: 3c294_src5.rmf source PHA: 3c294_src6.pi ARF: 3c294_src6.corr.arf RMF: 3c294_src6.rmf source PHA: 3c294_src7.pi ARF: 3c294_src7.corr.arf RMF: 3c294_src7.rmf source PHA: 3c294_src8.pi ARF: 3c294_src8.corr.arf RMF: 3c294_src8.rmf The following files were created: 3c294_combined_src.pi 3c294_combined_src.arf 3c294_combined_src.rmf
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
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:
Proceed to the Fitting section for instructions on how to analyse the data.
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 sherap In : %matplotlib Using matplotlib backend: TkAgg sherpa In : load_pha("multi/3c294_combined_src.pi") read ARF file multi/3c294_combined_src.arf read RMF file multi/3c294_combined_src.rmf sherpa In : load_pha(2, "field/combined.pi") read ARF file field/combined.warf read RMF file field/combined.wrmf sherpa In : notice(0,5, 7) sherpa In : group_counts(1, 20) sherpa In : group_counts(2, 20) sherpa In : plot_data(1) sherpa In : plot_data(2, overplot=True) sherpa In : last_curve = plt.gca().lines[-1] sherpa In : last_curve.set_mfc("none") sherpa In : last_curve.set_markersize(8) sherpa In : last_curve.set_mec("red") sherpa In : last_curve.set_marker("s") sherpa In : plt.xscale("log") sherpa In : plt.yscale("log")
Figure 2: Both spectra plotted in Sherpa
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.
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.
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.
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.
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.
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)
|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 188.8.131.52 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).|
|23 Dec 2014||Review for CIAO 4.7.|
|30 Jan 2017||Updated for CIAO 4.9. Replaced the combine=yes step in specextract with a separate run of combine_spectra.|
|07 Dec 2018||Updated dmmakereg output for ds9 v4.1 syntax used in CIAO 4.11. Updated fluximage to create the PSF map.|
|08 Apr 2019||Updated sherpa commands for matplotlib plotting.|