Last modified: 15 Feb 2022

URL: https://cxc.cfa.harvard.edu/ciao/threads/tgextract2/

Applying Customized Background Regions to LETG/HRC-S Observations

CIAO 4.16 Science Threads


Overview

Synopsis:

This thread demonstrates the use of tgextract2 to apply customized background regions to an LETG/HRC-S observation of a field with two, close, bright sources. With LETGS, it is important to include backgrounds in the analysis since the HRC-S instrumental background is significant.

Each source in the sample data (ObsID 29) falls into one of the default background regions for the other. The sources themselves are far enough apart that their spectra can be extracted cleanly. Creating a custom region will avoid background contamination by either of the sources.

Some detailed and interactive analysis is required, involving some scientific judgement by the analyst, to determine the background region to apply.

Related Links:

Last Update: 15 Feb 2022 - Review for CIAO 4.14. No updates needed.


Contents


Introduction to tgextract2

tgextract2 is an alternative to tgextract, the CIAO tool which filters and bins spectra. tgextract2 provides support for variable-shaped background regions which do not need to have a constant width ratio with the source region. tgextract2 computes a vector "backscale"—the scaling factor which relates the background count rate at any wavelegth to the source extraction region width at that wavelength. tgextract2 also sums the background extractions from each side of the spectrum, whereas tgextract writes two counts arrays, one from the "_UP" and one from the "_DOWN" sides of the spectrum.

Since tgextract2 is not backward-compatible with tgextract, it is not used in standard processing. It is instead recommended for special cases where some detailed and interactive analysis is required, involving some scientific judgement by the analyst.


Getting Started

Download the sample data: 29 (LETG/HRC-S, Alpha Cen)

The standard HRC grating region file, found in CalDB, is also used:

ln -s $CALDB/data/chandra/hrc/tgmask2/letgD1999-07-22regN0002.fits .

Open the event file in DS9:

unix% ds9 hrcf00029N006_evt2.fits &

Next, we set the height and width to something convenient for grating spectra, and adjust the colorscale and colormap in the DS9 menu.

  • Frame → Frame Parameters → Set Display Size → X=1024, Y=400
  • Scale → Log
  • Color → BB

Figure 1 shows the resulting image.

Figure 1: Source Image

[Thumbnail image: CCD image of source after color parameters have been changed]

[Version: full-size]

[Print media version: CCD image of source after color parameters have been changed]

Figure 1: Source Image

Source image with BB colormap and log-scale selected

Next, we set the binning parameters in the DS9 menu Bin → Binning Parameters (Figure 2):

Figure 2: Grating Coordinates Binning Parameters

[DS9 menu settings for binning the source image and converting it to grating dispersion coordinates to provide a grating spectrum.]
[Print media version: DS9 menu settings for binning the source image and converting it to grating dispersion coordinates to provide a grating spectrum.]

Figure 2: Grating Coordinates Binning Parameters

Binning parameter settings applied in DS9 by going to the menu and choosing Bin → Binning Parameters.

The binning was iteratively chosen such that the field fits the screen and it has enough counts/bin to give good statistics. Fiddling with parameters, we start off with knowing that the tg_r coordinate range is -1 to 1 degrees in about 10000 bins, or 0.0002 deg/bin. Binning by a factor of 10, we have a tg_d range of about -0.002 to 0.002 deg/bin.

After applying the binning parameters, the field image in grating dispersion coordinates shows two spectra (Figure 3):

Figure 3: Grating Spectrum of the Two Field Sources

[Thumbnail image: The grating spectrum of the two field sources after binning parameters are applied.]

[Version: full-size]

[Print media version: The grating spectrum of the two field sources after binning parameters are applied.]

Figure 3: Grating Spectrum of the Two Field Sources

This is the image in grating dispersion coordinates with tg_r (in-dispersion) on the horizontal axis, and tg_d (cross-dispersion) on the vertical axis (both in degrees). We can clearly see the two sources zero-orders near physical coordinates (0,0) and their dispersed and line-rich spectra, with negative-orders to the left and positive-orders to the right. The vertical dark regions are gaps between the HRC-S microchannel plates. Overplotted is the standard LETG extraction region, where, in this case, the lower background region includes the spectrum of the second source and therefore should not be used.


Customizing the Extraction Regions

The standard HRC/LETG extraction region file must be edited for the current field. This is most easily done in DS9 as long as the information about which region represents source and up and down background is preserved. The tgmask2reg and reg2tgmask scripts allow the use of DS9 to define custom extraction regions, and use the DS9 region "tagging" feature to encode the grating information along with the extraction region coordinates.

The first step in the editing process requires running tgmask2reg to convert the CALDB tgmask file into a DS9 format. The additional (non-coordinate) columns in the CALDB files are added as "tags" so that later tools know each regions' purpose.

unix% tgmask2reg infile=CALDB outfile=tg.reg
unix% more tg.reg 
# Region file format: DS9 version 4.1
# Filename:
global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=0 include=1 source=1
physical
polygon(-0.0057804,0.000531,-0.29867,0.000531,-1.1549,0.002114,-1.1549,-0.002114,-0.29867,-0.000531,-0.0057804,-0.000531,-0.0057804,0.000531) # tag={row_idx=1} tag={rowid=SOURCE} tag={tg_part=3} tag={tg_srcid=1} tag={tg_m=-1} tag={backscal=1.0}
polygon(-0.0057804,0.00731,-0.29867,0.00731,-1.1549,0.02314,-1.1549,0.002,-0.29867,0.002,-0.0057804,0.002,-0.0057804,0.00731) # tag={row_idx=2} tag={rowid=BACKGROUND_UP} tag={tg_part=3} tag={tg_srcid=1} tag={tg_m=-1} tag={backscal=5.0}
polygon(-0.0057804,-0.00731,-0.29867,-0.00731,-1.1549,-0.02314,-1.1549,-0.002,-0.29867,-0.002,-0.0057804,-0.002,-0.0057804,-0.00731) # tag={row_idx=3} tag={rowid=BACKGROUND_DOWN} tag={tg_part=3} tag={tg_srcid=1} tag={tg_m=-1} tag={backscal=5.0}
polygon(0.0057804,0.000531,0.29867,0.000531,1.1549,0.002114,1.1549,-0.002114,0.29867,-0.000531,0.0057804,-0.000531,0.0057804,0.000531) # tag={row_idx=1} tag={rowid=SOURCE} tag={tg_part=3} tag={tg_srcid=1} tag={tg_m=1} tag={backscal=1.0}
polygon(0.0057804,0.00731,0.29867,0.00731,1.1549,0.02314,1.1549,0.002,0.29867,0.002,0.0057804,0.002,0.0057804,0.00731) # tag={row_idx=2} tag={rowid=BACKGROUND_UP} tag={tg_part=3} tag={tg_srcid=1} tag={tg_m=1} tag={backscal=5.0}
polygon(0.0057804,-0.00731,0.29867,-0.00731,1.1549,-0.02314,1.1549,-0.002,0.29867,-0.002,0.0057804,-0.002,0.0057804,-0.00731) # tag={row_idx=3} tag={rowid=BACKGROUND_DOWN} tag={tg_part=3} tag={tg_srcid=1} tag={tg_m=1} tag={backscal=5.0}

Once the region is loaded into DS9 (Region → Load Regions → tg.reg), it can be customized interactivelty as needed either by moving points or adding new ones both for the source and the background areas. Note that polygon regions cannot added or deleted. They can only be modified. This ensures that the needed meta-data is preserved.

Note that in DS9 regions can be inspected to see the extra information associated with each region by going to Region → Groups . This will provide a list of the different tags in use. Selecting one of the tags will show and select the corresponding region(s) in the main ds9 window.

The final customized set of regions must be saved in DS9 format in physical coordinates (Region → Save Regions → ds9.reg) (see Figure 4).

Figure 4: Customized Extraction Regions

[Thumbnail image: Customized LETG extraction regions]

[Version: full-size]

[Print media version: Customized LETG extraction regions]

Figure 4: Customized Extraction Regions

This is the modified extraction region for source and background: the lower background regions have been moved well away from the second source spectrum while the source region has been left unaltered.

Finally the reg2tgmask tool converts the modified region ("ds9.reg") back into a FITS format ("tg_mask.fits") compatible with tgextract2.

unix% reg2tgmask infile=ds9.reg srcfile=CALDB out=tg_mask.fits 


Running tgextract2

With the background regions appropriately defined, we are now ready to run tgextract2. The opt parameter is set to "pha1" so that the tool creates separate files for the +1 and -1 spectra (Type 1 PHA file), rather than both spectra be written as rows in the same file (Type 2 PHA file); see the grating spectra data products guide for more information these formats. The files will also contain the background extractions and the "areas" (BACKSCAL) of the source and background regions for each spectral bin.

The outfile parameter value (a_) is used as the root name of the output files. The files will be named as "root<srcid><part><order>_pha2.fits", one for each source ID, part, and order combination.

unix% punlearn tgextract2
unix% pset tgextract2 infile=hrcf00029N006_evt2.fits
unix% pset tgextract2 outfile=a_ 
unix% pset tgextract2 region_file=tg_mask.fits

unix% pset tgextract2 opt=pha1

unix% tgextract2 
Input event file (output event file from L1.5 processing) (hrcf00029N006_evt2.fits):
Output pha file (for pha2, enter full file name or '.'; for pha1, enter rootname.) (a_):
Source ID's to process: 'all', comma list, @file (1):
Grating parts to process: HETG, HEG, MEG, LETG, header_value (HETG|HEG|MEG|LETG|header_value) (header_value):
Grating diffraction orders to process: 'default', comma list, range list, @file (default):

This produces two new files, a_01Lm1_pha2.fits and a_01Lp1_pha2.fits, which are -1 order and +1 order LEG spectra, respectively.

You can check the parameter file that was used with plist tgextract2.

After we run tgextract2, we will copy the background counts to a separate file, one for each order. This is because various analysis softwares typically need to have the background spectrum formatted the same way as the source (i.e. having the same column names).

unix% dmcopy 'a_01Lm1_pha2.fits[cols CHANNEL,BIN_LO,BIN_HI,COUNTS=bg_counts,BACKSCAL=bg_area]' a_01Lm1_pha2_bg.fits
unix% dmcopy 'a_01Lp1_pha2.fits[cols CHANNEL,BIN_LO,BIN_HI,COUNTS=bg_counts,BACKSCAL=bg_area]' a_01Lp1_pha2_bg.fits
unix% ls a_*
a_01Lm1_pha2.fits     a_01Lm1_pha2_bg.fits  a_01Lp1_pha2.fits     a_01Lp1_pha2_bg.fits

Inspecting the Counts and Background Spectra with Sherpa

We will now inspect the counts and background spectra. Using Sherpa, we will create a strip-chart displaying the plus- and minus-order counts spectra, with background over-plotted. We begin by loading the data and background with the minus-order spectrum associated with dataset ID=1 and positive-order spectrum associated with dataset ID=2.

unix% sherpa

sherpa In []: load_pha(1,"a_01Lm1_pha2.fits")
sherpa In []: load_pha(2,"a_01Lp1_pha2.fits")
sherpa In []: load_bkg(1,"a_01Lm1_pha2_bg.fits")
sherpa In []: load_bkg(2,"a_01Lp1_pha2_bg.fits")

Next we set Sherpa plotting preferences. The key parameter is set_analysis, which plots the spectrum as counts vs. wavelength.

sherpa In []: set_analysis("wave",type="counts")

Loading a two plot strip chart, we first add the minus-order counts spectrum to the upper plot:

sherpa In []: %matplotlib
sherpa In []: import matplotlib.pylab as plt
sherpa In []: plt.subplot(2,1,1)
sherpa In []: plt.subplots_adjust(hspace=0.0 )
sherpa In []: plt.plot(get_data_plot(1).x,get_data_plot(1).y,marker="None", color="black")

Next, we over-plot the minus-order count spectrum with the minus-order background in red:

sherpa In []: bs1=get_bkg_scale(1)
sherpa In []: ct1=get_counts(1,bkg_id=1)

sherpa In []: plt.plot(get_data_plot(1).x,(ct1*bs1),marker="None", color="red")

and add some labels to plot

