Phase-binning a Spectrum

Get Started

Download the sample data: 133 (ACIS-I, PSR B0540-69)

unix% download_chandra_obsid 133 evt2,eph1

For PSR B0540-69, we use the frequency of 19.7988900 Hz from Kaaret et al. (2001, ApJ, 546, 1159).

It is important to note that you should not fold a spectrum on a period that is close to or smaller than the temporal resolution of the data. The majority of ACIS observations are done in timed exposure (TE) mode, collecting photons and reading them out periodically. Each readout is called a "frame" and every event in a given frame has the same time assigned to it. The default resolution is 3.2 s, plus 0.04104 seconds frame transfer time; the exact value used, including the frame transfer time, is stored in the TIMEDEL header keyword:

unix% dmkeypar acis_timed_evt2.fits TIMEDEL echo+
3.24104

Trying to analyze this TE-mode observation on the 50 ms timescale used in this thread would result in meaningless output.

The benefit of continuous clocking (CC) mode observations, such as the one used here, is that they have a much shorter time resolution since the data is continually read out:

unix% dmkeypar acis_cc_evt2.fits TIMEDEL echo+
0.00285

Apply a Barycenter Correction

We begin by applying a barycenter correction 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:

unix% axbary infile=acisf00133N003_evt2.fits orbitfile=orbitf051580864N002_eph1.fits \
      outfile=acis_133_bary_evt2.fits ra=85.04689 dec=-69.33194

This command is explained in detail in the Apply Barycenter Correction thread.

Define a Source Region

In subsequent steps of this thread, we are only interested in results from the source region, not the full field of view. Here we show how to define the source region.

Load the file into ds9:

unix% ds9 acis_133_bary_evt2.fits &

Since this is a CC-mode observation, there is only one spatial dimension (see the Continuous Clocking Mode why topic for details), so the image looks like Figure 1.

Figure 1: Image of CC-mode observation

[Version: full-size]

Figure 1: Image of CC-mode observation

In CC-mode observations, there is only one spatial dimension. The other is sacrificed to increase the temporal resolution of the data.

Outline the source on the image. Figure 2 shows the rotated box used in this example; for information on how to rotate shapes, see this FAQ.

Figure 2: Source region overlaid on the event data

[Version: full-size]

Figure 2: Source region overlaid on the event data

Due to the one-dimensional nature of CC-mode data, the source region is a box instead of a circle.

Then save the region from ds9:

1. Region → Save Regions... → Save As "source.reg".
2. 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 will look something like this:

unix%  more source.reg 
# Region file format: CIAO version 1.0
rotbox(4095.4248,4055.2905,1.3894415,3.0144702,246.69322)

This region will be used in the Run dmextract step.

Create PHASE column

When working with large files, you can save some time by combining the two dmtcalc steps that are used in this section. Here we give the short-and-sweet example; for a full explanation of the steps, refer to the Convert Time to Phase step.

unix% cat dmtcalc.expr 
GPHASE=(time-TSTART)*19.7988900
PHASE=GPHASE-(long)GPHASE

unix% punlearn dmtcalc
unix% dmtcalc infile=acis_133_bary_evt2.fits outfile=acis_133_PHASE.fits \
      expression=@dmtcalc.expr

Note that this step does not take in account the frequency derivative; see the end of the Convert Time to Phase section for information on how to do so.

After completing this step, you can go right to the Extract the Spectrum section of this thread.

unix% punlearn dmtcalc
unix% pset dmtcalc infile=acis_133_bary_evt2.fits
unix% pset dmtcalc outfile=acis_133_GPHASE.fits
unix% pset dmtcalc expression="GPHASE=(time-TSTART)*19.7988900"
unix% dmtcalc 
Input file (acis_133_bary_evt2.fits): 
Output file (acis_133_GPHASE.fits): 
expression(s) to evaluate (GPHASE=(time-TSTART)*19.7988900): 

unix% pset dmtcalc expression="GPHASE=(time-TSTART)*(f0 +(fdot * (time-TSTART)))"

