Last modified: 09 Dec 2024

URL: https://cxc.cfa.harvard.edu/ciao/threads/marx_sim/

Using MARX to Simulate an Existing Observation

CIAO 4.17 Science Threads


Overview

Synopsis:

Simulating rays with MARX is very similar to the steps needed to reproject ChaRT rays, and is only a matter of setting a couple additional parameters. The MARX software is used to generate and project the rays onto the detector plane.

Note that MARX was originally intended as a PSF simulator and not an observation simulator.

Using MARX allows the user to take into account all changes to the photon distribution emerging from the HRMA due to the detector response. In particular, the detector QE & QEU and the roll are accounted for. In addition to simulating the detector response, MARX uses the ray weights to account for the mirror effects, i.e. different efficiency of different shells at different angles/energies.

Related Links:

Last Update: 09 Dec 2024 - Added note about MARX installation; happens by default with ciao-install in CIAO 4.17 or can be installed with conda.


Contents


Get Started

This example uses the point source spectrum of τ Canis Majoris—an eclipsing spectroscopic binary, which is extracted near the detector's aimpoint for an observation of the open cluster, NGC 2362—and used to simulate a point spread function for this particular observation.

Download the sample data: 4469 (ACIS-I, NGC 2362)

Re-process it with chandra_repro:

unix% cd 4469
unix% chandra_repro . repro/ check_vf_pha=yes

The files we are going to use are in the repro/ directory, shown in bold (note that only the events file needs to be explicitly specified since the remaining files can be found automatically by specextract).

unix% cd repro/
unix% ls

In order to simulate a PSF of the point source, a distribution of incoming photons is needed to interact with the telescope's optics before being projected on to a detector. Generally, when simulating a PSF for an observation, the incident source spectrum should be background subtracted. τ Canis Majoris' coordinates are RA=07:18:42.457 [109.6769 deg], Dec=-24:57:15.78 [-24.954383 deg], and specextract is used to extract the source spectrum in a circular region; a background with an annulus; and calculate unweighted response functions. Since the background will not be modeled, responses for the background need not be made.

unix% specextract \
  infile=acisf04469_repro_evt2.fits"[sky=circle(07:18:42.457,-24:57:15.78,0.0666667')]" \
  outroot=4469 \
  bkgfile=acisf04469_repro_evt2.fits"[sky=annulus(7:18:42.457,-24:57:15.78,0.0666667',0.15')]" \
  bkgresp=no weight=no correctpsf=yes grouptype=NONE verbose=0 mode=h

unix% ls spec
4469.arf  4469_bkg.pi  4469.corr.arf  4469.pi  4469.rmf 

MARX requires the spectrum to be in units of flux density, photon/s/cm2/keV, input as an ASCII table with the first column as energies in keV and flux densities in the second column. The observed spectrum is mapped from channel- to energy-space by the RMF, with the units as counts instead of photons. To analytically obtain the "true" source spectrum—which will be in photon units—requires deconvolving the responses from the observed spectrum. Numerically, this can be done by dividing the observed spectrum by the detector responses. In Sherpa, the plot_ratio command—which divides the observed spectrum with a model spectrum, the model spectrum being a source spectrum convolved with the ARF and RMF—using a constant source model set to unity, so that the model spectrum is a discrete representation of just the detector responses, may be used for this purpose.

sherpa> load_data("4469.pi")
sherpa> subtract()
sherpa> notice("0.3:11.0")

sherpa> set_source(polynom1d.truespec)
sherpa> truespec.c0 = 1

sherpa> photonflux = get_ratio_plot().y # ph/s/cm**2/keV
sherpa> energy = get_ratio_plot().x # keV

sherpa> photonflux[photonflux < 0] = 0 # replace negative fluxes, from background subtraction, with zeroes

sherpa> save_arrays("tCanMaj.dat",[energy,photonflux],["keV","ph_s_cm2_keV"],ascii=True)

The save_arrays command writes the energy and photonflux arrays between ~0.3-11 keV, rounded to the closest energy bin, to the ASCII file tCanMaj.dat.


Setting up MARX

