# Applying Customized Background Regions to LETG/HRC-S Observations

## 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.

Last Update: 8 Apr 2019 - Update to use matplotlib for plotting.

## 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

unix% download_chandra_obsid 29 evt2

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. Next, we set the binning parameters in the DS9 menu Bin → Binning Parameters (Figure 2): 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): ## 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). 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.

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.

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.

## Scripting It

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



##
## 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_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.