{'fraction_non_alphanumeric': 0.036408465173618246, 'fraction_numerical': 0.013149784261351962, 'mean_word_length': 4.529082096334444, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte-Carlo Simulator and Ancillary Response Generator of Suzaku XRT/XIS System for Spatially Extended Source Analysis

We have developed a framework for the Monte-Carlo simulation of the X-Ray Telescopes (XRT) and the X-ray Imaging Spectrometers (XIS) onboard Suzaku, mainly for the scientific analysis of spatially and spectroscopically complex celestial sources. A photon-by-photon instrumental simulator is built on the ANL platform, which has been successfully used in ASCA data analysis. The simulator has a modular structure, in which the XRT simulation is based on a ray-tracing library, while the XIS simulation utilizes a spectral"Redistribution Matrix File"(RMF), generated separately by other tools. Instrumental characteristics and calibration results, e.g., XRT geometry, reflectivity, mutual alignments, thermal shield transmission, build-up of the contamination on the XIS optical blocking filters (OBF), are incorporated as completely as possible. Most of this information is available in the form of the FITS (Flexible Image Transport System) files in the standard calibration database (CALDB). This simulator can also be utilized to generate an"Ancillary Response File"(ARF), which describes the XRT response and the amount of OBF contamination. The ARF is dependent on the spatial distribution of the celestial target and the photon accumulation region on the detector, as well as observing conditions such as the observation date and satellite attitude. We describe principles of the simulator and the ARF generator, and demonstrate their performance in comparison with in-flight data.


Introduction
A Monte-Carlo simulator is useful in characterizing a detector, and can relatively easily take into account many of the parameters which affect observations. Since the ultimate goal of X-ray data analysis is to estimate the true time, energy and position of the incoming X-ray photons, it is quite important to predict precisely how the photons interact with the telescope and detector. A good simulator is therefore strongly required not only for instrumental calibration and proposal planning, but also for scientific analysis. Chandra and XMM-Newton also have good simulators, named MARX 1 and SciSim 2 , respectively.
The X-ray observatory Suzaku (formerly known as Astro-E2) is the fifth Japanese X-ray astronomy satellite (Mitsuda et al. 2006). It has been developed under a Japan-US international collaboration, and was launched [Vol. , on 2005 July 10. Five X-Ray telescopes are present, sensitive to soft X-rays below ∼ 10 keV (XRTs; Serlemitsos et al. 2006). At the foci of four of the XRTs (XRT-I) are charge-coupled devices (CCD), known as the X-ray Imaging Spectrometers (XIS; Koyama et al. 2006); one (XRT-S) is combined with an X-ray calorimater known as the X-Ray Spectrometer (XRS; Kelley et al. 2006; XRS quit operation ∼ 1 month after the launch).
The Suzaku XRT is characterized by large collective area and relatively short focal lengths, compared with those of Chandra and XMM-Newton. In combination with these features, the low-earth orbit of Suzaku, where the particle background is low and stable, makes the non-Xray background (NXB) of Suzaku much lower compared with Chandra and XMM-Newton. In addition, the XIS achieves good spectral resolution, especially at the low energy range below ∼ 1 keV with the backside-illuminated (BI) CCD for XIS1. The front-illuminated (FI) CCDs for XIS0, XIS2, and XIS3 exhibit about half of the NXB rate than XIS1 (and less at energies > ∼ 8 keV), so they are complementary. Therefore, Suzaku has a unique advantage for spectroscopic observations of spatially extended sources (Mitsuda et al. 2006).
To achieve large collective area within the tight weight budget (1706 kg), the Suzaku XRT adopts the conical approximation of Wolter type I optics with 175 layers of the thin-foil-nested reflectors per quadrant (Serlemitsos et al. 2006). In return for the high throughput, it provides a moderate imaging capability of 2 ′ half power diameter with a complex point spread function (PSF), as well as the energy-dependent vignetting effects common to X-ray telescopes. In addition, there exists spatially-dependent contamination on the optical blocking filters (OBF) of the XIS (Koyama et al. 2006). These XRT and XIS characteristics often make extended source analysis complicated, so it is crucial to prepare a tool in order to precisely evaluate the effect of complex telescope and detector responses.
We developed a Monte-Carlo simulator of the Suzaku XRT/XIS system, which is incorporated into two practical tools, the XIS simulator "xissim" and the "Ancillary Response File" (ARF) generator "xissimarfgen". The simulator is constructed on the "ANL" platform ( § 2.1), which is used for almost all of the processing and analysis software of Suzaku.
While these tasks provide vast flexibility to the Suzaku XIS users, it is rather difficult to utilize them efficiently and appropriately. For example, there are more than 90 parameters for both xissim and xissimarfgen. There are several issues and limitations that one should be aware of in running these tasks. This paper is aimed to clarify these things by explaining principles of the software and by demonstrating performance with practical examples. We have also tried to separate the 'calibration issues', which can be changed (usually improved) by calibration updates, from those originated in the design of the software itself. The quality of the calibration is out of scope for this paper, although some aspects are discussed briefly in § 6. This paper is organized as follows. In § 2, we briefly show the strategy of the Suzaku software development, focusing on the ANL platform and simulators. In § 3 and § 4, we describe principles of xissim and xissimarfgen, respectively. In § 5, several notes on these tasks are described.
In § 6, we demonstrate these tasks with three distinct examples, the Crab Nebula, the North Ecliptic Pole (NEP) field, and Abell 1060. Finally, a summary is given in § 7. We also added three appendices which describe the coordinates definition, structures, parameters, and the output file formats, in detail.

The ANL Platform
When ASTRO-E software development started in 1995, the goal was a common software framework/platform which is used by realtime quick-look, data processing, and scientific analysis, both during pre-launch phase and after the launch. To that end, it was necessary to provide a common programming environment where instrument team members can easily develop, maintain and updates softwares that they need. This framework/platform also must allow end-users to share these softwares. The framework must be easy to learn for instrument team members, who, spending most of the time in calibrating the instruments, do not necessarily have extensive programming experience. Also, from the end users' point of view, it is desirable that those software tools developed based on this framework are maximally flexible and have an Ftoolslike simple interface which is familiar to most X-ray astronomers.
A software platform called "ASCA ANL", which had been developed for the ASCA satellite (Tanaka et al. 1994), fulfills these requirements. The ASCA ANL platform mandates modular design of the analysis software to be built upon it, and makes the software products easily configurable and reusable in components, so that software developers and end-users can share the same components for different purposes. This feature not only reduces code duplication, but also helps to quickly mature and refine the software.
Indeed, ASCA ANL fostered many practical tools including the instrument simulator SimASCA and the response generator SimARF. The advantages of the ASCA ANL platform are demonstrated by original scientific research which would have been difficult without SimASCA and SimARF; e.g., spatial-spectral analysis of clusters of galaxies (Ikebe et al. 1996;Honda et al. 1996), systematic analysis of large volumes of X-ray surveys and the cosmic X-ray background (Ueda et al. 1998;Ueda et al. 1999;Kushino et al. 2002). The SimASCA and SimARF were very helpful to realize a specialized analysis method in the analysis of spatially extended sources (which is not supported by standard analysis software), and to accurately compute complicated instrument responses.
On the other hand, however, ASCA ANL and other relevant software were based on the functions and libraries used for realtime quick-look software which had been de-veloped by the instrument teams, independently from the official ASCA analysis software (Ftools). This resulted in two independent streams to calculate basic physical values from the raw data, such as the pulse height corrected for the detector gain changes (known as "pulse invariant", or PI, corresponding to the detected photon energy), and the sky and detector coordinates of events, which caused confusion in the scientific analysis of the ASCA data.
Based on the ASCA experience, we adopted the "ANL", i.e. a generalized version of the ASCA ANL, as the software development platform for ASTRO-E and Suzaku. A brief history and concept of the ANL are described in Ozaki et al. (2006). At the same time, we developed a mechanism to convert ANL software directly to Ftools, to ensure that the ANL tools used for calibration by instrument teams are equivalent to Ftools used for pipeline processing and scientific data analysis. Common FITS-read and -write ANL modules and functions were also developed to handle photon event files and the calibration files in FITS format. 3 Now, almost all of the Ftools for Suzaku including xissim and xissimarfgen, released from the Guest Observer Facility at NASA/GSFC, are developed in the ANL framework.

History of Suzaku Simulators
The development of the Suzaku simulator had started before the failure launch of the ASTRO-E on 10 Feb 2000, especially for the bilinear 16 × 2 pixel XRS detector (Kelley et al. 1999). The detector size was comparable with the angular resolution of the XRT (Kunieda et al. 2001;Shibata et al. 2001), and so the XRS PSF was undersampled. The energy resolution of the XRS was also very dependent on the count rate of each calorimeter pixel. Therefore, the XRT/XRS system simulator, xrssim, was required to estimate the flux coming to each XRS pixel, which was critical for proposal planning. This was also the case for the Suzaku XRS. The xissim task subsequently developed by replacing the XRS component of the simulator with the XIS component.
The XRT ray-tracing part of the simulator has been significantly updated from the ASCA era. The code had been rewritten, by R. L. Fink (NASA/GSFC), from Fortran into C++, and the structure had been re-designed to utilize the mirror geometry and reflectivity files as separate calibration FITS files. This ray-tracing code is now supplied as the "xrrt" library by the XRT team. It has been utilized for the performance improvement of the XRT (Misaki et al. 2005) and the design of the pre-collimator to suppress the stray-light (Mori et al. 2005).
At present, xissim and xissimarfgen are publicly released in the Suzaku Ftools and all the source code is available including the ANL itself and the xrrt library. The latest version of the xissim package is 2006-08-26, which will be included in the next official release of the Suzaku Ftools for version 2.0 processing of Suzaku archival data scheduled in late 2006. All the calibration 3 The calibration FITS files are released to the public from the NASA/GSFC guest observer facility, as part of the official calibration database (CALDB; George et al. 1991).
information currently available is taken into account via the CALDB calibration database. The mkphlist, xissim, xissimarfgen, and xiscontamicalc tasks described in this paper are based on this version of the xissim package. The latest information on the xissim package is available at http://www-x.phys.metro-u.ac.jp/˜ishisaki/xissim/.
We also note that there is another Monte Carlo simulator for Suzaku, based on the Geant4 toolkit (Geant4 Collaboration et al. 2003) with ANL++ (Ozaki et al. 2006). This can simulate interactions of cosmic rays (both X/γ-rays and particles) with the satellite materials, such as satellite structures, shielding around detectors, and the detectors themselves. The main purpose of this simulator is to study response of the Hard X-ray Detector (HXD; Takahashi et al. 2006;Kokubun et al. 2006), and the NXB models for both HXD and XIS. See Terada et al. (2005) and Ozaki et al. (2006) for details.

Simulator: xissim
The xissim task simulates the interaction of the incident X-ray photons with the XRT/XIS system, using the XRT ray-tracing library and a spectral "Redistribution Matrix File" (RMF; see § 4.1) for the XIS, and generates a simulated event file. The format of the generated event file is a stripped-down version of that created by the pipeline processing of a real observation, so that users can analyze the simulated data in the same manner as the real data. To perform the simulation, users need to take three steps. First, the spatial distribution and the energy spectrum of the celestial source to be simulated must be specified. Second, a list of incident photons from the source needs to be prepared as FITS file(s). An auxiliary tool, mkphlist, may be used for this purpose. Third, the photon FITS file is passed to xissim, which then performs a photon-byphoton simulation, and creates a file of events detected by the XIS. In the following subsections, we describe how xissim performs the simulation.

Photon Generation
An auxiliary task mkphlist generates a list of faked photons of an X-ray source from the model spectral energy and spatial distribution, photon flux (in photons cm −2 s −1 ) in an arbitrary energy band, and the geometrical area of the XRT (in cm 2 ), provided by users. A model spectral distribution file (which is specified by the qdp spec file parameter) must be in units of photon flux (photons cm −2 s −1 keV −1 ) that can be easily produced with standard software packages such as Xspec (Arnaud 1996). mkphlist requires celestial coordinates of the point source or a surface brightness map (FITS image) on the sky for the spatial distribution of the source. Either the number of photons or exposure time is needed to determine how many photons are to be generated. Users can also specify equal or random interval steps for the photon arrival time. The structure of mkphlist is explained in Appendix 2.1, and a list of major parameters and the format of the photon file are summarized in table 5 and 8, respectively. Note that, by preparing an appropriate pho-ton file, users can in principle simulate any source with any energy spectrum and/or any spatial distribution.

Photon-by-Photon Simulation
By taking into account the XRT and XIS response, xissim performs photon-by-photon simulation for given input photon file(s). It has the capability to read up to eight photon files simultaneously. Figure 1 shows the schematic structure of the simulation implemented in xissim. Since understanding the coordinate systems is essential, we include the definitions in Appendix 1.
First of all, the ra and dec values in the photon file need to be converted to (θ, φ), i.e. offset angle ( ′ ) and azimuth angle ( • ), with respect to the XRT optical axis. This requires the satellite Euler angles (ea 1 , ea 2 , ea 3 ) ( • ), the observation date for an aberration correction (or parallax correction, see Appendix 1), and the alignment parameters in the telescope definition (teldef ) file. Users can supply an attitude file (set of Euler angles as a function of time) and a good time interval (GTI) file to take into account the wobbling of the spacecraft (See also § 5.2 for the attitude wobbling). The photon time column in the photon file usually starts from 0.0 s unless otherwise specified, and it is treated as the time offset relative to the GTI. Alternatively one may specify a fixed set of Euler angles and/or a fixed date. The aberration correction can be disabled by setting the hidden parameter aberration=no (hidden parameters are not required when invoking an Ftools task).
In the second stage, the geometrical area for a given photon is reduced by a factor of cos θ due to the slanted incidence to the XRT. This factor is usually very close to unity, and had been neglected in the older version of xissim. This behavior can be controlled with the parameter aperture cosine, and is set to 'yes' by default in the present version. The photon flux is further reduced due to transmission through the thermal shield on the top of the XRT. Xissim then assigns a random location for each photon at the top surface of the XRT, where the precollimator is placed. The task traces the path of each photon inside the XRT (pre-collimator, primary and secondary mirrors), using the XRT ray-tracing library, xrrt (Misaki et al. 2005;Mori et al. 2005), using the XRT geometry and reflectivity as described in the ray-tracing code and the calibration files. After the ray-tracing, some photons may be absorbed and disappear, while others reach the focal plane.
A fraction of the photons that have reached the focal plane are absorbed by the contamination on the OBF. The thickness of the contamination is time-and detectorposition-dependent (Koyama et al. 2006), and their dependence is given by a calibration file supplied by the XIS team. Xissim computes the transmissivity at a given time and position using this calibration file. The position of the photon on the detector is again calculated by the alignment parameters in the teldef file.
Finally, the simulated photons reach the detector (including both OBF and CCD), where the detection probability is determined using the RMF of the XIS. The XIS RMF contains the transmission of the OBF and the quantum efficiency of the CCD, as well as the spectral redistribution matrix from energy to PI. The line response function of the XIS CCDs is primarily a Gaussian distribution but it also includes other features such as escape ratios and tails that deviate from a Gaussian. Photons that have passed the test for detection are recorded as X-ray events, 4 and their PI values are determined from the incident photon energy by random choices according to the energy redistribution probability in the RMF. The Suzaku XIS detectors do not exhibit significant positional dependence in the energy resolution after the CTI (charge transfer inefficiency) correction while the energy resolution is known to degrade with time. Users should supply an appropriate RMF corresponding to the observation date, which can be generated by a separate task, xisrmfgen.
Note that the current version of xissim does not consider the NXB, bad CCD columns, event pile-up, event grade, nor CCD exposure frames. Although the output event files contains the same major columns as the event files of the real data, the status and grade columns are filled with 0, while the pha column has the same value as the pi column.
The ANL module structure and parameters of xissim is explained in Appendix 2.2 in detail. If one has the ANL programing environment available, he/she may add his/her own modules to the simulator. It is also easy to replace a module, e.g. if a module is available that more precisely simulates the CCD detection process, then this module can be substituted for the SimASTE XISRMFsim module which utilizes a ready-made RMF. This is one of the great benefits of the ANL. Table 1 summarizes the list of calibration files used by xissim. The file specified by the leapfile parameter is the leap second file, and is required to compute the mission time (or Suzaku time), defined as accumulative seconds since 2000 January 1, 00:00:00 (UTC). In fact, the default value of the leapfile parameter is set to a special keyword of "caldb" and is a hidden parameter. By installing CALDB and properly setting the environmental variables, the xissim task automatically searches the most recent calibration file for this category, i.e. 'Content Name' = leapsecs. This is also applicable to other parameters in table 1.

Calibration Files
The shieldfile, mirrorfile, reflectfile, and backproffile parameters are used for the XRT simulation. There are four FITS extensions in the mirrorfile to describe the geometry of each XRT, and three extensions are present in the reflectfile corresponding to materials of the reflection surface. As described in Appendix 1, teldef is used to describe the mutual alignments between XRT and XIS, as well as among the XIS sensors and the spacecraft. The contamifile describes the energy, time, and position dependence of the contamination on the XIS OBF.

Ancillary Response Generator: xissimarfgen
The xissimarfgen task generates a Suzaku XIS ARF based on user-defined conditions, such as an arbitrary shape of the X-ray emitting region and event extraction regions. Xissimarfgen does so by simulating photon detections at each energy. It then calculates the detection efficiency in a user-defined event accumulation region. 5 Since it utilizes a Monte-Carlo simulation, users need to simulate enough photons to avoid counting statistics errors. It can refer to the attitude file to reflect the change of effective area due to the attitude wobbling. The final ARF is in the standard FITS format, so that users can use Xspec or other standard fitting packages for spectral analysis.

Principle of ARF Calculation and Limitations
The ARF is utilized for spectral fitting combined with an RMF. See George et al. (1992) for detailed format 5 We shall use 'accumulate' for the simulated events, and 'extract' for the real observed events. of these files. The RMF is represented by an (m × n) matrix R(E i , PI j ), where E (keV) denotes the energy and PI (channel; hereafter chan) denotes the pulse invariant, with 1 ≤ i ≤ m and 1 ≤ j ≤ n. Regarding the XIS, m = 7900, n = 4096, E 1 = 0.201 keV, E m = 15.999 keV, PI 1 = 0 chan, and PI n = 4095 chan for the nominal RMF. The ARF is represented by an m-dimensional vector which we denote as S A(E i ) (cm 2 ), where S = 1152.41 cm 2 represents the geometrical area of the XRT. The goal of the spectral fitting is to find a model spectrum, M (E i ) (photons cm −2 s −1 keV −1 ), which fits the observed spectrum, D(PI j ) (count cm −2 s −1 ). The response and model spectrum are convolved, i.e., where ∆E (keV) is the energy bin width, and M(PI j ) and D(PI j ) are compared. As one can see easily from represents an expected spectrum for the monochromatic X-ray of E = E i keV, and Thus, calculating the ARF is reduced to the computation of the detection efficiency at each energy step, E i , of the RMF, a job well-suited for a Monte-Carlo simulation. For a given input N in counts of monochromatic X-ray photons at E = E i keV, the simulator predicts N det detected events and then the detection efficiency is simply A(E i ) = N det /N in . However, one must be very careful because the detection efficiency, namely N det , is influenced by many factors: first of all, the accumulation region of the event on the detector, and the spatial distribution of the celestial sources assumed on the sky. It is also affected by the observational conditions, such as the satellite Euler angles, the date of the observation due to the thickness of the XIS contamination and the parallax correction, etc. The quality of the calibration and/or the Poisson statistics can also impact N det . It is therefore important that one must reproduce the user-selection and the observational conditions of the real data as much as possible in the simulation. One must also take care to perform a simulation such that the photon statistics are sufficiently better than the statistics of the real observation.
In fact, the spatial distribution on the sky is sometimes complex and/or extended on a scale larger than the telescope FOV. Thus the accuracy of the spatial model can become a major cause of systematic error in the estimation of the detection efficiency, which leads to uncertainty in the source flux. For example, if one assumes a more core-concentrated image than in reality, more photons will be simulated to arrive at the detector, which will overestimate the detection efficiency. One can test the assumed spatial distribution on the sky by comparing the real observation image and the simulated one.
There is also another limitation due to the spectral fitting procedure itself. In the conventional spectral fitting package (e.g., Xspec v11 or before), one can choose only a single response matrix (ARF + RMF) for an observed spectrum in the spectral fitting. 6 For example, an observed spectrum may contain thermal emission which obeys an oval surface brightness profile, as well as the cosmic X-ray background (CXB) spectrum of a Γ ∼ 1.4 power-law which extends nearly uniformly on the sky. The ARF response for the oval surface brightness is different from that for the uniform-sky emission, hence one cannot fit the observed spectrum with the thermal model + power-law model in a usual way. Strictly speaking, the energy spectrum should be the same at every point in the assumed spatial distribution on the sky in order to conduct spectral fitting with a single ARF + RMF response.

Implementation of ARF Calculation
As described in § 3.2 and shown in figure 1, the XIS RMF takes care of the OBF transmission and the quantum efficiency of the CCD, hence the XIS ARF should consider other factors for the detection efficiency, namely, the thermal shield transmission, XRT effective area, transmission of the OBF contaminant, etc. Detailed explanation of structure, parameters, and the output ARF format are given in Appendix 2 and 3 It reads a number of parameters which specify the simulation conditions (table 7), and (1) determines energy steps to calculate detection efficiency; (2) generates monochromatic photons (or quasi-monochromatic within the narrow energy range) until the user-specified condition on the photon statistics is fulfilled at each energy step; (3) conducts the ray-tracing simulation for each photon; (4) counts up the number of detected events at each energy; (5) records the detection efficiency at each RMF energy bin to the output ARF(s) by interpolating the simulation result; (6) continues to the next energy step and loops to step (2).
Note that the energy step determined in step (1) is usually not same as the RMF energy bin, because the com-6 This restriction no longer holds in the latest release of Xspec v12, which allows different model components to have their own response. putation time would be very long to conduct photon-byphoton simulations in standard XIS RMF 2 eV steps up to 16 keV. Interpolation is therefore required in step (5). In addition, the XRT effective area usually changes only gradually with energy except for several characteristic energies such as the Au-M, Au-L, and Al-K edges, 7 so that we may often choose sparse energy steps. This feature can save the computation time effectively.
In the calculation of the detection efficiency, the three major factors of (i) transmission of the XRT thermal shield, (ii) effective area (cm 2 ) of the XRT, and (iii) transmission of the XIS OBF contaminant are treated separately. They are also written in separate columns in the resultant ARF as shield transmis, xrt effarea, and contami transmis (table 10). The resultant detection efficiency times the geometrical area, S A(E i ) (cm 2 ), is written in the specresp column, i.e., specresp = shield transmis × xrt effarea × contami transmis. Note that (i) and (iii) are supplied in the calibration files (table 1) in fine energy steps of ∼ eV, whereas (ii) is usually calculated in more sparse energy step. By separating these factors, one can obtain a good quality ARF even in a sparse energy step for the simulations, and moreover, one may remove, scale, or multiply the contami transmis factor afterwards. The xiscontamicalc task is provided to do this kind of the ARF manipulation.
Note that the thickness of the OBF contaminant is positionally dependent. It is therefore required to know the spatial distribution of photons 8 falling on the OBF at each RMF energy bin in order to evaluate the contami transmis factor. This energy dependence of the photon distribution is also determined by interpolation, which incurs additional calculation time when the simula-7 The front surface of the XRT reflector is coated with gold and its substrate is made of aluminum. The pre-collimator is made of aluminum, too, so that the Al-K edge appears in the largeoffset-angle response of the XRT.

8
This distribution is approximated in xissimarfgen by a DET coordinate image binned prior to applying absorption due the XIS contamination. tion energy step is much wider than the RMF energy bin. It is not easy to estimate the true photon distribution from the real observation data, because the observed image is affected by the XRT vignetting and the OBF contaminant, both of which are energy dependent. Vignetting is a more severe effect in the higher energy band, and the OBF contamination is severe in the lower energy band. In addition one must subtract background to utilize the observed image. The combined energy and spatial dependence of the XIS contamination is considered in the ARF generator rather than the RMF generator, for this reason.
At each simulation energy, in fact, the A(E) value is calculated using the weighted sum of events, N w , instead of N det , as, (2) in which w k (E) denotes the weight value (see Appendix 2.2) of each simulated photon at the energy of E keV. As mentioned above, the resultant A(E i ) values at the RMF energy bin are calculated by interpolation, complicated somewhat by when the transmission of the OBF contaminant is considered. Here, we de- The transmission of the OBF contaminant is calculated for each of the simulated photons, as τ k (E i , photon time k , detx k , dety k ). Note that the energy of each simulated photon, E ′ l = photon energy k , has been replaced by the energy of the RMF bin, E i . Thus xissimarfgen computes the final detection efficiency, A(E i ), with contamination, as by an interpolation. The definitions of s i and t i are given in eqs. (A2) and (A3). It also calculates the relative error of A(E) at each simulation energy, and the interpolated values are stored in the relerr column of the output ARF (table 10). This column is useful to judge the photon statistics is sufficient for the ARF calculation. The relative error is calculated as, if N in N det (N in − N det ) = 0, otherwise relerr = 1.0. The derivation of this formula is a little tricky, because we know the detected count N det and the undetected count N in − N det in the simulation, and both are considered to follow the Poisson statistics. Since A(E) is expressed as We therefore defines the relative error as

Notes
In this section, we describe several notes on xissim and xissimarfgen. § 5.1 applies to both tasks, and others apply mainly to xissimarfgen.

Notes on Random Numbers
The quality of the random number generator to be used can affect the quality of the Monte-Carlo simulation results. A good random number generator should include a very long cycle, fast computation, and wide significant bits. xissim and xissimarfgen use an internal random number generator in the astetool library, 9 utilized by all modules. This generates double precision floating point values in the range of 0 ≤ r < 1.0 based on the Tausworthe method (Tausworthe 1965). The generated random number has 62 significant bits (≃ 4.6 × 10 18 ) and its cycle is estimated to be about 2 250 ≃ 10 75 . These parameters are significantly wider and longer than the usual random number function, int rand(void), implemented in the standard C library.
Its code is machine independent, and it reproduces exactly the same series of random numbers as long as the rand seed and rand skip parameters are the same. It is recommended to set a prime number (except 2) to rand seed for good randomization. The default value of rand seed for the simulation tasks is 7. They record the number of random numbers generated in the simulation to the output event file, as the randngen keyword in the FITS header. One may re-continue the simulation with the same series of random numbers by setting the rand skip parameter to its value. However this code is not multi-thread compliant, which may need to be upgraded in the future for faster (i.e., distributed) simulations.

Notes on Accumulation Region
There is a difference between specifying the accumulation region in SKY coordinates versus DET coordinates. This may be ignored only when the attitude wobbling and the parallax correction are negligible. The accumulation region is fixed on the CCD when it is specified in DET coordinates. On the other hand, it moves around the CCD when specified in SKY coordinates, according to the attitude wobbling. In both cases, the celestial target moves around the CCD and is affected by the vignetting effect of the XRT, also due to the attitude wobbling. Xissimarfgen can treat both situations correctly, as far as the supplied attitude file is reliable, so that one should select the region mode parameter to match the extraction method of the real observation spectrum.
It is known that there is an unexpected attitude wobbling of ∼ 0.5 ′ due to thermal distortion of supporting structure (Serlemitsos et al. 2006), however this effect is not included in the present attitude file. This situation will be improved in near future by a dedicated Ftool, the aeattcor task. Until then, it is recommended to avoid 9 See also Appendix 1 for the astetool library. [Vol. , using too small of an accumulation radius (r < ∼ 3 ′ ). One may check this effect by changing the accumulation radius and test whether the fit results are affected significantly. Alternatively, one may track the position of the PSF core on the CCD for bright point-like source targets.
It is also notable that the background files for the XIS currently released, which are a collection of events when the XRT was pointed to the night (non-sunlit) Earth, do not support event extraction in SKY coordinates. This situation will be improved in the future. One may extract the background from the outer ring of the target, however, this region also contains the outskirt of the PSF of the main target, CXB, and the instrumental background, which have a small dependence on the detector position. The former two effects can be evaluated by xissimarfgen, and the latter can be tested with the released background file.

Notes on Flux Normalization
Because the detection efficiency defined in eq. (2) is considered for all the input photons coming from everywhere in the supplied source image, the normalization of the flux in spectral fitting gives the value integrated over the whole region of the source image. Therefore, if one generates a uniform-sky ARF with source rmax = 20 ′ to fit the CXB spectrum, the fit gives the flux from the π · source rmax 2 = 1257 arcmin 2 sky area, then the user needs to divide the flux by this area to convert it to a surface brightness.
Other cases can similarly be complex, e.g., an analysis of a cluster of galaxies. Extracting spectra from annular rings centered on the cluster core is frequently performed in the cluster analysis. Here, we assume that the cluster emission spectrum is identical everywhere on the sky, and only the normalization of the flux decreases as the distance from the cluster core increases. We also assume that the spatial distribution of the cluster on the sky can be perfectly predicted, which has been supplied to xissimarfgen as the source image. Then the fit results for each ring should give the same flux, while the observed count per unit area decreases as the ring radius increases, since the flux for the whole cluster is calculated for each fit.
If one gets different fluxes for each fit, then this is the 1st order approximation of the correction factor to the assumed source image at each ring. It is often desired to derive the flux only coming from each ring. To help with this kind of task, there is a keyword, source ratio reg, written in the output ARF (table 11). This keyword holds the ratio of the source image inside the specified accumulation region for the ARF, which has been calculated during the simulation. By multiplying this factor by the obtained flux, the user can calculate the flux in that ring.

Notes on Computation Time and Memory
The code of xissimarfgen is designed to conduct the computation of an ARF as efficiently as possible in both time and memory, although it still requires a significant amount of both. The simulation code has been tuned for speed; it reads all the required information including the attitude into memory before the simulation. Searches of tables such as reflectivity, transmission, spatial and spectral distributions are accelerated by adding an index. Several functions cache previous values to skip redundant calculations especially when the photon energy and/or time is similar to the previous ones. In addition, the binary distribution of the xissim/xissimarfgen package is compiled with fast C compilers using the highest optimization option. The required memory is usually around 130 MB, hence recent machines can easily run xissimarfgen task in memory. The actual calculation time is very dependent on the simulation energy step, the photon statistics, as well as the computer platforms. Table 2 shows examples of computation time on the AMD Athlon TM 64 2.4 GHz CPU with a 64-bit Linux OS. The example (A) is the ARF shown in § 6.1, in which full observational features are taken into account, and a Chandra image (1800 × 1800 pixels) was supplied for the source image with source mode=skyfits. Parameters of num photon=100000 and estepfile=sparse (55 energy steps) were chosen, so that 5.5 × 10 6 photons were simulated for the ARF calculation.
The time difference between (A) and (B) indicates that significant fraction of time (∼ 70 s) was consumed in the calculation of the XIS contamination. However, this time is only proportional to m × N det and does not depend on the simulation energy step, hence it should be acceptable. The parallax (aberration) correction also needs the nonnegligible cost of ∼ 9 s, which is proportional to N in = num photon. Similar time is needed for the randomization in the spatial distribution of the Crab nebula, as seen in (C) → (D). The consumed time in the XRTsim module is also displayed by the ANL, and it was 75.9 s for (D). This indicates that the ray-tracing code can perform the simulation of a single X-ray photon in less than 15 µs on this machine.

Demonstration
In this section, we demonstrate how xissim and xissimarfgen work using three distinct examples: the Crab nebula as a calibration source and a quasi-point-like source in § 6.1, the NEP field as a "blank sky" in § 6.2, and Abell 1060 as an spatially extended source in § 6.3.

Crab Nebula
First, we present the case of simulating the Crab nebula which is the main X-ray calibration source for effective area calibration (Toor & Seward 1974;Seward 1992;Kirsch et al. 2005). Small in angular scale, it has a complex spatial structure as seen in figure 2 (a) of the Chandra image (Weisskopf et al. 2000). With respect to the surface brightness map, we adopted this image, because Chandra's X-ray telescope, HRMA, has much superior angular resolution of ∼ 0.5 ′′ . We further compensated it manually for a point-like emission from the neutron star (K. Mori priv. comm.).
We made a photon list by supplying the image to mkphlist, and ran xissim with it. The simulated Suzaku image of the Crab nebula is shown in figure 2 (c). For comparison, we also present the simulated image for a pointlike source in figure 2 (b). The simulated Crab image appears as a smoothed PSF with the extent of the complex surface brightness profile of the Crab nebula. Figure 2 (d) shows the real observation image taken with the Suzaku XIS0 detector. The global extent of the Crab image is consistent with that of the simulated image. The anisotropy in the azimuth direction in the real image is due mainly to the complex PSF shape of the actual XRT, which will be more accurately reproduced by future improvement of the mirror geometry file (mirrorfile). Once the calibration file is updated, xissim can reflect it automatically via the CALDB. A narrow groove crossing the central area from east to west is due to a bad CCD column. The out-oftime events, which broadly spread on both the east and west sides, are also seen along the direction of the signal transfer from imaging area to frame-store region. These features are not implemented in the current version of xissim.
Using the Chandra image, we also generated an ARF for XIS0, and it is plotted in figures 3 (a) and (b). The specresp (black) and xrt effarea (green) columns are plotted in figure 3 (a), and the contami transmis (black) and shield transmis (green) columns are plotted in figure 3 (b). The 90% confidence range of the specresp is also drawn by cyan lines in figure 3 (a). Full observational conditions, namely attitude, gtifile, aberration, and contamifile, are considered in the ARF generation with num photon=100000 and estepfile=sparse (55 energy steps). The accumulation radius is 6 mm = 250 pixel ≃ 4.34 ′ in the SKY coordinate.
For comparison, we plot the nominal ARF in CALDB without contamination, ae xi0 xisnom6 20060615.arf, in red line. In fact, the nominal ARF was also generated by an older version of xissimarfgen, and the calibration files were not changed between the two versions. However, the nominal ARF is calculated with much denser energy step (2 eV steps below 4 keV, and at most 10 eV steps above 4 keV, with 3450 energy steps), and 4 times higher photon statistics (num photon=400000). Although slight jerks are seen in black and green lines in figure 3 (a), these two ARFs are quite consistent. Discrepancy in the lower energy range is due to the XIS contamination, which is plotted by a black line in figure 3 (b). Therefore, in the spectral fitting, the Crab nebula can be treated as a point-like source with Suzaku XIS, if the extraction radius is large enough (r ∼ 6 mm) and the spacecraft attitude is stable. We also note that the nominal ARFs give flux consistent with that obtained by Toor & Seward (1974) within ∼ 2% for all the XIS sensors (Serlemitsos et al. 2006).

NEP Field
The NEP field is an archetypal "blank field", where no X-ray bright objects exist. In such a region, the X-ray background, including both the extra-galactic (Brandt & Hasinger 2005) and the Galactic components (Snowden et al. 1995), is the dominant X-ray source. The X-ray background can be treated as almost uniform distribution, hence we tested the uniform-sky ARF generated by xissimarfgen in this field.
We created an ARF assuming a uniform distribution for the source from a circular region with a radius of 20 ′ (see caption of figure 4 for details of the parameters), and fitted the observed spectrum with it. In figure 4 (a), the effective area, specresp, of the obtained ARF is displayed in comparison with that for a point source. One can see that the effective area is relatively smaller in the higher energy band ( > ∼ 7 keV), due to the vignetting effect of the XRT. After subtracting the NXB contribution estimated using the night Earth database ( § 5.2), the spectrum can be well fitted with a power-law model representing the CXB and one or two thin-thermal plasma models representing local Galactic thermal components, as shown in figure 4 (b). The photon index and the surface brightness of the power-law component are consistent with the parameters reported so far (Gendreau et al. 1995;Kushino et al. 2002). See Fujimoto et al. (2006) for the details of the analysis. This result demonstrates that xissimarfgen properly generates the ARF for the uniform-sky emission.

Abell 1060 Cluster of Galaxies
Finally, we present an example of the Abell 1060 cluster of galaxies observed with Suzaku. Scientific results will be published by K. Sato et al. in preparation. Abell 1060 is a circular and nearly isothermal (∼ 3 keV) cluster of galaxies (Tamura et al. 2000;Furusho et al. 2001;Hayakawa et al. 2004;Hayakawa et al. 2006) and is suitable for testing the ARF for extended sources. There were two pointings performed with Suzaku at the central region and the ∼ 20 ′ east offset region, as shown in figure 5 (a). These observations were conducted at the end of November 2005, when the XIS contamination was already significant and was starting to saturate.
The observed spectrum is assumed to contain (a) thin thermal plasma emission from the intra cluster medium (ICM), (b) local Galactic emission, (c) CXB, and (d) NXB. We can estimate (d) using the night Earth database mentioned in § 5.2, and can subtract it from the observed spectrum. As demonstrated in § 6.2, the spectrum of (b) can be represented by the (apec + apec) model with 1 solar abundance, and that of (c) has a shape of absorbed power-law with Γ ≃ 1.4. However, we cannot  fit the observed spectrum directly by a sum of (a) + (b) + (c), because the spatial distribution of these three are different on the sky, as described in § 4.1. We therefore adopted the following method for the spectral analysis of Abell 1060. We extracted several spectra from annular regions centered on the cluster core, and here we show two samples of the innermost region at the projected radius of 0-2 ′ from the central observation, and the outermost region of 17-27 ′ from the offset observation, as representatives. We generated two different ARFs for each spectrum, A u (E i ) and A b (E i ), which respectively assume the uniform-sky emission and the ICM surface brightness profile obeying an analytical model obtained with the XMM-Newton data.
As described in § 4.1, it is important that the assumed spatial distribution on the sky well agrees with the actual data in the calculation of the ARF response. We therefore compared the observed images with the simulated ones in figure 5 (b) and (c). The 1-4 keV energy range was chosen so that the distortion of the image due to the XIS contamination and the XRT vignetting was not severe. In this energy range, the Galactic component (b) is almost negligible, whereas the CXB and NXB components cannot be neglected especially in the offset observation. The NXB component (red line) is estimated from the night Earth database. A small ACTY (= DETY −1 for XIS0, see figure 8) dependence of the NXB intensity is seen, reflecting the dwell time at the frame-store region of the CCD. The CXB component (blue line) is estimated by the xissim simulation, assuming the uniform sky and the previous ASCA results of the CXB intensity (Kushino et al. 2002). The vignetting effect is seen in the CXB counts, hence the count rate slightly drops at the CCD rim. After the subtraction of the estimated CXB and NXB components, the observed distribution of the cluster (black crosses) is fairly well reproduced by the xissim simulation of the cluster emission (green line), although a small asymmetry is observed for the real cluster in the central observation. Although the accumulation area is smaller for (a) than (b), the calculated effective area is much larger for (a) than (b) as seen in figure 6 (c), plotted in black and red lines, respectively, due to the assumed surface brightness profile. One can see the position dependence of the XIS contamination (thinner towards the CCD edge) is treated appropriately as seen in figure 6 (d).
Denoting the spectra of (a), (b), (c), and (d) as M icm (E i ), M gal (E i ), M cxb (E i ), and M nxb (PI j ), the observed spectrum can be expressed by a sum of, where the operator ⊗ denotes the transformation defined by eq. (1). It is known that the CXB spectrum, M cxb , is fairly constant over the sky except for difference in the neutral hydrogen column density, N H , for absorption, whereas the local Galactic emission, M gal , may vary from field to field by more than an order of magnitude (Kushino et al. 2002).
Considering this situation, we assumed a power-law spectrum for the CXB with the values by Kushino et al. (2002), Γ = 1.4 and S X = 5.97 × 10 −8 erg cm −2 s −1 sr −1 (2-10 keV), 10 measured with the ASCA GIS (Ohashi et al. 1996;Makishima et al. 1996). The neutral hydrogen column density was fixed to N H = 4.9 × 10 20 cm −2 (Dickey & Lockman 1990). We calculated the estimated contribution of the CXB, A u ⊗ M cxb , using the fake command in Xspec. This contribution for each region is indicated by blue crosses in figures 7 (a) and (b) for XIS1. We subtracted the CXB contribution from the observed spectrum as well as the estimated NXB spectrum. The XIS1 (red) and FI (XIS0+XIS2+XIS3; black) spectra in figures 7 (a) and (b) denote those after the CXB and NXB subtraction. We then fitted the XIS1 and FI spectra simultaneously for the offset observation, where the Galactic component (b) is prominent, with the [ apec (cyan) + apec (orange) + phabs × vapec (magenta) ] model, using the ARF response, A b (E i ).
In a strict sense, using A b for the component (b) is not correct, which should be A u instead. We made this choice due to the limitation of the Xspec v11 ( § 4.1), however, it does not matter practically if we only notice the shape of the spectrum. The absolute surface brightness of the Galactic component was evaluated separately using the Xspec fake command and the A u response. We then fitted the central region, fixing the shape of the Galactic component, but with its normalization scaled so that the surface brightness is preserved between the two different sky regions. Xspec v12 can handle this situation more straightforwardly.
The released RMF, ae xi[0-3] 20060213.rmf, was used for the spectral analysis. The ARFs were generated by xissimarfgen, and were convolved with the RMFs and added for three FI sensors, using the marfrmf and addrmf tasks in Ftools. As for the photon statistics of the simulation, limit mode=mixed was chosen with num photon=100000 and accuracy=0.005. As seen in figures 7 (a) and (b), both the observed spectra can be well fitted by one temperature plasma emission model for ICM, and (apec + apec) model with 1 solar abundance for the local Galactic emission. The surface brightness and the spectral shape of the Galactic emission is kept constant between both regions. This is confirmed by the fact the ratio of the Galactic components to the CXB is almost equal between figures 7 (a) and (b). Note that the fit appears equally good to the BI (XIS1) and FI sensors, which are different in the thickness of the contamination. So far, it has been confirmed that the temperature for each ring derived from the spectral fitting with this method is quite consistent with the previous results with XMM-Newton and Chandra (Hayakawa et al. 2004;Hayakawa [Vol. , et al. 2006). See K. Sato et al. in preparation for details of the results.

Summary
• We have developed a Monte-Carlo simulator of the Suzaku XRT/XIS system taking into account full calibration results.
• We adopted the ANL platform that provides us a flexible and comprehensive environment for the Suzaku software production.
• There is a dedicated task named mkphlist which generates a photon file to feed to the simulator.
• The task xissim reads the photon file, and conducts the instrumental simulation using the XRT ray-tracing library and the RMF of the XIS, and generates an event file, which is consistent with that for real observation, so that users can analyze the simulated data in the same manner as real data.
• The simulator-based ARF generator is named xissimarfgen, which can compute up to 200 ARFs corresponding to different accumulation regions by a single batch of simulations.
• The combination of xissim and xissimarfgen enables users to analyze spatially extended and spectroscopically complex celestial sources.
• Since one of the Suzaku's unique features is the low and stable particle background, these simulators are crucial for producing scientific results with low signal-to-noise data from extended sources.
• The latest public version is 2006-08-26, which will be included in the next official release of the Suzaku Ftools scheduled in late 2006.
Thanks are given to the referee, Dr. K. Arnaud, for useful comments which improved the original manuscript. We express sincerely thanks to H. Honda for early stage work on the ASCA ANL, and SimASCA. We also show R. L. Fink our appreciation for early stage work on the xrrt ray-tracing library. We acknowledge to L. Angelini and I. Harrus for useful discussion on the CALDB file format. We thank K. Mori to provide us a Chandra image of the Crab nebula corrected for pile-up photons. We also acknowledge D. McCammon, R. Smith, Y. Takei, C. Matsumoto, and N. Ota for testing and giving valuable comments on the xissim and xissimarfgen tasks. Part of this work was financially supported by the Ministry of Education, Culture, Sports, Science and Technology of Japan, Grant-in-Aid for Scientific Research No. 14079103, 15001002. Default image region. X/Y values can be outside of the region. ♯ Angular scale at the center. Outer pixels are slightly different due to the tangential projection.

Item
Ideal Value Focal length 4750 mm Optical axis location in DET (512.5, 512.5) Size of the DET pixel 0.024 mm/pixel Offsets between DET and FOC (0.0, 0.0) Roll angle between DET and FOC 0.0 deg Alignment matrix for FOC → SKY * 3 × 3 identity matrix * Alignment matrix is common to all sensors.

Appendix 1. Definition of the Coordinates
The following coordinates are defined to describe event locations in the telemetry, on the detector, or on the sky. RAW coordinates: Original digitized values in the telemetry to identify the pixels of the events. This may not reflect physical locations of the pixels on the sensor. For example, XIS RAWX (or RAWY) coordinate will have values from 0 to 255 (or 1023) on each CCD segment. Each of the four XIS sensors has a single CCD chip, and a single chip is divided into four segments. ACT coordinates: The ACTX/Y values are defined to represent actual pixel locations in the CCD chips. ACTX/Y will take 0 to 1023 to denote the 1024 × 1024 pixels in the chip. The XIS RAW to ACT conversion depends on the observation modes (such as Window Options) and will require housekeeping information. The XIS ACT coordinate is defined by looking down on the sensors, hence the ACTX/Y to DETX/Y conversion needs a flip in the Y-direction. DET coordinates: Physical positions of the pixels within each sensor, XIS0-3. Misalignments between the sensors are not taken into account. The DETX/Y coordinate are defined by looking up the sensor, such that the spacecraft (S/C) +Y direction becomes the −DETY direction (the same convention as with ASCA ). The S/C Z-axis points in the telescope direction, and +Y direction is toward the solar paddle. For XIS, the DETX and DETY values take 1 to 1024. FOC coordinates: Focal plane coordinate common to  (a) (b) Fig. 7. Example spectra of the Abell 1060 cluster of galaxies, (a) for the central observation, (b) for the offset observation. In both figures, red or black crosses represent the observed spectrum with the XIS1 (BI) or XIS0+XIS2+XIS3 (FI) sensor(s), respectively, for the upper panels, and the fit residuals for the lower panels. The estimated CXB + NXB spectrum has been subtracted from each observed spectrum, and the estimated CXB spectra for XIS1 are indicated by blue crosses. The spectra are fitted by the [ apec (cyan) + apec (orange) + phabs × vapec (magenta) ] model in Xspec 11.3.2t indicated by green line for XIS1 and yellow line for the FI sensors. The model components are only plotted for the XIS1 spectrum. all the sensors. Misalignments between the sensors are taken into account so that the FOC images of different sensors can be superposed. The origin of the FOC coordinate corresponds to the XIS nominal position for pointing observations. FOC is calculated from DET by linear transformation to represent the instrumental misalignment, i.e., the offset and the roll angle. SKY coordinate: Positions of the events on the sky. For each XIS event, the equatorial coordinate of the pixel center projected on a tangential plane are given. The aberration correction due to parallax (i.e., the revolution of the Earth around the Sun) is also considered.
XRT coordinate: This is given by (XRTX, XRTY) in mm on the focal plane, or (θ, φ) corresponding to the offset angle ( ′ ) and the azimuth angle ( • ) with respect to the optical axis of each XRT. The location of the optical axis on the DET coordinate is defined so that effective area of the XRT is maximized. The RAW, ACT, DET, FOC and SKY coordinate are written in the Suzaku XIS event files. Relations between RAWX/Y, ACTX/Y, DETX/Y among the four XIS sensors are summarized in figure 8. The DETX/Y pixel sizes correspond to the physical pixel size of the XIS CCD, while the X/Y pixel size corresponds to the angular scale of a single CCD pixel at the reference pixel. To allow rotation of the image and some shift of the pointing direction during the observation, the X/Y range is taken slightly bigger than √ 2 × 1024. The minimum value, maximum value, origin of the coordinate (reference pixel location), and pixel size are summarized in table 3.
There is a file called teldef (namely, telescope definition) for each sensor. In the primary header of each teldef file, alignment data for the individual sensors (DET→FOC, FOC→SKY, and DET→XRT) are given. The alignment parameters in the teldef file are summarized in table 4. In the extensions of the teldef files, sensor-dependent additional calibration information may be written. For example, the 1st extension of the XRS teldef file has measured positions and sizes of the XRS pixels.
In this scheme, the conversion from RAW to DET does not depend on the misalignments between the sensors. Therefore, DETX/Y, as well as RAWX/Y, can be writ-ten in the event files without having the calibration information. The DET to FOC conversion requires the sensor misalignment data. The conversion from FOC to SKY is made using the satellite Z-Y-Z Euler angles (ea 1 , ea 2 , ea 3 ) in the attitude file and the 3×3 alignment matrix given in the teldef file. One must be careful because this conversion is dependent on the observation date and direction due to the parallax (aberration) correction. The magnitude of the correction is about ±20.5 ′′ at maximum.
All the conversions between these coordinates are supplied in the form of the C functions in the astetool library. These functions make use of the information given by the teldef file, and it is strongly recommended to use them for the coordinate conversions. They are built on the ISAS-made mission-independent library named at-Functions, which includes basic routines to handle 3dimensional vectors and rotation matrices. There is also a frontend of the coordinate conversions in the Suzaku Ftools, named aecoordcalc.

A.2.1. mkphlist
The mkphlist task consists of three ANL modules as listed in table 5. The SimASTE Root (we will omit SimASTE hereafter in the main text) module is a root module for the Suzaku simulators, that handles initialization of random numbers and common CALDB files. The PhotonGen module generates photons according to the parameters set by a user, and caches the photon parameters in an internal storage area called BNK (Ozaki et al. 2006). The PhotonFitsWrite module retrieves the photon data from the BNK and writes the data to the photon file. By splitting these functions into dedicated ANL modules, it is easier to understand the structure of the task, and furthermore we can share the modules among several tasks. For example, the Root module is used for all the SimASTE tasks, and the PhotonGen modules is shared with xissim.
The parameters of the PhotonGen module (table 5) is classified into the following five groups: (1) to determine the X-ray flux, photon flux, flux emin, flux emax, and geometrical area; (2) to determine the spectral shape of inci-dent X-rays, spec mode, qdp spec file, and energy; (3) to determine the spatial distribution on the sky, image mode, ra, dec, sky r min, sky r max, fits image file; (4) to determine the photon arrival time to be equal or random interval steps, time mode; (5) to determine how many photons are to be generated, limit mode, nphoton, and exposure. Table 6 summarizes the ANL modules and major parameters for xissim. It consist of eight modules, the first two modules of which are common to mkphlist.

A.2.2. xissim
In the Root module, the simulation mode, instrume, teldef, and leapfile parameters are added (which are ignored in mkphlist) when compared with table 5. The simulation mode parameter determines the default mode of the simulation, and the two defined modes are discard and weight. In the discard mode, each absorbed photon is discarded, for example, by absorption in the XRT thermal shield. In contrast, the weight of the photon is decreased by multiplying the transmission probability of the thermal shield in the weight mode. The final value of the weight is written to the weight column of the output event file. This feature enables efficient simulation when most of photons disappear during the simulation, however one needs to use care in the handling of the simulation results. The default simulation mode is discard for xissim, whereas simulation mode=weight for xissimarfgen to treat the thermal shield transmission separately ( § 4.2).
The PhotonGen module enables on-the-fly photon generation without input photon files, and is usually deactivated (enable photongen=no). The same parameters in table 5 are usable in this mode. The PhotonRead module reads up to eight photon files, as well as the GTI file and the attitude file, and puts the photon data (ra, dec, photon time, photon energy) and the Euler angles at photon time into BNK. By mixing multiple photon files, it is capable of simulating an observation, e.g. hot and widely extended emission from a cluster of galaxies with cool emission from the core region.
The ECStoXRTIN module takes care of the pre-XRT component. It retrieves the photon data and the Euler angles, and converts the photon positions into (θ, φ). The parallax (aberration) correction and the cos θ effect are also considered here. XRTsim conducts the ray-tracing by calling the xrrt library, and the XRTOUTtoDET compute the detector position hit by the photon. XISRMFsim simulates the XIS using the RMF, and XISevtFitsWrite write the final output (table 9) into the event file. Table 7 summarizes the structure and parameters of xissimarfgen. It consists of five ANL modules, and three out of which are common to xissim. The two dedicated modules for xissimarfgen are XISarfPhotonGen and XISarfBuild, and they closely cooperate to calculate and generate the resultant ARF(s) by driving the XRT part of the simulator, XRTsim and XRTOUTtoDET.

A.2.3. xissimarfgen
In table 7, parameters of common modules to xissim long as the WCS (world coordinate system; ) keywords are correctly assigned, one may use an image for source image. When source mode=skyxy, things are a little complicated. If pointing=user, the ref alpha, ref delta, and ref roll parameters are utilized for skyref. If point-ing=auto, which is the default, skyref is read from the header keywords of the attitude file, ra nom and dec nom, and roll of skyref is always set to 0 • , unless attitude=none. If pointing=auto and attitude=none, skyref is calculated from the specified Euler angles, ea1, ea2, and ea3, as ra = ea1, dec = 90 • − ea2, and roll = 0 • . Group (b) parameters decide the accumulation region(s) of the detected events. The num region parameter specify the number of regions to be considered in the ARF calculation. Up to 200 regions may be specified in a single batch of simulations. One may specify FITS image(s) or DS9-style region file(s) in SKY-or DET-coordinate (region mode=skyfits/detfits/skyreg/detreg). In supplying a FITS image, an unbinned image (1536 × 1536 for SKY, 1024 × 1024 for DET) is needed to avoid ambiguity. The skyref is adopted in the same way when source mode=skyxy, and the header keywords in the FITS image(s) are always ignored. One may optionally set the detmask parameter to specify a mask image in DET or ACT coordinates, which is automatically judged by the ctype1 and ctype2 header keywords. The detmask image is commonly applied to the all of the specified regions, so that this feature is useful in excluding the calibration source regions, bad CCD columns, and hot/flickering pixels from the accumulation regions.
In supplying a FITS image with the regfileN parameter, the pixel values are interpreted as follows. After the simulation of each photon, the pixel location on the image is determined. If the pixel value is zero or negative, the photon is discarded as a non-detection. If it is positive, the weight of the photon is multiplied by the pixel value. Therefore, one should normally supply a binary Column Name Format * Unit Description photon time 1d s arrival time photon energy 1e keV X-ray energy ra 1e deg right ascension of incidence dec 1e deg declination of incidence * The column format in FITS convention. '1d' is a double precision floating point, and '1e' is a single precision floating point.
(0/1) mask image, while a gray-scale image can represent non-uniform exposure areas.
Group (c) parameters specify conditions for the photon statistics at each simulation energy.
When limit mode=num photon, num photon count of photons are generated at each simulation energy, regardless of the number of detected photons. When limit mode=accuracy, photons are generated until the relative error of the detection efficiency becomes less than the accuracy parameter for all the accumulation regions. The relative error is calculated by eq. (4). It may take so many photons when the detection efficiency is quite low in this mode that it is recommended to specify limit mode=mixed, in which it generates photons until either of the two conditions is fulfilled. If num photon=0 or accuracy=0 in limit mode=mixed, the former or the latter is ignored, respectively.
Group (e) parameters are utilized to specify the simulation energy step. The XIS RMF specified by the rmffile parameter is only used to retrieve the energy bin, E i , for the output ARF(s), so that an out-of-date RMF may do so far as the energy bin is the same. The estepfile parameter should be one of the preset keywords full/dense/medium/sparse, or point to a file that contains three decimal numbers, E min , E max , E bin , on each line. When estepfile=full, the simulation is conducted at every RMF energy bin (m = 7900 energy steps), randomizing each photon energy within a ±1 eV range, which can result in a very long computation time. When estep-file=dense/medium/sparse, the simulation is conducted at built-in fixed energies of 2303/157/55 steps, respectively. These energies are optimized considering the edge energies, so that even estepfile=sparse can produce an acceptable quality ARF for the scientific spectral fitting analysis with a featureless and moderate-statistics spectrum.

Appendix 3. Output File Formats
A.3.1. Output of mkphlist Table 8 denotes the format of the output photon file from mkphlist, which is also the input to xissim. The time, energy, and direction of the incident photons are contained in the photon time, photon energy, and (ra, dec) columns, respectivly. This file is basically mission independent, except for the geomarea keyword in the FITS header, which contains the geometrical area of XRT (cm 2 ) specified by the geometrical area parameter of mkphlist.  Column Name Format * Unit Description Primary extension 'wmap' [image] 1d count image of detected events 1st extension 'specresp' energ lo 1e keV lower energy bin energ hi 1e keV higher energy bin specresp 1e cm 2 computed effective area resperr 1e cm 2 error of specresp resprerr 1e relative error of specresp xrt effarea 1e cm 2 XRT only effective area shield transmis 1e thermal shield transmission contami transmis 1e contamination transmission index 1j index of simulated energy s 1e interpolation coefficient t 1e interpolation coefficient input 1e count number of input photons detect 1e count sum of detected events weisum 1e count weighted sum of events relerr 1e relative error of detect * The column format in FITS convention. '1d' is a double precision floating point, '1e' is a single precision floating point, '1j' is a 32-bit integer.

Keyword Name
Description geomarea geometrical area of XRT (cm 2 ) teldef teldef file name leapfile leap second file name creator xissimarfgen credit and version source ratio reg source image ratio inside selected region mask ratio ccd detmask ratio in the whole CCD area mask ratio reg detmask ratio inside selected region n photon number of input photons generated n detect number of events detected * n weisum weighted sum of events detected * randseed random number seed, = rand seed randskip random number skip count, = rand skip randngen number of random numbers generated * Events fallen on a pixel with a positive value for at least one of the accumulation regions are treated as "detection", here. Table 9 denotes the format of the output event file from xissim. The output is a standard FITS event file with events and gti extensions, plus information on the faked input photons. If the attitude file is given, position on the sky (x and y) is corrected for the wobbling of the spacecraft, as well as the parallax (aberration). For the calculation of the x and y columns, the sky reference positions have a one-to-one correspondence to sky directions, however these direction do not necessarily coincide with the incident direction of the photon, (ra, dec), due to blurring by the PSF of the XRT. Table 10 summarizes the format of the output ARF generated by xissimarfgen. In comparison with the minimum set of the ARF (energ lo, energ hi, and specresp columns), it has several additional columns and the primary image extension to record the information of the simulation. All of the parameters for xissimarfgen are written in the history of the output ARF, and several important simulation parameters and results are written to the FITS header keywords listed in table 11.

A.3.3. Output of xissimarfgen
The energ lo and energ hi columns, both in units of keV, are copied from the input RMF, and xissimarfgen assumes E i = (energ lo + energ hi) / 2 and ∆E i = energ hi − energ lo for the i-th row (i = 1 ∼ m, and m = 7900 for the nominal RMF). The specresp column holds the final result of the detection efficiency times the geometrical area, S A(E i ), in unit of cm 2 . The value of S (cm 2 ) assumed in the calculation is written to the FITS header keyword of geomarea (table 11). The resperr and resprerr columns hold the absolute and relative errors, respectively, estimated for the specresp column, i.e. resperr = resprerr × specresp. Details are described in § 4.2 for the xrt effarea, shield transmis, and contami transmis columns.
The input, detect, and weisum columns hold interpolated values of N in , N det , and N w , respectively. The weisum represents the weighted sum of events, and is usually equal to detect unless one supplies a gray-scale region file. The relerr is an interpolated value of the relative error calculated by eq. (4). The coefficients of the interpolation, s i and t i , are defined as, where E ′ indexi denotes the simulated energy indexed by an integer value of index i , V (E i ) and V (E index i ) denotes a value at the i-th row and the simulation value at E index i keV, and s i , t i , and index i denote the i-th row values of the s, t, and index columns, respectively. As one can see easily from these formulae, a simple linear interpolation is adopted in xissimarfgen. Note that these column values in this paragraph do not include the effect of the XIS contamination, so that relerr is slightly different from resprerr.
There is a primary extension image written in the output ARF, too. This image holds the collection of all the simulated photons detected in the specified accumulation region for the ARF calculation. The coordinates of the image are dependent on the region mode; SKY coordinate (1536 × 1536) for region mode=sky*, and DET coordinate (1024 × 1024) for region mode=det*. The weight value of each photon without contamination is filled in the image, so that the image usually contains integer pixel values, unless one supplies an gray-scale region file. Pixels out of the accumulation region are filled with a value of −1, which is useful in checking that the accumulation region is correctly assigned. One can also examine whether the celestial target is correctly placed in the FOV.