This thread uses MARX v5.5.3; if you do not already have the software, download and install MARX before continuing. [New] MARX is installed automatically when CIAO is installed using the ciao-install script or it can be installed using conda by adding marx to the list of packages to install. Alternatively, the install_marx script may also be used to download and install the latest available version of MARX.

Set up your environment to access all the files and directories which MARX needs:

---(t)csh---
unix% setenv MARX_ROOT <local_marx_directory>/marx-5.5.3
unix% setenv MARX_DATA_DIR ${MARX_ROOT}/share/marx/data
unix% setenv PATH ${PATH}:${MARX_ROOT}/bin
unix% setenv PFILES ${PFILES}:${MARX_ROOT}/share/marx/pfiles

---bash---
unix% export MARX_ROOT=<local_marx_directory>/marx-5.5.3
unix% export MARX_DATA_DIR=${MARX_ROOT}/share/marx/data
unix% export PATH=${PATH}:${MARX_ROOT}/bin
unix% export PFILES=${PFILES}:${MARX_ROOT}/share/marx/pfiles

If you installed MARX from the CIAO conda channel then you should set:

---(t)csh---
unix% setenv MARX_ROOT $CONDA_PREFIX

---bash---
unix% export MARX_ROOT=$CONDA_PREFIX

and the rest of the parameters should be set using that value for MARX_ROOT.

Alternatively, the MARX DataDirectory parameter can be defined to point to the data files. By default, this parameter is set to the contents of the $MARX_DATA_DIR environment variable. At this point, your system is ready to run simulations.


Simulating the Source

Before running the simulation, it is necessary to determine the nominal position of the detector during the observation. These values are stored in the header of the level=2 events file:

unix% dmlist acisf04469_repro_evt2.fits header |grep _NOM
0055 RA_NOM                     109.6834360799           Real8        Nominal RA
0056 DEC_NOM                    -24.9551940942           Real8        Nominal Dec
0057 ROLL_NOM                    22.3127528406           Real8        Nominal Roll

Additionally, the source coordinates, in decimal degrees, are needed to produce an appropriate PSF, which is (RA, Dec) = (109.6769 deg, -24.954383 deg) as noted earlier.

Other information needed for the simulation is the observation date, which affects the contamination model applied to the detector, and exposure time.

unix% dmkeypar acisf04469_repro_evt2.fits EXPOSURE echo+
97883.798919246

unix% dmkeypar acisf04469_repro_evt2.fits TSTART echo+
188535838.52712
[NOTE]
TSTART Date Format

The MARX TSTART parameter can be in terms of "mission time" or decimal years. Values <2100 are interpreted as decimal years, otherwise the value is treated as seconds on the spacecraft clock.

Additionally, the appropriate telescope and detector configuration should be set. First, copy over the parameter file locally:

unix% cp $MARX_ROOT/share/marx/pfiles/marx.par ./marx5sim.par

and set the necessary parameters:

unix% pset ./marx5sim ExposureTime=97883.798919246
unix% pset ./marx5sim TStart=188535838.52712

unix% pset ./marx5sim OutputDir=tCanMaj_marxsim.dir

unix% pset ./marx5sim SourceRA=109.6769
unix% pset ./marx5sim SourceDEC=-24.954383

unix% pset ./marx5sim RA_Nom=109.6834360799
unix% pset ./marx5sim Dec_Nom=-24.9551940942
unix% pset ./marx5sim Roll_Nom=22.3127528406

unix% pset ./marx5sim DetectorType=ACIS-I
unix% pset ./marx5sim GratingType=NONE
unix% pset ./marx5sim MirrorType=HRMA

There are a few things to note in the pset commands:

  • Make sure to use the correct detector for the DetectorType parameter; valid options are: HRC-S, ACIS-S, HRC-I, or ACIS-I. Valid options for the GratingType parameter are: HETG, LETG, or NONE. The ACIS array-types are specific, defined sets of CCDs:

    • ACIS-I = ACIS-0123
    • ACIS-S = ACIS-456789
  • The Exposuretime is set to ensure that the simulated output has the desired exposure length in units of seconds. The "Simulation Control: Exposure Time" section of the MARX manual has more information on setting this parameter.

  • If simulating a source directly on-axis, then SourceRA = RA_NOM and SourceDEC = DEC_NOM. If an off-axis simulation is being performed, they should not be the same.