unix% pset dmtcalc expression="GPHASE=(time-TSTART)/(P0 +(Pdot * (time-TSTART)))"

unix% punlearn dmtcalc
unix% pset dmtcalc infile=acis_133_GPHASE.fits
unix% pset dmtcalc outfile=acis_133_PHASE.fits
unix% pset dmtcalc expression="PHASE=GPHASE-(long)GPHASE"
unix% dmtcalc
Input file (acis_133_GPHASE.fits): 
Output file (acis_133_PHASE.fits): 
expression(s) to evaluate (PHASE=GPHASE-(long)GPHASE):

Make a Lightcurve (optional)

At this point, it is possible to make a lightcurve of the phase-folded data. Use dmextract to bin on the PHASE column while simultaneously filtering on the source region:

unix% punlearn dmextract
unix% dmextract "acis_133_PHASE.fits[sky=region(source.reg)][bin PHASE=0:1:0.1]" \
      acis_133_lightcurve.fits opt=generic

The resultant lightcurve can be displayed in ChIPS:

unix% chips
-----------------------------------------
Welcome to ChIPS: CXC's Plotting Package
-----------------------------------------
CIAO 4.4 ChIPS version 1 Friday, December 2, 2011

chips> make_figure("acis_133_lightcurve.fits[cols phase,counts,stat_err]")
chips> set_curve("err.style=cap line.style=none")

Figure 3 shows the plot, which clearly illustrates the counts as a function of phase.

Figure 3: Phase-binned lightcurve

[Version: full-size, PNG]

Figure 3: Phase-binned lightcurve

The counts as a function of phase are plotted with error bars in ChIPS.

Quit ChIPS before continuing:

chips> quit

Important: make sure that any interesting feature that may appear in your lightcurve is not the result of a false period; for more information, see the Timing Analysis with Lightcurves why topic. An unfortunate choice in the width of the phase bin may lead to features that are not real. This is due in part to the fact that events need not be distributed uniformly in the phase histogram even if they are distributed uniformly in phase from the source.

Also be aware the Good Time Intervals are not taken into account when binning on PHASE. If there are gaps in the data due to telemetry saturation then the phase light curve extracted this way will be incorrect.

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

Extract the Spectrum

1. Run dmextract

Now a source spectrum for any phase interval, in this case between 0.2 and 0.4, can be created.

Using the region we defined previously:

unix% punlearn dmextract
unix% pset dmextract \
      infile="acis_133_PHASE.fits[sky=region(source.reg),phase=0.2:0.4][bin pi]"
unix% pset dmextract outfile=acis_133_phase_bin_spectrum.fits
unix% dmextract 
Input event file (acis_133_PHASE.fits[sky=region(source.reg),phase=0.2:0.4][bin pi]): 
Enter output file name (acis_133_phase_bin_spectrum.fits): 

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

---

2. Plot the Spectrum

The resulting spectrum may be plotted in Sherpa:

unix% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.4 Sherpa version 1 Friday, December 2, 2011

sherpa> load_pha("acis_133_phase_bin_spectrum.fits")
sherpa> set_analysis("channel", type="counts")
sherpa> plot_data()
sherpa> set_curve("symbol.style=none line.style=solid")

which will create the plot shown Figure 4. The value along the y-axis is the sum of all the counts that occurred in the chosen 20% of the observation's cycle. In this example, we did not include "f dot" in the Convert Time to Phase step.

Figure 4: Phase-resolved spectrum

[Version: full-size, PNG]

Figure 4: Phase-resolved spectrum

This plot shows the spectrum of the source when it has a phase between 0.2 and 0.4.

Quit Sherpa when you are done:

sherpa> quit

Caveats

1. This thread follows the usual guidelines: filtering on one column doesn't update any information in the other columns. This means that it does not correctly account for the data in the GTI tables when binning on phase. That is, the EXPOSURE keyword in the header of the new file is not the sum of the new time bins.