Simulating two overlapping sources

In each run of marx the user can choose one and only one source from the list in The spectrum and the spatial shape of the X-ray source. Most of these sources are simple geometric shapes like a point or a disk. The IMAGE Source allows the user to specify an itensity image as input to define a complicated morphology on the sky, but the spectrum is the same in every position.

If spectrum and intensity change with position, the user can write his/her own C-code to generate the photons as described in USER Source. Some of those codes are available for download, see The spectrum and the spatial shape of the X-ray source for links.

However, in most cases it is easier to run marx several times and use marxcat to combine the simulations.

In this example we want to simulate a binary star, where both components are separated by only a few arcseconds. In the center of the field-of-view that is enough to separate the images of component A and B, but what happens when we insert a grating? Will the dispersed spectra overlap, or can they be separated in the data?

The parameters in the following example are inspired by the observations of the M-dwarf binary EQ Peg (ObsID 8453). An analysis is presented by Liefke et al. (2008). We will run two marx simulations, one for EQ Peg A and one for EQ Peg B. The parameters of the observations (poining coordinates, roll angle, exposure time, grating) have to be the same for both simulations, but the parameters for the sources (coordinates, count rate, spectrum) are different.

Creating a file with the input spectra

First, we need to generate the input spectra. For this example, we use Sherpa, but XSPEC or ISIS work very similar. We need two 2-column text files that tabulate the model flux [photons/sec/keV/cm^2] (second column) as a function of energy [keV] (first column).

Stars are generally well described by an optically thin, collisionally excited plasma. We use a model with only a single temperature component with a different temperature for EQ Peg A and EQ Peg B and slightly non-solar abundances as suggested by Liefke et al. (2008).

# set source properties
set_source(xsvapec.a1)
a1.Ne = 1.2
a1.Fe = 0.8

# get source
my_src = get_source()

# set energy grid
bin_width = 0.01
energies = np.arange(0.15, 12., bin_width)

### EQ Peg A ###
# check that out! Download and take number from real data.
a1.norm = 0.001
a1.kT = 0.95
# evaluate source on energy grid
flux = my_src(energies)

# Sherpa uses the convention that an energy array holds the LOWER end of the bins;
# Marx that it holds the UPPER end of a bin.
# Thus, we need to shift energies and flux by one bin.
# Also, we need to divide the flux by the bin width to obtain the flux density.
save_arrays("EQPegA_flux.tbl", [energies[1:], flux[:-1] / bin_width], ["keV","photons/s/cm**2/keV"], ascii=True)

### EQ Peg B ###
a1.norm = 0.0003
a1.kT = 0.7
flux = my_src(energies)

save_arrays("EQPegB_flux.tbl", [energies[1:], flux[:-1] / bin_width], ["keV","photons/s/cm**2/keV"], ascii=True)

More details about the format of the marx input spectrum can be found at SpectrumFile.

Running marx

We run marx twice - once for each star.

marx SourceFlux=-1 SpectrumType="FILE" SpectrumFile="EQPegA_flux.tbl" \
       ExposureTime=5e5 TStart=2006.9 \
       OutputDir=EQPegA GratingType="HETG" DetectorType="ACIS-S" \
       DitherModel="INTERNAL" RA_Nom=352.966 Dec_Nom=19.935 Roll_Nom=290 \
       SourceRA=352.96851 SourceDEC=19.937133

marx SourceFlux=-1 SpectrumType="FILE" SpectrumFile="EQPegB_flux.tbl" \
       ExposureTime=5e5 TStart=2006.9 \
       OutputDir=EQPegB GratingType="HETG" DetectorType="ACIS-S" \
       DitherModel="INTERNAL" RA_Nom=352.966 Dec_Nom=19.935 Roll_Nom=290 \
       SourceRA=352.97011 SourceDEC=19.937225

The SourceFlux parameter may be used to indicate the integrated flux of the spectrum. The value of -1 means that the integrated flux is taken from the file. The results of the simulation will be written to sub-directories called EQPegA and EQPegB, as specified by the OutputDir parameter. We use marxcat to combine both simulations in the directory EQPeg_both:

marxcat EQPegA EQPegB EQPeg_both

In the simulation we know exactly which photon comes from which star, so we can generate three fits files, one for EQ Peg A only, one for EQ PEg B only, and one that contains both sources, similar to the observed data.

marx2fits EQPegA EQPegA.fits
marx2fits EQPegB EQPegB.fits
marx2fits EQPeg_both EQPeg_both.fits

The fits files can be further processed with standard CIAO tools. As some of these tools require the aspect history, the marxasp program is used to create an aspect solution file. Since both simulations used the same pointing and dither, we can use the same asol file for all three fits files.

marxasp MarxDir="EQPegA" OutputFile="EQPeg_asol1.fits"

Analyzing the results

It is interesting to look at the event file with a viewer such as ds9. The grating spectra of EQ Peg A and EQ Peg B are close to each other on the detector, and we’ll now test if the spectrum of EQ Peg A contaminates the spectrum of EQ Peg B.

Two spectra on a detector very close to each other.

Small section of the detector image of the combined simulation

The grating spectra of both sources are located fairly close to each other. On the left is an MEG arm, the fainter signal of an HEG arm is seen on the right. The green lines mark the extraction regions for the EQ Peg B MEG and HEG spectra. tgextract subdivides this green region into source and background regions for spectral extraction.