To define the source being simulated, the ASCII spectrum file created earlier will be spatially distributed as a point source:

unix% pset ./marx5sim SpectrumType=FILE
unix% pset ./marx5sim SpectrumFile=tCanMaj.dat
unix% pset ./marx5sim SourceType=POINT
unix% pset ./marx5sim SourceFlux=-1

The SourceFlux parameter is used to set the overall normalization of the input spectrum. Setting SourceFlux to a positive, non-zero value will cause MARX to re-normalize the spectrum read from the ASCII file to the specified total flux; a negative value will cause MARX to normalize the input spectrum by the specified factor. Using SourceFlux=-1 will use the normalization inherent to the input spectrum.

The telescope dithering can also be included—and is necessary for sub-pixel-scale ACIS analysis.

unix% pset ./marx5sim DitherModel=FILE
unix% pset ./marx5sim DitherFile=pcadf04469_000N001_asol1.fits

When simulating a PSF of an observation, the aspect solution from the observation may be used; the standard Chandra Lissajous dither pattern may also be used instead by setting DitherModel=INTERNAL—which will be used when simulating an unobserved source, for example during proposal planning—and the corresponding dither amplitudes to DitherAmp_RA=8 and DitherAmp_Dec=8 for ACIS observations or DitherAmp_RA=20 and DitherAmp_Dec=20 for HRC observations (in arcsec).

[TIP]
Correcting for SIM Offset and Drift in an Observation

Observations are often made with the science instrument module (SIM) is often offset from its nominal position. If these offsets are not included in the MARX run, the simulated PSF will not be placed at the correct location on the detector. The nominal MARX SIM positions are at:

Detector SIM_X [mm] SIM_Y [mm] SIM_Z [mm]
ACIS-I -0.78234819833843 0 -233.5924630914
ACIS-S -0.68426746699586 0 -190.1325231040
HRC-I -1.0402925884 0 126.9854943053
HRC-S (spec) -1.4295857921 0 250.4559758190
HRC-S (image) -1.5333365632 0 250.4559758190

and can be also looked up in the $MARX_DATA_DIR/caldb/telD1999-07-23aimptsN0002.fits file or by running the MARX ray projection with the default offset values of zero.

If there is a SIM offset from the nominal location, first, compare the SIM values in the original Chandra events file to the ACIS-I values in the table.

unix% dmlist acisf04469_repro_evt2.fits header | grep SIM_
0039 SIM_X                   -0.78090834371673           Real8        SIM focus pos (mm)
0040 SIM_Y                                   0           Real8        SIM orthogonal axis pos (mm)
0041 SIM_Z                     -233.5874344608           Real8        SIM translation stage pos (mm)

MARX SIM_X: -0.78234819833843
MARX SIM_Y: 0
MARX SIM_Z: -233.5924630914

Some observations are done with an offset of several millimeters. The difference, calculated as event_file_value - marx_value, is small in this case:

SIM_X: -0.78090834371673 - (-0.78234819833843) =  0.00143985462170

SIM_Z: -233.5874344608   - (-233.5924630914)   = 0.0050286306

The SIM offset generally has a negligible effect, but often manifests itself by placing the source at a different position angle relative to the detector aimpoint. The determined SIM offsets can be then be set:

unix% pset ./marx DetOffsetX=0.001434985462179
unix% pset ./marx DetOffsetY=0
unix% pset ./marx DetOffsetZ=0.0050286306

At the start of the mission, the SIM was located on-axis, but has since drifted due to changes in telescope geometry, with effects that can no longer be neglected. Fortunately, if an aspect solution is used for the simulation's dither model, MARX will account for this drift using the available positional information at all points in time from the aspect solution file.

If MARX's internal dither model is used instead, then the SIM drift can be accounted for, in addition to the SIM offset, in the DetOffset parameters. The drift is characterized in the DY_AVG and DZ_AVG header keywords (in units of millimeters) for an events file that has gone through Repro IV in the archives; has been reprocessed with chandra_repro; or has had r4_header_update applied to it.

