Phase Resolved Spectroscopy
CIAO 4.15 Science Threads
Overview
Synopsis:
Creating a phase-resolved spectrum for a periodic source allows the user to compare the spectrum from different portions of the cycle and look for changes that may occur as the source rotates or orbits.
This thread has been completely re-written. Previous versions of this thread used a different technique to filter the event file on PHASE that left the user to compute the correct exposure time. This new technique provides a way to filter the data using good time intervals (GTIs) so that the Datamodel can compute the exposure time automatically. This version of the thread also provides some guidance on how to determine the period.
Purpose:
To create a phase-resolved spectrum for a source.
Related Links:
- MIT's SITAR package
- Timing Analysis with Lightcurves
- Basic Lightcurves
- Search for Variability in a Source
- Apply Barycenter Correction
- Determining Exposed Region Area
- Good Time Intervals
- A discussion on barycentric time corrections.
Last Update: 13 Feb 2023 - Added a note about change to default dither parameters for ACIS.
Contents
- Get Started
- Apply Barycenter Correction
- Determining the Period
- Compute Phase
- Extract Phase Resolved Spectrum
- Summary
- History
-
Images
- Figure 1: Image of OBS_ID 2795
- Figure 2: Light Curve for RXJ0806.3+1527
- Figure 3: Gregory-Loredo Lightcurve for RXJ0806.3+1527
- Figure 4: Power Spectrum of Lightcurve
- Figure 5: Period Fold Results
- Figure 6: Period Fold Lightcurve
- Figure 7: dmextract: 321.85 seconds period folded lightcurve
- Figure 8: GTI Alignment Corrections
- Figure 9: GTIs shown on lightcurve
- Figure 10: Spectrum for Phase 0.4:0.9 with 321.85 second period.
Get Started
Download the sample data: 2795 (ACIS-7, RXJ0806.3+1527)
unix% download_chandra_obsid 2795
Users should not fold a spectrum on a period that is close to or smaller than the temporal resolution of the data. The resolution is the TIMDEL value in the header of the event file.
unix% dmkeypar acisf02795_repro_evt2.fits TIMEDEL echo+ 0.44104
Users working with gratings data should pay extra attention to the dmcopy opt=all bug when working with TIME/GTI filters.
For the remainder of this thread it is assumed that the data have been reprocessed using chandra_repro and that all the files have been copied into the current working directory.
The Orbit Ephemeris file, orbitf121694700N001_eph1.fits, is included in the primary/ directory. Users with recently obtained datasets may only have access to the Level 0 ephemeris file (eph0.fits). This is the predicted ephemeris. The level 1 ephemeris (eph1.fits) is the definitive ephemeris which is generally available within about a month of the observation.
The data for this observation are shown in Figure 1. To simplify later commands, the event file is filtered with the region file and on energy
unix% dmcopy "acisf02795_repro_evt2.fits[sky=region(ciao.reg),energy=500:7000]" RXJ0806.3+1527.evt
[Version: full-size]
Figure 1: Image of OBS_ID 2795
The filtered version of the event file is used in the remainder of this thread.
Apply Barycenter Correction
The times in the data products need to be corrected to barycenter following the applying a barycenter correction thread. This is to account for the difference in photon arrival times as Chandra orbits the Earth and the Earth moves around the Sun. This barycenter correction should precede any phase calculation when trying to align phase between observations.
A summary of the steps in the thread are shown below.
Correct the event file and GTIs
unix% axbary RXJ0806.3+1527.evt orbitf121694700N001_eph1.fits RXJ0806.3+1527_bary.evt ra=121.5956933655199 dec=15.45862251788717
Correct the aspect solution file
unix% axbary pcadf121941470N004_asol1.fits orbitf121694700N001_eph1.fits RXJ0806.3+1527_bary.asol ra=121.5956933655199 dec=15.45862251788717
If the observation has multiple aspect solution files, this correction must be applied to each individually. The new ASOL file name is now added to the header of the event file so that later tools can automatically locate it.
unix% dmhedit RXJ0806.3+1527_bary.evt file= op=add key=asolfile value=RXJ0806.3+1527_bary.asol
Correct the Exposure Statistics file
For ACIS observations, there is one additional data product that must be barycenter corrected: the exposure statistics, stat1 file:
unix% axbary acisf02795_000N003_stat1.fits orbitf121694700N001_eph1.fits.gz RXJ0806.3+1527_bary.expstats ra=121.5956933655199 dec=15.45862251788717
This file will be used to correctly align the good time intervals (GTIs) with ACIS exposure frame time boundaries.
Determining the Period
Users who know the period of their source can skip directly to the Compute Phase section of this thread.
CIAO provides some tools to perform rudimentary timing analysis that can be used to determine the periodicity of a variable source. This section highlights those tools and provides some insight into how to interpret the outputs.
There are several papers in the literature about this source which all cite a period around 321 seconds:
which provide independent verification of the results presented here.
Light Curves
dmextract can be used to extract a light curve binned at a fixed time resolution. Below is an example using a bin size of 20 seconds. The output lightcurve is shown in Figure 2. Remember that the event file used here has already been filtered spatially and by energy range.
unix% dmextract "RXJ0806.3+1527_bary.evt[bin time=::20]" RXJ0806.3+1527.lc op=ltc1 clob+
Background has been ignored in this example due to the brightness of the source and the small region size.
[Version: full-size]
Figure 2: Light Curve for RXJ0806.3+1527
This lightcurve shows a strong periodic, what might be eclipsing, behavior. Users can learn more about lightcurves in the Basic Lightcurves thread.
The choice of bin size can sometimes mask variability or accentuate artificial/instrumental frequencies. The next tool uses a Bayesian blocks approach to determine the optimal bin size.
Gregory-Loredo Bayesian Blocks
The glvary tool creates a probability weighted, optimally binned lightcurve using the Gregory-Loredo Bayesian blocks algorithm.
unix% glvary RXJ0806.3+1527_bary.evt outfile=glvary.prob lcfile=glvary.lc effile="" clob+
The data are shown in Figure 3.
[Version: full-size]
Figure 3: Gregory-Loredo Lightcurve for RXJ0806.3+1527
Figure 3 closely resembles Figure 2. The glvary curve is more smooth since it is computed by summing probability weighted lightcurves at different resolutions. Users will also note that the curve does not go to zero which is a byproduct of the GL algorithm.
The Search for Variability in a Source thread shows more examples of how to use the glvary tool.
Having established that the source is periodic, attention now turns to determining the frequency.
Power Spectrum
The standard approach to determine the frequency for a periodic signal is to compute the power spectrum (amplitude of the Fourier transform). Users should in general be cautious when trying to interpret the power spectrum of a lightcurve for low count data; Poisson noise can easily mask faint signals. This dataset has a significant number of counts per bin so the power spectrum is meaningful.
The dmextract lightcurve should be used to create the power spectrum as the glvary lightcurve has been too heavily processed. Also, although not a concern here, dmextract can produce a background subtracted lightcurve while glvary cannot (at least not directly).
First, the mean count rate is subtracted from the lightcurve. This is to suppress the large, constant, spike at 0 frequency that can obscure low frequency components.
unix% dmstat "RXJ0806.3+1527.lc[cols count_rate]" COUNT_RATE[count/s] min: 0 @: 1 max: 0.5513 @: 265 mean: 0.072791138986 sigma: 0.10489306672 sum: 85.037711491 good: 1169 null: 0 unix% set meanval = `pget dmstat out_mean` unix% dmtcalc RXJ0806.3+1527.lc RXJ0806.3+1527_mean_offset.lc exp="count_rate=count_rate-${meanval}" clob+
Then the power spectrum is computed using the apowerspectrum tool.
unix% apowerspectrum infilereal="RXJ0806.3+1527_mean_offset.lc[cols time,count_rate]" infileimag=none outfile=RXJ0806.3+1527.power clob+ unix% dmlist RXJ0806.3+1527.power cols -------------------------------------------------------------------------------- Columns for Table Block POWERSPECTRUM -------------------------------------------------------------------------------- ColNo Name Unit Type Range 1 FREQUENCY (s)**-1 Real8 -Inf:+Inf 2 DATA Real4 -Inf:+Inf
Since the lightcurve is entirely real valued, use infileimag=none. The output file is simply Frequency versus Power (amplitude of the Fourier transform). The data are shown in Figure 4.
[Version: full-size]
Figure 4: Power Spectrum of Lightcurve
The power peaks at a frequency between 0.00310 and 0.00314 Hz (period between 318.5 and 322.6 seconds), with some higher frequency harmonics present. Those values agree well with the 321 second period reported in the literature.
To provide a better estimate of the frequency using the power spectrum approach would require using a smaller bin size with dmextract. However, the bin size needs to be significantly larger than the native time resolution (TIMEDEL=0.44 seconds) to avoid aliasing effects so this is about at the practical limit for this dataset.
Another common approach to locate the frequency of a periodic signal is to perform period folding (may also be called epoch folding). This is performed with the pfold tool.
Period Fold
Another way to refine the period is to search for periods in some specified range. The basic algorithm is to iterate over a grid of expected periods, create a phased-binned lightcurve for each period, and the locate which period shows the strongest deviation from a uniform/flat signal. Folding lowers the statistical uncertainty in the phase bins as the number of counts increase.
There are various metrics that can be used to measure strongest, ie the L-statistic used in the SITAR package. pfold currently only outputs the standard deviation of the count rate, SIGMA_RATE, together with the average count rate MEAN_RATE.
To ensure that all the products are phase aligned (using the same time offset for PHASE=0), the TSTART value from the event file header is explicitly being used.
unix% dmkeypar RXJ0806.3+1527_bary.evt tstart echo+ 121941081.557
The power spectrum identified a periodic signal with a period of ~320 second, so a grid of times around that period will be searched. This example uses a very generous range of periods from 300 to 400 seconds with a step size of 0.01 seconds.
unix% pfold RXJ0806.3+1527_bary.evt RXJ0806.3+1527.fold periodgrid=300:400:0.01 tzero=121941081.557 clobber=yes unix% dmlist RXJ0806.3+1527.fold cols -------------------------------------------------------------------------------- Columns for Table Block RXJ0806.3+1527.fold -------------------------------------------------------------------------------- ColNo Name Unit Type Range 1 PERIOD sec Real8 -Inf:+Inf Folded Period 2 PHASE_MID[10] Real8(10) -Inf:+Inf PHASE GRID 3 COUNTS[10] Int4(10) - Counts in each phase bin 4 RATE[10] counts/sec Real8(10) -Inf:+Inf Folded rates 5 ONTIME[10] sec Real8(10) -Inf:+Inf Exposure time 6 DTF[10] Real8(10) -Inf:+Inf Deadtime correction 7 MEAN_RATE counts/sec Real8 -Inf:+Inf Average rate 8 SIGMA_RATE counts/sec Real8 -Inf:+Inf Stdev of rate
The pfold output contains the phase binned lightcurve (PHASE_MID vs. RATE) for each of the input periods. Additionally the MEAN_RATE (average count rate) and SIGMA_RATE (standard deviation of count rate values) is reported; these values shown in Figure 5.
[Version: full-size]
Figure 5: Period Fold Results
The peak in the SIGMA_RATE column occurs at a period of 321.85 seconds. The phase binned lightcurve for that period is shown in the Figure 6 inset.
[Version: full-size]
Figure 6: Period Fold Lightcurve
Analysis of the pfold output shows that a period of 321.85 seconds should be used for to compute the PHASE.
The pfold array columns can be converted into normal scalar arrays using the dmtype2split tool. While this tool is normally used to convert gratings type II PHA files into standard type I files, it can be used more generically to convert any file with array columns into scalar columns.
unix% dmtype2split RXJ0806.3+1527.fold"[period=321.85]" 321.85.fits unix% dmlist 321.85.fits cols -------------------------------------------------------------------------------- Columns for Table Block 321.85.fits -------------------------------------------------------------------------------- ColNo Name Unit Type Range 1 PHASE_MID Real8 -Inf:+Inf 2 COUNTS Int4 - 3 RATE counts/sec Real8 -Inf:+Inf 4 ONTIME sec Real8 -Inf:+Inf 5 DTF Real8 -Inf:+Inf
This file format may be easier for other applications to use.
Compute Phase
In the previous section several types of analysis were performed to come to the conclusion that RXJ0806.3+1527 has a period of 321.85 seconds. This number will now be used to create a set of Good Time Intervals (GTIs) and to extract a spectra from the brightest part of the cycle.
Period fold expression
This thread will be using dmtcalc to compute the phase (fraction of the period) of times in various data products. dmtcalc can compute arbitrary expressions via its rich expression syntax.
The expression to compute the phase is shown below and saved into a file named dmtcalc.expr
unix% cat dmtcalc.expr .period=321.85 .tzero=121941081.5572 .freq=(1.0/.period) GPHASE=((time-(.tzero))*.freq) PHASE=(GPHASE-((long)GPHASE))
The leading dot (period) for the columns .period, .tzero, and .freq designate them as temporary columns which are not written to the output file.
If the frequency was known, instead of the period, then the .freq column can be set directly.
Folding the observation entails associating all the times that happen at a given point in the phase cycle with one another, with the phase traditionally being defined as a fraction of the orbital or rotational cycle < 1. To do this, the integer portion of the GPHASE must be subtracted off so that all values are greater than or equal to 0 and less than 1 (e.g. phases of 0.5, 1.5, and 2.5 all fall at identical points in the cycle):
Previous versions of this thread used the header keyword TSTART directly
GPHASE=((time-TSTART)*.freq)
This is a problem because the TSTART value is different in different data products. Be sure to use the same tzero value with all the files.
This dmtcalc.expr file can now be used as input to dmtcalc. First, the phase for the times in the event file are computed.
unix% dmtcalc RXJ0806.3+1527_bary.evt RXJ0806.3+1527_phase.evt exp=@dmtcalc.expr clob+ unix% dmlist RXJ0806.3+1527_phase.evt cols ... 17 GPHASE Real8 -Inf:+Inf User defined column 18 PHASE Real8 -Inf:+Inf User defined column ...
The expression for phase can be modified to include a frequency shift, fdot, by replacing the GPHASE expression with a more general form:
GPHASE=(time-(.tzero))*(.freq +(.fdot * (time-(.tzero))))"
and then dmextract is used to compute the phase binned lightcurve
unix% dmextract "RXJ0806.3+1527_phase.evt[bin phase=0.0:1.0:0.1]" phase.lc op=generic clob+
This example bins the events into 10 phase bins (going from 0 to 1.0 in 0.1 (10%) period bins). Note that the data are not being binned on time, they are being binned on phase, so the output format must be set as op=generic, not op=ltc1. The phase binned lightcurve is shown in Figure 7.
[Version: full-size]
Figure 7: dmextract: 321.85 seconds period folded lightcurve
While the shape of the curve is the same as the inset in Figure 5, we see that the COUNT_RATE is an order of magnitude lower. This is because the exposure time dmextract uses is the full length of the observation, instead of the per-phase-bin time. dmextract does not know that binning on PHASE implies a binning on TIME, nor does it know how to convert between the two.
For things to work out correctly, any expression involving phase needs to be expressed differently; it must be expressed as a function of time. This will be done by creating Good Time Interval (GTIs).
Create Good Time Intervals (GTIs)
The correct way to filter/bin on phase is to convert PHASE into TIME. This will be accomplished using dmgti. The input to dmgti is a file with a TIME column. For phase binning the time resolution should be smaller than the time resolution of the event file. The aspect solution has a time resolution of 0.256 seconds which meets this need.
First, the PHASE for the times in the barycenter corrected aspect file are computed using dmtcalc as was done earlier with the event file.
unix% dmtcalc RXJ0806.3+1527_bary.asol RXJ0806.3+1527_phase.asol exp=@dmtcalc.expr clob+
The next step is to determine the good times. Suppose for this example the goal was to anlyze the data when the source is brightest. Looking at Figure 5 this corresponds to a PHASE between 0.4 to 0.9. These values are then used with dmgti to determine what times (good times) have 0.4 <= PHASE < 0.9
unix% punlearn dmgti unix% pset dmgti infile=RXJ0806.3+1527_phase.asol unix% pset dmgti outfile=high.gti unix% pset dmgti userlimit="((0.4<=phase)&&(phase<0.9))" unix% pset dmgti clobber=yes unix% dmgti Input MTL file (RXJ0806.3+1527_phase.asol): Output GTI file (high.gti): User defined limit string (((0.4<=phase)&&(phase<0.9))):
dmgti uses the same mathematical expression syntax as dmtcalc.
Examining the output there are 67 time intervals where the PHASE is within the specified range.
unix% dmlist high.gti blocks -------------------------------------------------------------------------------- Dataset: high.gti -------------------------------------------------------------------------------- Block Name Type Dimensions -------------------------------------------------------------------------------- Block 1: PRIMARY Null Block 2: FILTER Table 1 cols x 0 rows Block 3: GTI Table 2 cols x 67 rows
The times in the high.gti file now correspond to times when the phase is between 0.4 and 0.9.
Align GTIs on ACIS exposure frame time boundaries
The next step is to align those GTI time boundaries (start and stop times) to match with the start and stop times of the actual ACIS exposure times. (This step is not necessary when working with HRC data.) This is a small correction to each GTI, just a fraction of the frame time; however, with many GTIs the cumulative adjustment can become significant.
The gti_align script is used to adjust the arbitrary phase times. The script uses the exposure statistics file, stat1.fits, which is why the barycenter time corrected were also applied to it.
unix% pset gti_align times=high.gti unix% pset gti_align statfile=RXJ0806.3+1527_bary.expstats unix% pset gti_align outfile=high_align.gti unix% pset gti_align evtfile=RXJ0806.3+1527_bary.evt unix% gti_align mode=h clobber=yes
Figure 8 shows a comparison of the GTI START and STOP times after the aligmnent. Since the frame time is approximately 0.44 seconds, the amount each GTI is shifted in either side is no more than 0.44 seconds.
[Version: full-size]
Figure 8: GTI Alignment Corrections
The difference in exposure time can be seen when applying two different GTIs
unix% dmkeypar "RXJ0806.3+1527_bary.evt[@high.gti]" ontime echo+ 10601.9740087092 unix% dmkeypar "RXJ0806.3+1527_bary.evt[@high_align.gti]" ontime echo+ 10633.0197870433
A thirty second difference in ONTIME is probably not significant for this observation. The effect, however, is cumulative. For other periodic sources with shorter periods, there will be more GTIs and the correction can become a significant fraction of the total exposure time.
Extract Phase Resolved Spectrum
First, the event file is filtered with the GTI file. This is necessary to correctly set the EXPOSURE (and releated) keywords.
unix% dmcopy RXJ0806.3+1527_bary.evt"[@high_align.gti]" RXJ0806.3+1527_high.evt clob+
Users working with gratings data should follow the work around in the dmcopy opt=all bug when working with TIME/GTI filters.
Then specextract is used to extract the spectrum.
unix% specextract "RXJ0806.3+1527_high.evt[sky=region(ciao.reg)]" RXJ0806.3+1527_high clob+ Running: specextract Version: date Using event file RXJ0806.3+1527_high.evt[sky=region(ciao.reg)] Aspect solution file RXJ0806.3+1527_bary.asol found. Bad-pixel file acisf02795_repro_bpix1.fits found. Mask file acisf02795_000N003_msk1.fits found. Setting bad pixel file Extracting src spectra Creating src ARF Creating src RMF Using mkacisrmf... Grouping src spectrum Updating header of RXJ0806.3+1527_high.pi with RESPFILE and ANCRFILE keywords. Updating header of RXJ0806.3+1527_high_grp.pi with RESPFILE and ANCRFILE keywords.
Even though the data have already been spatially filtered, specextract checks the infile for a spatial filter. Otherwise, applying the same filter twice has no effect.
As before with the lightcurve, background has been ignored when extracting the spectrum.
The spectrum can then be loaded into sherpa for modeling and fitting. The phase resolved spectrum is shown in Figure 10.
[Version: full-size]
Figure 10: Spectrum for Phase 0.4:0.9 with 321.85 second period.
Since the PHASE was selected via GTIs, the EXPOSURE time has automatically been recomputed. The spectrum can now be modeled and fit.
Summary
This thread analyzed the variable, periodic source RXJ0806.3+1527. The data products for the observation were barycenter corrected (axbary). Several timing tools were used to create lightcurves (dmextract and glvary) , search for variability (apowerspectrum), and determine the period of the variation (pfold). The period of 321.85 seconds agrees well with what is reported in the literature for this source in this observation. The data products were then phase folded using dmtcalc and a set of Good Time Intervals (GTIs) were created for the brightest part of the variation using dmgti. These GTIs were aligned to ACIS exposure frame time boundaries using the gti_align script. Finally these GTIs were used to filter the event file and specextract was used to extract the spectrum.
History
03 Jan 2005 | reviewed for CIAO 3.2: no changes |
19 Dec 2005 | updated for CIAO 3.3: default value of dmextract error and bkgerror parameters is "gaussian"; output filenames include ObsID |
01 Dec 2006 | updated for CIAO 3.4: ChIPS and Sherpa versions |
22 Jan 2008 | updated for CIAO 4.0: updated ChIPS and Sherpa syntax; kernel parameter removed from dmtcalc |
25 Jun 2008 | updated image display to place figures inline with text |
09 Jan 2009 | updated for CIAO 4.1: provided both Python and S-lang syntax for ChIPS and Sherpa |
08 Feb 2010 | updated for CIAO 4.2: ChIPS and Sherpa version |
19 Jul 2010 | the S-Lang syntax has been removed from this thread as it is not supported in CIAO 4.2 Sherpa v2. |
13 Jan 2011 | reviewed for CIAO 4.3: no changes |
05 Mar 2012 | reviewed for CIAO 4.4: frequency has been updated to 19.7988900 Hz (instead of 19.794885 Hz) |
24 Apr 2012 | added punlearn step to extracting lightcurve |
03 Dec 2012 | Review for CIAO 4.5. Added note about GTIs and link to pfold. |
10 Dec 2013 | Review for CIAO 4.6; added note about early data. |
22 Dec 2014 | Reviewed for CIAO 4.7; no changes. |
13 Apr 2015 | This thread has been completely rewritten. It now uses an approach that created good time intervals (GTIs) to filter the data and provides some guidance on how to determine the target period. This makes use of the new gti_align script. The bary center correction has been extended to include the aspect solution and the exposure statistics files. |
11 Dec 2015 | Updated for CIAO 4.8; removed work around in the axbary section. |
23 Mar 2016 | Added a warning about REGION block with gratings data. |
05 Apr 2019 | Updated to use matplotlib for plotting. |
30 Nov 2020 | Fix typo in T0 value used in plots. |
07 2022 | Review for CIAO 4.14. Updated for Repro5 and CALDB 4.9.6. |
13 Feb 2023 | Added a note about change to default dither parameters for ACIS. |