sherpa In []: plt.ylabel("Counts") 
sherpa In []: plt.figtext(0.93,0.75,"Negative-Order",rotation=270)
sherpa In []: plt.title('Counts and Background (red) for Negative- and Positive-Orders')
sherpa In []: plt.xlim(0, 180)
sherpa In []: plt.gca().set_xticks([]) # hide x-ticks

Now we populate the lower plot with the plus-order counts and background spectra in the same manner

sherpa In []: plt.subplot(2,1,2)
sherpa In []: plt.plot(get_data_plot(2).x,get_data_plot(2).y, marker="None",color="black")
sherpa In []: 
sherpa In []: bs2=get_bkg_scale(2)
sherpa In []: ct2=get_counts(2,bkg_id=1) # 1st (only) background associated with dataset 2
sherpa In []: plt.plot(get_data_plot(2).x,(ct2*bs2),marker="None", color="red")
sherpa In []: 
sherpa In []: plt.xlabel(r"$\lambda (\AA)$")
sherpa In []: plt.ylabel("Counts")
sherpa In []: plt.figtext(0.93,0.35,"Positive-Order",rotation=270)
sherpa In []: plt.xlim(0, 180)
sherpa In []: plt.ylim(0,35)

The results of this plot can be seen in Figure 5.

Figure 5: Spectra of Field Sources with Background Spectra over-plotted in red

[Thumbnail image: minus- and plus-order counts and background spectra of field]

[Version: full-size]

[Print media version: minus- and plus-order counts and background spectra of field]

Figure 5: Spectra of Field Sources with Background Spectra over-plotted in red

Strip chart of the spectra of counts (white) and background (red) with negative orders in the upper frame and postive orders in the lower frame.

Before exiting Sherpa, you may want to write the plot to a file. In this example, we use the interactive print dialog.

sherpa In []: plt.savefig("letgs_spec_bg.png")

Events in Wavelength Coordinates

To confirm that the edits made to the region were well-chosen, we inspect the events in wavelength coordinates. With the event file displayed in DS9:

unix% ds9 hrcf00029N006_evt2.fits &

we change the parameters to bin on wavelength coordinates (tg_lam,tg_d) (Bin → Binning Parameters). Again, we have to fiddle with our binning parameters, but we have some starting points. We can keep the tg_d binning of 0.00015 since the cross-dispersion ranges between -1 to 1 degree over about 10000 bins. The nominal wave-band coverage for the LETGS when using HRC-S is ≅150-175 Å. In wavelength, LETGS has 0.015 Å/bin; binning up by a factor of 10, we can get a good plot.

Figure 6: Wavelength Coordinates Binning Parameters

[DS9 menu settings for binning the source image and converting it to wavelength coordinates.]
[Print media version: DS9 menu settings for binning the source image and converting it to wavelength coordinates.]

Figure 6: Wavelength Coordinates Binning Parameters

Binning parameter settings applied in DS9 by going to the menu and choosing Bin → Binning Parameters.

Then load the region which is stored in the PHA file, a_01Lp1_pha2.fits[REGION]. Figure 7 shows the events and extraction regions in wavelength coordinates.

Figure 7: Grating Spectrum of the Two Field Sources in Wavelength Coordinates

[Thumbnail image: The grating spectrum of the two field sources after new binning parameters are applied, region from pha2 file overlayed.]

[Version: full-size]

[Print media version: The grating spectrum of the two field sources after new binning parameters are applied, region from pha2 file overlayed.]

Figure 7: Grating Spectrum of the Two Field Sources in Wavelength Coordinates

This is the field image in wavelength coordinates with tg_lam (wavelength) on the horizontal axis, and tg_d (cross-dispersion) on the vertical axis. The region from a_01Lp1_pha2.fits is overlaid, showing the source and background.


Scripting It

The file, tgextract2.py, performs the Sherpa commands used above.



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


##
## create 1D spectrum(a) table file(s) from the L1.5 output event list
##
        infile = hrcf00029N006_evt2.fits Input event file (output event file from L1.5 processing)
       outfile = a_               Output pha file (for pha2, enter full file name or '.'; for pha1, enter rootname.)
#
# tg_srcid_list parameter explanation:
#  - "all, blank or none" processes all the sources id's found in the event list.
#  - a comma list is a comma separated string list of all the sources to process. 
#          eg.  "1,2,5,7"
#  - @file is a pointer to an ascii file which contains the id's to process. 
#    Limit one id per line. no comma list is allowed in @file.
#
 tg_srcid_list = 1                Source ID's to process: 'all', comma list, @file
  tg_part_list = header_value     Grating parts to process: HETG, HEG, MEG, LETG, header_value
#
# tg_order_list parameter explanation :
#  - "default, blank or none" is set to process the following :
#          for ACIS :  1, 2, 3, -1, -2, -3
#          for HRC :   -1, 1
#  - a comma list is a comma separated string list of the orders
#    the user wants to process.
#          eg.  "-5, -1, 1, 3"
#  - a range list sets the min and max of the orders to process.
#    all the orders in between, will be processed. 
#          eg.  "-1..5"  will process orders from -1 to +5th order
#    a range list can be mixed with comma separated list
#  - @file is a pointer to an ascii file which contains 
#    the orders to process. Limit one order per line.
#    a range list is allowed in @file, but no comma list.
#
 tg_order_list = default          Grating diffraction orders to process: 'default', comma list, range list, @file
          (opt = pha1)            Output file type: pha1 (single spectrum) or pha2 (multiple spectra)
#
# If wav_grid is set to a binning specification (other than "use_header"), then 
# the following grating-specific parameters 'wav_grid_*' will be ignored.
#
     (wav_grid = use_header)      grid specification (use_header| min:max:step| min:max:#bins|pre-defined standard grid)
 (wav_grid_heg = 1.0:21.48:0.0025) pre-defined standard HEG grid
 (wav_grid_meg = 1.0:41.96:0.0050) pre-defined standard MEG grid
 (wav_grid_leg = 1.0:205.80:0.0125) pre-defined standard LETG grid
(wav_grid_leg_acis = 1.0:103.4:0.0125) pre-defined standard LETG/ACIS grid
   (evt_filter = none)            Filter to apply to events for counts spectra
(ignore_source_id = yes)             match source id between regionFile and evtFile?
        (error = gaussian)        Method for error determination (gaussian|gehrels)
  (region_file = tg_mask.fits)    Input region file ( NONE | none | <filename>)
      (geompar = geom)            Parameter file for Pixlib Geometry files
     (min_tg_d = default)         Minimum tg_d for source spectrum, in degrees
     (max_tg_d = default)         Maximum tg_d for source spectrum, in degrees
(extract_background = yes)             Extract the local background spectrum?
(min_upbkg_tg_d = default)         Minimum tg_d for upper background spectrum, in degrees
(max_upbkg_tg_d = default)         Maximum tg_d for upper background spectrum, in degrees
(min_downbkg_tg_d = default)         Minimum tg_d for down background spectrum, in degrees
(max_downbkg_tg_d = default)         Maximum tg_d for down background spectrum, in degrees
(backscale_method = events)          Use events and region, or region-only for backscale
(backscale_resolution = 64)              Amount in pixels to bin in wavelength to improve statistics in backscale
      (clobber = no)              overwrite any existing output file?
      (verbose = 0)               Verbosity level (0 = no display)
         (mode = ql)              



History

20 Apr 2009 New for CIAO 4.1
01 Feb 2010 reviewed for CIAO 4.2: no changes
19 Jul 2010 the S-Lang syntax has been removed from this thread as it is not supported in CIAO 4.2 Sherpa v2.
12 Jan 2011 reviewed for CIAO 4.3: no changes
10 Jan 2012 reviewed for CIAO 4.4: minor update to tgextract2_crates for the updated pycrates module.
03 Dec 2012 Review for CIAO 4.5. Updated sherpa command to use the get_bkg_scale routine instead of getting bscale and expsoure separately.
24 Apr 2013 Updated after the release of the tgmask2reg and reg2tgmask scripts
11 Dec 2013 Review for CIAO 4.6; no changes.
23 Dec 2014 Review for CIAO 4.7; no changes.
08 Apr 2019 Update to use matplotlib for plotting.
15 Feb 2022 Review for CIAO 4.14. No updates needed.