We use CIAO to extract the grating spectrum from EQPegB.fits and EGPeg_both.fits.

#!/bin/sh

X=4068.3;    # determined by hand in ds9
Y=4112.7;    # determined by hand in ds9

# We'll use the same mask for both extractions, because we want the same region.
tg_create_mask infile="EQPegB.fits" outfile="EQPegB_reg1a.fits" \
  use_user_pars=yes last_source_toread=A \
  sA_id=1 sA_zero_x=$X sA_zero_y=$Y clobber=yes grating_obs=header_value \
  width_factor_hetg=10


### Extract EQ Peg B from the EQ PEg B - only simulation ###

tg_resolve_events infile="EQPegB.fits" outfile="EQPegB_evt1a.fits" \
  regionfile="EQPegB_reg1a.fits" acaofffile="EQPeg_asol1.fits" \
  osipfile=CALDB verbose=0 clobber=yes

tgextract infile="EQPegB_evt1a.fits" outfile="EQPegB_pha2.fits" mode=h clobber=yes


### Extract EQ Peg B from the combined simulation ###

tg_resolve_events infile="EQPeg_both.fits" outfile="EQPegB_both_evt1a.fits" \
  regionfile="EQPegB_reg1a.fits" acaofffile="EQPeg_asol1.fits" \
  osipfile=CALDB verbose=0 clobber=yes

tgextract infile="EQPegB_both_evt1a.fits" outfile="EQPegB_both_pha2.fits" mode=h clobber=yes


# make response matrices for the MEG +1 and -1 order
# We are extracting the exact same region in both cases above for EXACTLY the same
# observational parameters (aimpoint, exposure time, ...).
# Thus the arf and rmf files will be identical, so we have to run this only once.
mkgrmf grating_arm="MEG" order=1 outfile="eqpeg_meg1_rmf.fits" srcid=1 detsubsys=ACIS-S3 \
   threshold=1e-06 obsfile="EQPegB_pha2.fits" regionfile="EQPegB_pha2.fits" \
   wvgrid_arf="compute" wvgrid_chan="compute"

mkgrmf grating_arm="MEG" order=-1 outfile="eqpegmeg-1_rmf.fits" srcid=1 detsubsys=ACIS-S3 \
   threshold=1e-06 obsfile="EQPegB_pha2.fits" regionfile="EQPegB_pha2.fits" \
   wvgrid_arf="compute" wvgrid_chan="compute"

fullgarf phafile="EQPegB_pha2.fits"  pharow=3 evtfile="EQPegB.fits" \
   asol="EQPeg_asol1.fits" engrid="grid(eqpeg_meg1_rmf.fits[cols ENERG_LO,ENERG_HI])" \
   maskfile=NONE dafile=NONE dtffile=NONE badpix=NONE rootname="eqpeg" ardlibqual=";UNIFORM;bpmask=0"

fullgarf phafile="EQPegB_pha2.fits"  pharow=4 evtfile="EQPegB.fits" \
   asol="EQPeg_asol1.fits" engrid="grid(eqpeg_meg1_rmf.fits[cols ENERG_LO,ENERG_HI])" \
   maskfile=NONE dafile=NONE dtffile=NONE badpix=NONE rootname="eqpeg" ardlibqual=";UNIFORM;bpmask=0"

Now, we want to make use of our marx simulation to see how much the spectrum of EQ Peg A contaminates the extracted grating spectrum of EQ Peg B. We compare the spectra in Sherpa. The difference between the two spectra in some of the lines is caused by photons from EQ Peg A that fall in the extraction region of EQ Peg B.

The spetra differ only at short wavelength.

Contaminated and clean spectrum extracted for EQ Peg B (MEG order +1).

The red spectrum shows the uncontaminated spectrum of EQ Peg B, extracted from the simulation that contained only one source. The black spectrum is extracted from the combined fits file that contains both EQ Peg A and EQ Peg B.

Looking at this figure, we see that both spectra are very similar, but not identical - that is the contamination by EQ Peg A.

tgextract uses only the innermost part of the green region in figure Small section of the detector image of the combined simulation for the source spectrum and apparently that is narrow enough to avoid most of the contamination. Still, if we want to measure or fit lines in the EQ Peg B spectrum we have to be careful.

marx simulations do not contain any background, but real data does and one common (although arguably not the best) approach is to subtract a background spectrum from the data. What would happen if we blindly subtracted the background from the EQ Peg B spectrum?

spectrum with negative bins

“Background subtracted” spectrum for EQ Peg B

Because the bright spectrum of EQ Peg A falls in the detector region where the background is measured, the background is overestimated and causes bins with negative fluxes here.

Figure “Background subtracted” spectrum for EQ Peg B shows that we end up with a spectrum with a significant number of negative bins because the much brighter spectrum of EQ Peg A falls in the region where tgextract extracts the background, leading to an overestimate of the background.

Summary

This example shows how marxcat can be used to build up a more complex distribution of sources on the sky from simple components. We can asses how much the sources overlap and compare the simulated spectra of each component and thus estimate the contamination in the real data.

In this example we simulated grating spectra of a binary star, but the same principle works for all cases where sources overlap on the detector, e.g. a cluster of galaxies with a foreground HMXB, binary AGN, or imaging of a dense star forming region.