Filtering Lightcurves

Get Started

Sample ObsID used: 1712 (ACIS-S, 3C 273)

File types needed: evt2

This thread uses the lc_sigma_clip routine from the lightcurves.py module. In CIAO 4.0 and earlier, the S-Lang routine was called analyze_ltcrv and the script analyze_ltcrv.sl (there was no Python version). The CIAO 4.1 version restores the diagnostic plots that were not present in the CIAO 4.0 version, as well as adding new capabilities (such as the ability to create a GTI file).

The scripts are part of the CIAO Scripts distribution. The CIAO scripts package should be the following version or newer:

unix% cat $ASCDS_CONTRIB/VERSION.CIAO_scripts
17 Apr 2009

Please check that you have at least this version of the scripts package installed before continuing. If you do not have the scripts installed or need to update to a newer version, refer to the Scripts page.

For this thread we shall restrict attention to the ACIS-S3 chip, and the 0.5 to 7 keV energy range:

unix% punlearn dmcopy
unix% dmcopy "acisf01712N003_evt2.fits[energy=500:7000,ccd_id=7]" evt2_c7.fits 

An image of the data is shown in Figure 1.

[Version: full-size]

Figure 1: ACIS-S3 observation of a field

ACIS-S3 observation in the 0.5 to 7 keV energy range.

Remove bright/variable sources from the dataset

To avoid any background variations due to sources in the field we first filter out regions of high/variable emission from the dataset. There are a number of ways to create a list of regions that should be excluded - for instance you may wish to use the output from one of the source detection programs. In this thread we shall use visual estimation, as an example.

Using ds9, we chose several regions and saved them in CIAO format as exclude.reg, which looks like:

unix% cat exclude.reg
# Region file format: CIAO version 1.0
rotbox(4200.3328,4137.9892,1129.5056,74.07019,24.22333)
circle(4076.5,4088.5,316)
circle(4296.5,5024.5,48)

The chosen regions are shown in Figure 2 (note that we ignored several obvious point sources in this example).

[Version: full-size]

Figure 2: "Source" regions

For this example, the regions were chosen to contain most of the source emission.

Create the lightcurve (dmextract)

We can now create a lightcurve of this background region using the CIAO dmextract tool. The best way to bin the lightcurve depends on the dataset and your scientific objectives; for this example we use a bin length of 200 seconds:

unix% punlearn dmextract
unix% pset dmextract "evt2_c7.fits[exclude sky=region(exclude.reg)][bin time=::200]" 
unix% pset dmextract outfile=lc_c7.fits
unix% pset dmextract opt=ltc1
unix% dmextract
Input event file  (evt2_c7.fits[exclude sky=region(exclude.reg)][bin time=::200]): 
Enter output file name (lc_c7.fits): 

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

Analyze the lightcurve using lc_sigma_clip()

The lc_sigma_clip() routine provided by the lightcurves.py script performs an iterative sigma-clipping algorithm, removing those points that fall outside +/-3 sigma from the mean at each iteration until all data points are within +/-3 sigma (the number of sigma used to clip the data defaults to 3 but can be changed). This algorithm is robust but not perfect; it can easily "overclean" a noisy lightcurve and should not be used blindly. The routine can print out the limits used and - new to CIAO 4.1 - can also create a GTI file that encodes these limits.

The script must be loaded into chips, sherpa, or ipython before it can be run (this step is only necessary once per session). In the examples below we load the script into a ChIPS session:

chips> from lightcurves import *

If called with one argument - the name of the lightcurve to analyze - then the routine will use a 3-sigma clip, print out the resultant filter, and create a plot of the data (Figure 3):

chips> lc_sigma_clip("lc_c7.fits")
Parameters used to clean the lightcurve are:
  script version = CIAO 4.1 - 1.0
  sigma          = 3
  minlength      = 3
  plot           = True
  rateaxis       = x
  color          = green

Total number of bins in lightcurve   = 153
Max length of one bin                = 197.467 s
Num. bins with a smaller exp. time   = 2
Num. bins with exp. time = 0         = 13
Number of bins with a rate of 0 ct/s = 13

Rate filter:  0.127776 <= count_rate < 1.53282 
Mean level of filtered lightcurve = 0.830298 ct/s
GTI limits calculated using 2 time blocks
    ((time >= 77379470.949928) && (time < 77399470.949928)) ; 20.00 ksec
    ((time >= 77404870.949928) && (time < 77406670.949928)) ;  1.80 ksec

[Version: full-size, PNG, postscript]

Figure 3: Plot from lc_sigma_clip

The top plot shows the light curve of lc_c7.fits, where the time axis is displayed on the x axis as a relative value (the first bin in the light curve is used as the 0 point). Points that are in green are those that passed the sigma-clipping and the horizontal dotted line shows the mean value of these points (the numerical value is given in the plot title).

The bottom plot shows a histogram of the count rate values; the green version again shows the data that passes the sigma-clipping algorithm.

We change the Y axis limits in the top plot to select the "good" count rates, using the approximate range of 0.1 to 1.6 count/s (since the filter limits are given as 0.127776 - 1.53282 in the screen output). Note that this ChIPS command also changes the X-axis range in the bottom plot since it has bound to the Y axis in the top plot, as shown in Figure 4:

chips> print info()
Window [lc_plot]
  Frame [lc_sigma_clip]
    Plot [rate]   (0.15,0.56)  .. (0.90,0.90)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      X Axis [ax1]
      Y Axis [ay1]
      Curve [bad]
      Curve [good]
      Line [hline]
    Plot [hist]   (0.15,0.15)  .. (0.90,0.44)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      X Axis [ax1]
      Y Axis [ay1]
      Histogram [everything]
      Histogram [good]
      Line [vline]
    Plot [labels]   (0.15,0.15)  .. (0.90,0.90)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      Label [filename]
      Label [method]
      Label [object]
      Label [obsid]

chips> current_plot("rate")
chips> limits(Y_AXIS, 0.1, 1.6)

[Version: full-size, PNG, postscript]

Figure 4: Viewing the "good" count rate periods

Changing the Y-axis limits to 0.1 to 1.6 count/s of the top plot ("hist") zooms both plots from Figure 3. The asymmetry of the distribution of count rates around the calculated mean value is easily seen here.

A number of options are available to control the behavior of lc_sigma_clip, as described in "ahelp lc_sigma_clip". Below we change the rateaxis option to swap the axes in the top plot, and set verbose to 0 to remove the screen output. The new plot is shown in Figure 5:

chips> lc_sigma_clip("lc_c7.fits", rateaxis="x", verbose=0)

[Version: full-size, PNG, postscript]

Figure 5: Plot from lc_sigma_clip with rateaxis="x"

As rateaxis was set to "x", the top plot has the count rate displayed on the X axis rather than the Y axis (Figure 3).

The count rate axes are still linked together, as can be seen if you enter the following commands (the top plot is named rate and the bottom plot hist):

chips> current_plot("rate")
chips> limits(X_AXIS, 0.1, 1.6)

Filtering the event file

The output time periods can be used to filter the event list. The lc_sigma_clip routine can create a GTI file directly, or you can use the limits it calculates to create the file yourself using the dmgti tool, or use them directly within a DM filter expression, as described below.

1. Using lc_sigma_clip

New to CIAO 4.1 is the ability of lc_sigma_clip to directly create a GTI file, by setting the outfile parameter:

chips> lc_sigma_clip("lc_c7.fits", outfile="lc_c7.gti")
Parameters used to clean the lightcurve are:
  script version = CIAO 4.1 - 1.0
  sigma          = 3
  minlength      = 3
  outfile        = lc_c7.gti
  plot           = True
  rateaxis       = x
  color          = green
  pattern        = crisscross
  pattern color  = red

Total number of bins in lightcurve   = 153
Max length of one bin                = 197.467 s
Num. bins with a smaller exp. time   = 2
Num. bins with exp. time = 0         = 13
Number of bins with a rate of 0 ct/s = 13

Rate filter:  0.127776 <= count_rate < 1.53282 
Mean level of filtered lightcurve = 0.830298 ct/s
GTI limits calculated using 2 time blocks
    ((time >= 77379470.949928) && (time < 77399470.949928)) ; 20.00 ksec
    ((time >= 77404870.949928) && (time < 77406670.949928)) ;  1.80 ksec
Creating GTI file
Created: lc_c7.gti

Using this option will change the plot from that shown in Figure 3, since the regions filtered out are now indicated by a solid red pattern as shown in Figure 6.

[Version: full-size, PNG, postscript]

Figure 6: Plot from lc_sigma_clip when outfile is set

The top plot now indicates the regions excluded by the GTI file created by lc_sigma_clip - here lc_c7.gti - using a solid red pattern.

If you do not want a plot created then you can set the plot option to False; you may also want to set the verbose option to 0 to remove the screen output: for example

chips> lc_sigma_clip("lc_c7.fits", outfile="lc_c7.lc.gti", plot=False, verbose=0)
chips>

See "ahelp lc_sigma_clip" for a full list of options.

The output, here lc_c7.lc.gti, can then be used to filter an event list:

unix% dmcopy "evt2_c7.fits[@lc_c7.lc.gti]" evt2_c7_lc.fits

---

2. Using dmgti

The time limits can be used by dmgti to create a GTI file, as shown below (note that the individual ranges are combined by "||"):

unix% punlearn dmgti
unix% pset dmgti infile=lc_c7.fits
unix% pset dmgti outfile=lc_c7.dmgti.gti
unix% pset dmgti \
userlimit="((time >= 77379470.949928) && (time < 77399470.949928))||((time >= 77404870.949928) && (time < 77406670.949928))"
unix% dmgti
Input MTL file (lc_c7.fits): 
Output GTI file (lc_c7.dmgti.gti): 
User defined limit string (((time >= 77379470.949928) && (time <
77399470.949928))||((time >= 77404870.949928) &&
(time < 77406770.949928))

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

The output of dmgti, here lc_c7.dmgti.gti, can then be used to filter an event list:

unix% dmcopy "evt2_c7.fits[@lc_c7.dmgti.gti]" evt2_c7_dmgti.fits

---

3. Using a DM filter

Alternatively, you can use filter the event list directly using the column filtering capabilities of the CIAO Data Model:

unix% dmcopy \
      "evt2_c7.fits[time=77379470.949928:77399470.949928,77404870.949928:77406670.949928]" \
      evt2_c7_dmcopy.fits

The final, filtered, exposure time can be found by looking at the EXPOSURE keyword of the filtered file:

unix% dmkeypar evt2_c7_lc.fits EXPOSURE echo+
21362.821808294

unix% dmkeypar evt2_c7_dmgti.fits EXPOSURE echo+
21362.821808294

unix% dmkeypar evt2_c7_dmcopy.fits EXPOSURE echo+
21362.821808294