unix% dmkeypar acisf00942N004_evt2.fits DY_AVG echo+
0.0
unix% dmkeypar acisf00942N004_evt2.fits DZ_AVG echo+
0.0

To obtain the appropriate DetOffset values to run the MARX projection, add the calculated SIM offset values to the average SIM drift values.

DetOffsetX: N/A + 0.00143985462170 = 0.00143985462170

DetOffsetY: 0  + 0 = 0

DetOffsetZ: 0.0 + (0.0050286306) = 0.0050286306

To avoid pixelization effects in the pseudo-events file generated from the simulation, there is an AspectBlur parameter which can be set. Users may experiment with the values of AspectBlur to better match the observation since it may deviate from the default value depending on the source spectral distribution and instrument used. In this particular example, the figures were generated with an aspect blur of ~0.2 arcsec.

unix% pset ./marx5sim AspectBlur=0.19
[TIP]
Need an Aspect Solution File for the Internal Dither Model?

Should an aspect solution be needed for your analysis, and DitherModel=INTERNAL be used with your MARX simulation, then the marxasp tool can be used.

unix% cp $MARX_ROOT/share/marx/pfiles/marxasp.par .
unix% marxasp MarxDir=tCanMaj_marxsim.dir OutputFile=tCanMaj_asol.fits

This asol file can be used in CIAO along with the generated pseudo-events file. Note that if marxasp is used with a simulation where DitherModelINTERNAL, an error message is produced.

Now run the tool. Note the syntax (@@) required to use the local parameter file; see Example 9 of ahelp parameter for details.

unix% marx @@./marx5sim.par
MARX version 5.5.3, Copyright (C) 2002-2023 Massachusetts Institute of Technology

... screen output omitted...

Initializing source type POINT...
Reading spectrum from:
	tCanMaj.dat
Your spectrum has a total flux of 5.459085e-04 photons/cm^2/sec
System initialized.

Starting simulation.  Exposure Time set to 9.788380e+04 seconds
Collecting 100000 photons...
	61571 collected.
Reflecting from HRMA
Detecting with ACIS-I

Writing output to directory 'tCanMaj_marxsim_test' ...
Total photons: 61571, Total Photons detected: 8382, (efficiency: 0.136136)
  (efficiency this iteration  0.136136)  Total time: 97884.252890

Using simulate_psf to run MARX

simulate_psf may be used instead of explicitly running MARX and setting all the parameters.

unix% pset simulate_psf infile=acisf04469_repro_evt2.fits
unix% pset simulate_psf outroot=tCanMaj_marxsim
unix% pset simulate_psf ra=109.6769
unix% pset simulate_psf dec=-24.954383
unix% pset simulate_psf spectrum=tCanMaj.dat
unix% pset simulate_psf blur=0.19
unix% pset simulate_psf readout_streak=yes
unix% pset simulate_psf pileup=no
unix% pset simulate_psf ideal=no
unix% pset simulate_psf extended=no
unix% simulate_psf 

This will replicate the behavior demonstrated in the steps throughout this thread, including the creation of the projected event file and FITS image of the PSF. The key parameters being the source RA and Dec and the reference event file. Many of the MARX parameters that are calculated are automatically determined and set by the script, taking the information from the input event file's header.

[WARNING]
Input Spectrum for simulate_psf

The input spectrum for simulate_psf has stricter requirements than the spectrum used directly with MARX. The simulate_psf spectrum must only contain energies between 0.3 and 10 keV; the flux at all the energies must be >0.

There are a couple approaches available to handle these restrictions. The most obvious is when starting from the extracted spectrum of an observation, restrict the noticed energy range to the imposed limits, and set any fluxes ≤0 to a very small positive value. For example:

 
sherpa> load_data("4469.pi")
sherpa> subtract()
sherpa> notice("0.3:10.0")

sherpa> set_source(polynom1d.truespec)
sherpa> truespec.c0 = 1

sherpa> photonflux = get_ratio_plot().y # ph/s/cm**2/keV
sherpa> energy = get_ratio_plot().x # keV

