Filtering out dark events by Pulse-Height Amplitude (PHA)

This tutorial outlines the steps needed to restrict the PHA filtering in COS datasets. This involves examining the datasets, modifying a reference file, and re-calibrating the data. The standard values are currently 2-23, meaning only events within this range are included in final processing. We will see how and why to restrict this below.

In [42]: 
import os
import glob
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from astropy.io import fits
In [43]: 
#-- iPython notebook magic command to make plots appear within the 
#-- browser instead of opening a separate plotting window.
%matplotlib inline

#-- Increase the default fontsize for the plots
mpl.rcParams['font.size'] = 16
In [44]: 
#-- CalCOS needs reference files to function.
#-- this will set the lref environment variable if you don't have it set.
#--
#-- This is the correct location for users internal to STScI, you may have to change it.
os.environ['lref'] = '/grp/hst/cdbs/lref/'

Pulse Height Distribution

To see what kind of filtering we may be able to do, we need to first look at the detected events in PHA space. To do this, we create histograms of the PHA column in the corrtag for FUVA and FUVB data. The range of PHAs in data is integer values from 0-31, so these will be our bins.

Below, we can see the PHD for each segment from this particular dataset. The standard filtering limits are plotted as red lines. What you see in the distribution is a strong feature from the external target itsself superimposed with a lower-level, mostly flat, distribution of dark events. These dark events come from a variety of sources: cosmic rays, double-counted events, and detector-intrinsic noise. Each of these sources contribute a different shape to the background distribution, but as a 1st order approximation the filtering boundaries can be moved in from their default values until they start to touch the source distribution. This will ensure that you are not including areas of Pulse-height space that are dark-dominated. The approximate regions that can additionally be filtered for this dataset are shown in grey.

In [45]: 
fig, (ax1, ax2) = plt.subplots(2, figsize=(7, 7), sharex=True)

#-- Create a histogram of the FUVA events in PHA space
hdu = fits.open('data/lbry01icq_corrtag_a.fits')
ax1.hist(hdu['events'].data['PHA'], bins=range(0, 32))
ax1.set_ylabel('FUVA Distribution')
ax1.axvspan(xmin=2, xmax=2, ymin=0, ymax=1, color='grey', alpha=.5)
ax1.axvspan(xmin=19, xmax=23, ymin=0, ymax=1, color='grey', alpha=.5)

#-- Create a histogram of the FUVB events in PHA space
hdu = fits.open('data/lbry01icq_corrtag_b.fits')
ax2.hist(hdu['events'].data['PHA'], bins=range(0, 32))
ax2.set_xlabel('PHA')
ax2.set_ylabel('FUVB Distribution')
ax2.axvspan(xmin=15, xmax=23, ymin=0, ymax=1, color='grey', alpha=.5)

for axis in [ax1, ax2]:
    axis.axvline(x=2, color='r', ls='--', lw=3, label='Lower Bound')
    axis.axvline(x=23, color='r', ls='--', lw=3, label='Upper Bound')

fig.subplots_adjust(hspace=0)
fig.savefig('filtering_by_segment.pdf', bbox_inches='tight')

Modifying the reference file

The PHA limits are set in the PHATAB reference file. We can modify this very simply.

In [46]: 
print hdu[0].header['phatab']
print
#-- Patch together LREF and the phatab filename to open
phatab = fits.open(os.path.join(os.environ['lref'], hdu[0].header['phatab'].split('$')[1]))

phatab.info()

lref$ucl1623gl_pha.fits

Filename: /grp/hst/cdbs/lref/ucl1623gl_pha.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      29   ()              
1                BinTableHDU     29   6R x 6C      [8A, 4A, 1J, 1J, 1E, 1E]
In [47]: 
phatab[1].data
Out[47]:
FITS_rec([('G130M', 'FUVA', 2, 30, 0.1, 1.9),
       ('G130M', 'FUVB', 2, 30, 0.1, 1.9),
       ('G160M', 'FUVA', 2, 30, 0.1, 1.9),
       ('G160M', 'FUVB', 2, 30, 0.1, 1.9),
       ('G140L', 'FUVA', 2, 30, 0.1, 1.9),
       ('G140L', 'FUVB', 2, 30, 0.1, 1.9)], 
      dtype=[('OPT_ELEM', 'S8'), ('SEGMENT', 'S4'), ('LLT', '>i4'), ('ULT', '>i4'), ('MIN_PEAK', '>f4'), ('MAX_PEAK', '>f4')])

The table is broken up by grating and segment. LLT is the lower threshold, while ULT is the upper threshold. Lets set all the FUVA dataset to filter from 2 to 19, and FUVB to filter from 2 to 15.

In [48]: 
index = np.where(phatab[1].data['SEGMENT'] == 'FUVA')
phatab[1].data['ULT'][index] = 19

index = np.where(phatab[1].data['SEGMENT'] == 'FUVB')
phatab[1].data['ULT'][index] = 15

Ok, we can see we modified the correct column, now we can just write it out to a new filename.

In [49]: 
phatab[1].data
Out[49]:
FITS_rec([('G130M', 'FUVA', 2, 19, 0.1, 1.9),
       ('G130M', 'FUVB', 2, 15, 0.1, 1.9),
       ('G160M', 'FUVA', 2, 19, 0.1, 1.9),
       ('G160M', 'FUVB', 2, 15, 0.1, 1.9),
       ('G140L', 'FUVA', 2, 19, 0.1, 1.9),
       ('G140L', 'FUVB', 2, 15, 0.1, 1.9)], 
      dtype=[('OPT_ELEM', 'S8'), ('SEGMENT', 'S4'), ('LLT', '>i4'), ('ULT', '>i4'), ('MIN_PEAK', '>f4'), ('MAX_PEAK', '>f4')])
In [50]: 
phatab.writeto('data/modified_phatab.fits')

Now you simply need to set the PHATAB header keyword in the raw data to point to this new file, and run CalCOS again.

In [50]: