Search for Variability in a Source
CIAO 4.16 Science Threads
Overview
Synopsis:
The glvary tool searches for variability using the Gregory-Loredo algorithm. It splits the events into multiple time bins and looks for significant deviations between them. The tool assigns a variability indexed based on the total odds ratio, the corresponding probability of a variable signal, and the fractions of the lightcurve which are within 3σ and 5σ of the average count rate.
For an in-depth discussion of the variability algorithm, refer to the Gregory-Loredo Variability Probability why topic in the Chandra Source Catalog website.
Purpose:
To determine whether a source is variable.
Last Update: 15 Feb 2022 - Review for CIAO 4.14. Updated for Repro-5/CALDB 4.9.6.
Contents
- Get Started
- Select a Source
- Compute the Fractional Area
- Run glvary
- Examining the Output
- Comparison: running glvary without dither_region
- Parameter files:
- History
- Images
Get Started
Download the sample data: 635 (ACIS-I, RHO OPH CORE)
unix% download_chandra_obsid 635 evt2,asol,msk,bpix
Select a Source
We need to define a source region to use in the analysis. To manually create the source region, first display the image:
unix% ds9 acisf00635N005_evt2.fits &
In this example, we chose a source on ACIS-I3 (ccd_id=3), as shown in Figure 1.
[Version: full-size]
Figure 1: Source region on the event file
To save the regions, follow these steps:
- Region → Save Regions... → src.reg
- After choosing "OK" in the region filename dialog, a format dialog is opened. Set the format to "CIAO" and the coordinate system to "Physical".
The resulting region file looks like:
unix% more src.reg # Region file format: CIAO version 1.0 circle(4161.3752,3573.2228,10)
Compute the Fractional Area
A normalized effective area file is needed for input to glvary. This file is created by the dither_region tool. dither_region takes a region on the sky and computes its location on the detector. When part of the region falls onto a bad-pixel, goes off the detector, or goes outside the bounds of the mask (optional) the area is decremented. Why this file is necessary is explained in further detail in the "Comparison: running glvary without dither_region" section.
dither_region also takes bad pixels into account. Users need to be sure that the ardlib parameter file is set before running dither_region
unix% punlearn ardlib unix% acis_set_ardlib acisf00635_000N005_bpix1.fits abs- Updated ardlib parameter file: /home/user/cxcds_param4/ardlib.par AXAF_ACIS0_BADPIX_FILE -> acisf00635_000N005_bpix1.fits[BADPIX0] AXAF_ACIS1_BADPIX_FILE -> acisf00635_000N005_bpix1.fits[BADPIX1] AXAF_ACIS2_BADPIX_FILE -> acisf00635_000N005_bpix1.fits[BADPIX2] AXAF_ACIS3_BADPIX_FILE -> acisf00635_000N005_bpix1.fits[BADPIX3] AXAF_ACIS4_BADPIX_FILE -> CALDB AXAF_ACIS5_BADPIX_FILE -> CALDB AXAF_ACIS6_BADPIX_FILE -> acisf00635_000N005_bpix1.fits[BADPIX6] AXAF_ACIS7_BADPIX_FILE -> CALDB AXAF_ACIS8_BADPIX_FILE -> CALDB AXAF_ACIS9_BADPIX_FILE -> CALDB
The tool requires the aspect solution for the observation and the source region as inputs. Note that if there are multiple asol1.fits file for your observation, all of them must be provided, either as a list or as a stack. The mask file is necessary for the tool to know which CCDs were used in the observation.
unix% punlearn dither_region
unix% pset dither_region infile=pcadf00635_000N001_asol1.fits
unix% pset dither_region region="region(src.reg)"
unix% pset dither_region maskfile=acisf00635_000N005_msk1.fits
unix% pset dither_region outfile=fracarea.fits
unix% dither_region
Input aspect solution or histogram file(s) (pcadf00635_000N001_asol1.fits):
Region specification (region(src.reg)):
Output file name (fracarea.fits):
The output file, fracarea.fits, contains a fractional area for each row in the input file. The STATUS column indicates why the fraction for that row is less than 1; see the help file for an explanation of the status bits.
unix% dmlist fracarea.fits cols -------------------------------------------------------------------------------- Columns for Table Block AREAFRACTION -------------------------------------------------------------------------------- ColNo Name Unit Type Range 1 TIME s Real8 72039163.6412879974: 72141147.5575280041 Time 2 EQPOS(RA,DEC) deg Real8 -360.0: 360.0 Sky Position 3 ROLL deg Real8 -Inf:+Inf Roll angle 4 FRACAREA Real8 -Inf:+Inf Fraction area 5 AREA_CHIP_FRAC[10] Real8(10) -Inf:+Inf Region Area Fraction per chip 6 DELTA_T s Real8 -Inf:+Inf Time 7 CHIP_FRAC_TIME[10] Real8(10) -Inf:+Inf Fraction of ontime per chip 8 STATUS Bit[8] Why fraction < 1
The contents of the parameter file may be checked with plist dither_region.
Figure 2: Fraction of Region Area vs. Time
Run glvary
To run glvary, just the event file and the fractional area file need to be provided. The source region file is the same one used to run dither_region; a ccd_id filter is added to ensure that the correct good time information (GTI) is used.
The FRACAREA column in the dither_region only accounts for the geometric area of the region. In order to produce a properly normalized lightcurve, the dead time correction must also be included in the efficiency factor. For ACIS, this is the DTCOR value in the header of the event file.
unix% dmkeypar acisf00635N005_evt2.fits DTCOR echo+ 0.98733739787229 unix% dmtcalc "fracarea.fits[cols time,fracarea]" dtf_fracarea.fits \ expression="dtf=(0.98733739787229*fracarea)" clob+
In this example, we choose to have the tool output a probability-weighted lightcurve (lc_prob.fits) along with the table of output probabilitities (gl_prob.fits):
unix% punlearn glvary
unix% pset glvary infile="acisf00635N005_evt2.fits[sky=region(src.reg),ccd_id=3]"
unix% pset glvary effile=dtf_fracarea.fits
unix% pset glvary outfile=gl_prob.fits
unix% pset glvary lcfile=lc_prob.fits
unix% glvary
Input file specification (acisf00635N005_evt2.fits):
Output: probabilities as a function of m (gl_prob.fits):
Output: resulting light curve (lc_prob.fits):
Input file efficiency factors (fracarea.fits[cols time, dtf=fracarea]):
The contents of the parameter file may be checked with plist glvary.
Examining the Output
There are two output files from glvary: a probability-weighted lightcurve (lc_prob.fits) along with the table of output probabilitities (gl_prob.fits).
The variability index is calculated as:
Var. Index | Condition | Result |
---|---|---|
0 | P <= 0.5 | Definitely not variable |
1 | 0.5 < P < 2/3 AND f3 > 0.997 AND f5 = 1.0 | Considered not variable |
2 | 2/3 <= P < 0.9 AND f3 > 0.997 AND f5 = 1.0 | Probably not variable |
3 | 0.5 <= P < 0.6 | May be variable |
4 | 0.6 <= P < 2/3 | Likely to be variable |
5 | 2/3 <= P < 0.9 | Considered variable |
6 | 0.9 <= P AND Odd < 2.0 | Definitely variable |
7 | 2.0 <= Odd < 4.0 | Definitely variable |
8 | 4.0 <= Odd < 10.0 | Definitely variable |
9 | 10.0 <= Odd < 30.0 | Definitely variable |
10 | 30.0 <= Odd | Definitely variable |
Refer to the glvary help file for details.
The probability, odds, and resulting variability index are all recorded in the header of the odds ratio file:
unix% dmlist gl_prob.fits header | grep -i variab 0008 ODDS 155.1249009276 Real8 Odds for variable signal 10Log 0009 PROB 1.0 Real8 Probability of variable signal 0014 VARINDEX 10 Int4 Variability index
The variability index value of 10 indicates that this source is definitely variable.
The lightcurve consists of the binnings weighted by the odd ratios and shows the optimal binning for the curve. The standard deviation is provided for each point on the lightcurve. The lightcurve can be plotted to visualize the data:
unix% python >>> from pycrates import read_file >>> import matplotlib.pylab as plt >>> >>> tab = read_file("lc_prob.fits") >>> xx = tab.get_column("Time").values >>> yy = tab.get_column("COUNT_RATE").values >>> ye = tab.get_column("COUNT_RATE_ERR").values >>> >>> plt.errorbar(xx,yy,yerr=ye, ecolor="red", color="black") >>> plt.xlabel("Time (s)") >>> plt.ylabel("Count Rate (count/s)") >>> plt.title("Rho Oph Core") >>> plt.show()
Figure 3 show the resulting plot.
[Version: full-size]
Figure 3: Lightcurve created by glvary
Comparison: running glvary without dither_region
When extracting a lightcurve in a given region using dmextract, the count rate is given by the total counts encircled by the region divided by the good time in the time bin. No account is taken of the fraction of the region that is on or off the chip at any given moment. Due to the dither motion of the spacecraft, the extraction region can pass over the edge of the chip or over bad pixels and columns. Thus the detected counts from the source can vary on dither time scales solely due solely to the changing effective area, rather than any intrinsic variations, as shown in Figure 4. Such purely instrumental variations will be reflected in lightcurves created with dmextract. When searching for variability and creating a lighcurve with the glvary tool, we wish to minimize such purely instrumental signatures.
[Version: full-size]
Figure 4: Lightcurves illustrating the effect of the fractional area file
To this end, the glvary tool utilizes information from the dither_region tool to create a simple correction for these instrumental effects. Rather than basing the variability estimate on the extracted counts per good time in a time bin, glvary uses the extracted counts per good time in a bin per normalized extraction area. The dither_region tool calculates the fractional area in the extraction region as a function of time, accounting for chip edges, bad pixels, and bad columns. Not only is this geometric correction used by glvary in calculating the significance of the variability in a lightcurve, it is also applied to correct the rate in the estimated lightcurve produced by the tool. Thus for sources that dither over chip edges or bad pixels, the glvary lightcurve will almost always show less variability than the comparable lightcurve for the same region extracted with dmextract.
Note that neither glvary nor the dither_region tool account from more subtle instrumental variations. For example, if the source dithers over chip regions with differing amounts of contamination, or over chip regions with other effective area changes not due to bad pixels or chip edges, these instrumental variations can be reflected in the extracted lightcurves with both dmextract and glvary.
Parameters for /home/username/cxcds_param/dither_region.par infile = pcadf00635_000N001_asol1.fits Input aspect solution or histogram file(s) region = region(src.reg) Region specification outfile = fracarea.fits Output file name (maskfile = acisf00635_000N005_msk1.fits) Mask file (psffile = ) PSF Image file (gtifile = ) GTI File (dtffile = ) DTF File (wcsfile = ) WCS File (imapfile = ) Stack of Instrument files (tolerance = INDEF) Tolerance of aspect solution [arcsec] (resolution = 1) Binning resolution of region [pixels] (maxpix = 1000) Maximum number of pixels regardless of resolution (convex = no) Use convex hull around aspect histogram? (geompar = geom) Parameter file for Pixlib Geometry files (ardlibpar = ardlib) Parameter file for ARDLIB files (detsubsysmod = ) Detector sybsystem modifier (verbose = 0) Tool verbosity (clobber = no) Remove outfile if it already exists (mode = ql)
Parameters for /home/username/cxcds_param/glvary.par infile = acisf00635N005_evt2.fits Input file specification outfile = gl_prob.fits Output: probabilities as a function of m lcfile = lc_prob.fits Output: resulting light curve effile = fracarea.fits[cols time, dtf=fracarea] Input file efficiency factors (probfile = NONE) Input probability file for background (frac = 1.0) Fraction of events to be included in subsample (seed = 1) Seed for random subsample selection (mmax = INDEF) Maximum number of model bins (mmin = INDEF) Minimum number of model bins (nbin = 0) Number of bins to use in light curve (mintime = 50) Range of binnings, maximum resolution in seconds (clobber = no) Overwrite output files if they exist? (verbose = 0) Tool chatter level (mode = ql)
History
20 Apr 2009 | New for CIAO 4.1.2 |
07 Jul 2009 | the mask file is necessary when running dither_region so that the tool knows which CCDs were used in the observation |
27 Jan 2010 | updated for CIAO 4.2: updated DS9 instructions |
12 Jan 2011 | reviewed for CIAO 4.3: no changes |
10 Jan 2012 | reviewed for CIAO 4.4: no changes |
03 Dec 2012 | Review for CIAO 4.5; updated file name versions |
13 Dec 2013 | Review for CIAO 4.6; added plot showing dither_region output. |
23 Dec 2014 | Review for CIAO 4.7; no changes. |
13 Jan 2016 | Review for CIAO 4.9. Added badpixel file to retrieval list and note to be sure ardlib is set before running dither_region. Also added step to scale FRACAREA by DTF to get correct normalize of the lightcurve. |
08 Apr 2019 | Updated to use matplotlib for plotting. |
15 Feb 2022 | Review for CIAO 4.14. Updated for Repro-5/CALDB 4.9.6. |