sherpa> photonflux[photonflux <= 0] = 1e-12 # replace negative and zero fluxes,
                                            # from background subtraction, with a 
                                            # tiny, positive value

sherpa> save_arrays("tCanMaj.dat",[energy,photonflux],["keV","ph_s_cm2_keV"],ascii=True)

On the other hand, suppose you already have an existing appropriately formatted ASCII input spectrum, from a model simulation for example, and want to work with it instead. Then omitting all the entries with fluxes ≤0 is a valid approach, and then restricting the energy range between 0.3 and 10 keV with various Numpy tools.

sherpa> import pycrates as pcr
sherpa> cr = pcr.read_file("tCanMaj.dat") # read ascii spectrum into a crate

sherpa> energy = cr.get_column("KEV").values # load the energy column to a variable
sherpa> flux = cr.get_column("PH_S_CM2_KEV").values # load the flux column to a variable

# get indicies where the fluxes are > 0 and omit energies <0.3 keV & >10 keV
sherpa> idx = (flux > 0) & (energy > 0.3) & (energy <= 10)
sherpa> flux2 = flux[idx]
sherpa> energy2 = energy[idx]

sherpa> save_arrays(tCanMaj_simulate-psf.dat,[energy2,flux2],["keV","ph_s_cm2_keV"],\
ascii=True,clobber=True) # write filtered spectrum to file 'tCanMaj_simulate-psf.dat'
[NOTE]
Default Parameters

Be aware that the default MARX parameters do not entirely match the default parameters of simulate_psf. In particular, MARX simulates the readout streak and the rays are projected on to a non-ideal detector with QE<1 and pileup affects are not accounted. Furthermore, the chips are not extended to visualize the rays that would otherwise fall onto chip gaps.

simulate_psf in contrast does not simulate the readout streak while extending the visualization into the chip gap and assumes an ideal detector (QE=1) with its default settings.


Create an Events File

MARX creates a number of ASCII files in the specified output directory (tCanMaj_marxsim.dir):

unix% ls tCanMaj_marxsim.dir
b_energy.dat  det_theta.dat  obs.par      sky_roll.dat  xpos.dat    zcos.dat
det_dy.dat    energy.dat     pha.dat      time.dat      ycos.dat    zpos.dat
det_dz.dat    marx.par       sky_dec.dat  xcos.dat      ypixel.dat
detector.dat  mirror.dat     sky_ra.dat   xpixel.dat    ypos.dat

These .dat files need to be converted to a FITS events file before they can be used in CIAO. The MARX tool marx2fits does this, given the directory name and the output filename:

unix% marx2fits --pixadj=EDSER tCanMaj_marxsim.dir tCanMaj.fits
Reading subpix file
Examining tCanMaj_marxsim.dir/time.dat
Examining tCanMaj_marxsim.dir/detector.dat
Examining tCanMaj_marxsim.dir/energy.dat
Examining tCanMaj_marxsim.dir/b_energy.dat
Examining tCanMaj_marxsim.dir/xpos.dat
Examining tCanMaj_marxsim.dir/ypos.dat
Examining tCanMaj_marxsim.dir/zpos.dat
Examining tCanMaj_marxsim.dir/xcos.dat
Examining tCanMaj_marxsim.dir/ycos.dat
Examining tCanMaj_marxsim.dir/zcos.dat
Examining tCanMaj_marxsim.dir/xpixel.dat
Examining tCanMaj_marxsim.dir/ypixel.dat
Examining tCanMaj_marxsim.dir/pha.dat
Examining tCanMaj_marxsim.dir/mirror.dat
Examining tCanMaj_marxsim.dir/sky_ra.dat
Examining tCanMaj_marxsim.dir/sky_dec.dat
Examining tCanMaj_marxsim.dir/sky_roll.dat
Examining tCanMaj_marxsim.dir/det_dy.dat
Examining tCanMaj_marxsim.dir/det_dz.dat
Examining tCanMaj_marxsim.dir/det_theta.dat
unix%

By default, data obtained from the data archive—and if using the default pix_adj values when reprocessing the data with chandra_repro—has the EDSER sub-pixel event repositioning algorithm applied to ACIS data, and no pixel randomization is applied to HRC data. The --pixadj switch is used to replicate this behavior, where --pixadj=EDSER is used for ACIS and --pixadj=NONE or --pixadj=EXACT for HRC simulations.

The file may be viewed in ds9:

unix% ds9 tCanMaj.fits &

Figure 1: The events file created by MARX

[Thumbnail image: The events file created by MARX includes the ACIS readout streak.]

[Version: full-size]

[Print media version: The events file created by MARX includes the ACIS readout streak.]

Figure 1: The events file created by MARX

MARX includes a simulation of the ACIS readout streak in the events file.

[TIP]
Include ACIS Pileup Effects?

If the objective is to simulate an observation as opposed to generating a model PSF for an observation, then you may want to include pileup effects. First you need to be sure that your simulation parameters exactly match the observation (for example the spectrum flux was not artifically increased to increase sampling of the PSF). If these are not correct, the observed pileup cannot be matched to the simulation.

To include pileup, you must run the marxpileup tool:

unix% cp $MARX_ROOT/share/marx/pfiles/marxpileup.par .
unix% marxpileup MarxOutputDir=tCanMaj_marxsim.dir

When creating the events file, to include the pile up effects, use the --pileup flag to process the pileup simulation with marx2fits.

unix% marx2fits --pixadj=EDSER --pileup tCanMaj_marxsim.dir/pileup tCanMaj_pileup.fits

Note: With pileup enabled, it is unclear how to apply an ad-hoc normalization to get the simulation to match the data. One can try to use a region away from the core of the PSF; however, if the observed data does have any extended emission (dust halo, PWN, etc) then the result will be incorrect.

*In MARX 5.0.0, marxpileup introduces an offset on the source position, placing the source between 0.2-0.4 arcsec at 40 degrees from the readout streak, in the direction of the streak. This bug was fixed in MARX 5.1.0.

This pseudo-events file can now be used, for example, to create an image of the PSF.


Sub-pixel Analysis

We can now go ahead and demonstrate an example of using the PSF for imaging analysis.

Since the observed and simulated events files had the EDSER algorithm applied, sub-pixelated images may be created. We will create a 128x128 image, binning by a fifth of the native pixel size. We will use the source's position as the image center, first obtaining its corresponding sky position with dmcoords and then determine the image binning factor.

unix% punlearn dmcoords
unix% dmcoords acisf04469_repro_evt2.fits asol=pcadf04469_000N001_asol1.fits \
? opt=cel celfmt=deg ra=109.6769 dec=-24.954383

unix% pget dmcoords x y
4140.048508882587
4102.357195442711

# x-direction:
unix% echo "4140.0-12.8" | bc -l
4127.2
unix% echo "4140.0+12.8" | bc -l
4152.8

# y-direction:
unix% echo "4102.4-12.8" | bc -l
4089.6
unix% echo "4102.4+12.8" | bc -l
4115.2

This provides the necessary information to create the sub-pixel images of the observed source and its PSF, 4469_bin0.2.img and tCanMaj_bkgsub_bin0.2.img, respectively.

unix% dmcopy infile="acisf04469_repro_evt2.fits[bin x=4127.2:4152.8:0.2,y=4089.6:4115.2:0.2]" \
? outfile="4469_bin0.2.img" 

unix% dmcopy infile="tCanMaj.fits[bin x=4127.2:4152.8:0.2,y=4089.6:4115.2:0.2]" \
? outfile="tCanMaj_bkgsub_bin0.2.img" 

The image of the source can be restored with the arestore tool by deconvolving the PSF from the observed source using the Lucy-Richardson method, which does a good job with bright point sources.

unix% arestore infile=4469_bin0.2.img \
? psffile=tCanMaj_bkgsub_bin0.2.img \
? outfile=4469_bkgsub_bin0.2_deconv.img \
? numiter=256

Figure 2: Deconvolved Images of τ Canis Majoris

[Thumbnail image: The deconvolved image of the source by the MARX-generated PSF.]

[Version: full-size]

[Print media version: The deconvolved image of the source by the MARX-generated PSF.]

Figure 2: Deconvolved Images of τ Canis Majoris

τ Canis Major as observed by the ACIS detector (left), the MARX simulated PSF with the background subtracted spectrum (center), and the image deconvolved by the PSF of the source (right). The three images are on the same angular scale.

To the south of the source in the deconvolved image, there appears to be an unresolved blob, which in all likelihood is not physically real. There is a likely chance that this is an optical artifact from the HRMA that has been identified by the Chandra calibration team and is under investigation. This artifact only affects images on sub-arcsecond scales. The make_psf_asymmetry_region tool can be used to create a region file to overlay over an image of where the artifact will appear on the detector.

unix% make_psf_asymmetry_region infile=acisf04469_repro_evt2.fits \
? outfile=4469_asym.reg x=4140.0485 y=4102.3572

Figure 3: HRMA Artifact Observed in a Deconvolved Image

[Thumbnail image: HRMA artifact in the deconvolved τ Canis Major image.]

[Version: full-size]

[Print media version: HRMA artifact in the deconvolved τ Canis Major image.]

Figure 3: HRMA Artifact Observed in a Deconvolved Image

The region file produced by make_psf_asymmetry_region encompasses the extended feature seen in the deconvolved image.

The region overlaid on the deconvolved image corresponds to the expected location of the HRMA artifact, indicating that the feature is most likely not physically related to the source.


Encircled Counts Fraction

To get an idea of how well the PSF simulates the source, it may be worthwhile to examine the distribution of counts by looking at the encircled counts fraction (ECF) profile.

unix% ecf_calc infile=acisf04469_repro_evt2.fits outfile=tCanMaj.ecf \
? radius=8.2 xpos=4140.0485 ypos=4102.3572 binsize=0.2 mode=h

unix% ecf_calc tCanMaj.fits outfile=psf.ecf \
? radius=8.2 xpos=4140.0485 ypos=4102.3572 binsize=0.2 mode=h

The ecf_calc tool is applied on the observed and simulated events files to obtain the ECFs from a circular region, 8.2 sky pixels in radius, centered on the source position in sky coordinates. The profile is binned to match the earlier analysis of 0.2 pixels.

Figure 4: Encircled Counts Fraction Profiles

[Thumbnail image: ECFs of observation and PSF]

[Version: full-size]

[Print media version: ECFs of observation and PSF]

Figure 4: Encircled Counts Fraction Profiles

The encircled counts fraction profiles of τ Canis Major (white); the MARX-generated PSF (red); and the average of one-thousand separate realizations of the PSF (green) as a function of radius in pixels. The error bars are derived from 3σ uncertainties.

In this particular realization of the PSF, the MARX generated rays produce a slightly wider core than the observed source within a one pixel radius, but the observed core becomes slightly broader within a three ACIS pixel radius. The wings of the PSF is roughly equivalent to that of the observed.


Summary

This thread demonstrated how to generate ray tracing simulations with MARX and produce pseudo-events files from the results.


History

31 Jul 2014 original version
12 Feb 2016 updated to use MARX 5.2.0.
06 Apr 2016 revised AspectBlur values.
11 Apr 2016 Added warning about installing MARX on OSX with default C compiler.
12 Apr 2016 added information on using simulate_psf.
10 Jun 2016 Revised information about handling the SIM offset and added information about the MARX AspectBlur parameter.
20 Jun 2016 Fixed typos for suggested values in link to AspectBlur why page.
10 May 2017 Updated running marx2fits with --pileup.
17 Aug 2017 fixed typos.
27 Oct 2017 added admonition about simulate_psf input spectrum restrictions.
30 Mar 2018 Removed warning about clang compiler on OSX.
01 Aug 2019 added bash setup instructions.
26 Jan 2021 added conda setup instructions.
31 Jan 2022 Review for CIAO 4.14. Updated for Repro5 and CALDB 4.9.6.
21 Mar 2023 Replaced references to ${MARX_DIR} to ${MARX_ROOT}.
27 Apr 2023 revised AspectBlur description.
23 Oct 2024 revised and simplified Python/NumPy filtering syntax in input spectrum restrictions admonition.
09 Dec 2024 Added note about MARX installation; happens by default with ciao-install in CIAO 4.17 or can be installed with conda.