Marx 5.0.0

Simulating a user-defined spectrum with Marx

The purpose of this example is to show how to use marx to simulate an HETG grating spectrum of a star whose spectrum is represented by an aped model.

Creating the spectral file for marx

The generation of the spectral file for marx is a bit more involved than that of the powerlaw example. As before, marxflux will be used to generate the spectral file from the parameter file describing the model. But in this case, the plasma model must be created before isis can compute the model flux. So, the idea is to put the definition of the model into a file called apedfun.sl, which not only defines the model but also creates the parameter file that marxflux will load. The contents of apedfun.sl are:

    plasma(aped);
    define create_aped_model ()
    {
       variable p = default_plasma_state ();
       p.temperature = [9.7e7, 3.2e7];
       p.norm = [ 0.016, 0.0156];
       p.elem = [Fe, Ne];
       p.elem_abund = [ 0.2, 2 ];

       create_aped_fun ("Aped", p);
       fit_fun("Aped(1)");
    }

    create_aped_model ();
    save_par ("aped.p");

Then the marxflux script may be used to create the spectral file via

  marxflux --script apedfun.sl -l '[1:30:#16384]' aped.p apedflux.tbl

Here, the --script apedfun.sl arguments instruct marxflux to execute the apedfun.sl file. The -l '[1:30:#16384]' parameters indicate that a 16384 point wavelength grid running from 1 to 30 angstroms is to be used for the spectrum--- marxflux will convert this to energy units. A plot of the spectrum in apedflux.tbl is shown below.

Running marx

The next step is to run marx in the desired configuration. For this example, ACIS-S and the HETG are used:

    marx SourceFlux=-1 SpectrumType="FILE" SpectrumFile="apedflux.tbl" \
       ExposureTime=80000 OutputDir=aped DitherModel="INTERNAL" \
       GratingType="HETG" DetectorType="ACIS-S"

This will run the simulation and place the results in a subdirectory called aped. The results may be converted to a standard Chandra level-2 fits file by the marx2fits program:

    marx2fits aped aped_evt2.fits

The resulting fits file (aped_evt2.fits) may be further processed with standard CIAO tools, as described below. As some of these tools require the aspect history, the marxasp program will be used to create an aspect solution file that matches the simulation:

    marxasp MarxDir="aped" OutputFile="aped_asol1.fits"

Creating a type-II PHA file

For a Chandra grating observation, the CIAO tgextract tool may used to create a type-II PHA file. Before this can be done, an extraction region mask file must be created using the CIAO tg_create_mask, followed by order resolution using tg_resolve_events. The first step is to determine the source position, which is used by tg_create_mask. While there are many ways to do this, perhaps the most flexible way is to use findzo as described on its web page. For this particular example, the position of the zeroth order in sky coordinates is (4017.88, 4141.43).

The following Bourne shell script runs the above set of tools to create a PHA2 file called aped_pha2.fits.

    #!/bin/sh
    evtfile="aped_evt2.fits"
    evt1afile="aped_evt1a.fits"
    reg1afile="aped_reg1a.fits"
    asol1file="aped_asol1.fits"
    pha2file="aped_pha2.fits"

    X=4017.88;    # determined by findzo
    Y=4141.43;    # determined by findzo

    tg_create_mask infile="$evtfile" outfile="$reg1afile" \
      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

    tg_resolve_events infile="$evtfile" outfile="$evt1afile" \
      regionfile="$reg1afile" acaofffile="$asol1file" \
      alignmentfile="$asol1file" osipfile=CALDB verbose=0 clobber=yes

    tgextract infile="$evt1afile" outfile="$pha2file" mode=h clobber=yes

The resulting file may be loaded into isis and displayed using

   isis> load_data ("aped_pha2.fits");
   isis> list_data;

with the result

Current Spectrum List:
 id    instrument   m prt src    use/nbins   A   R     totcts   exp(ksec)  target
  1     HETG-ACIS  -3  1   1    8192/ 8192   -   -  3.5100e+02    78.988  POINT
  2     HETG-ACIS  -2  1   1    8192/ 8192   -   -  1.0040e+03    78.988  POINT
  3     HETG-ACIS  -1  1   1    8192/ 8192   -   -  1.4086e+04    78.988  POINT
  4     HETG-ACIS   1  1   1    8192/ 8192   -   -  1.2569e+04    78.988  POINT
  5     HETG-ACIS   2  1   1    8192/ 8192   -   -  6.7800e+02    78.988  POINT
  6     HETG-ACIS   3  1   1    8192/ 8192   -   -  2.2800e+02    78.988  POINT
  7     HETG-ACIS  -3  2   1    8192/ 8192   -   -  2.6500e+03    78.988  POINT
  8     HETG-ACIS  -2  2   1    8192/ 8192   -   -  6.9800e+02    78.988  POINT
  9     HETG-ACIS  -1  2   1    8192/ 8192   -   -  2.9684e+04    78.988  POINT
 10     HETG-ACIS   1  2   1    8192/ 8192   -   -  2.9691e+04    78.988  POINT
 11     HETG-ACIS   2  2   1    8192/ 8192   -   -  5.8600e+02    78.988  POINT
 12     HETG-ACIS   3  2   1    8192/ 8192   -   -  2.0430e+03    78.988  POINT

A plot of the MEG-1 spectrum, which corresponds to id=9 in the above list may be created using

   isis> xrange (1, 25);
   isis> group_data (9, 2);
   isis> plot_data (9);

In making the above plot, the group_data function was used to rebin the data by combining adjacent wavelength channels.

This page was last updated Jan 25, 2012 by John E. Davis.
Technical questions should be addressed to marx-help at space mit edu.