{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"# Two Clusters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the SOXS Python interface, this example shows how to generate photons from two thermal spectra and two $\\beta$-model spatial distributions, as an approximation of two galaxy clusters. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib\n",
"\n",
"matplotlib.rc(\"font\", size=18)\n",
"import soxs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate Spectral Models"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to generate thermal spectra, so we first create a spectral generator using the ``ApecGenerator`` class:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"emin = 0.05 # keV\n",
"emax = 20.0 # keV\n",
"nbins = 20000\n",
"agen = soxs.ApecGenerator(emin, emax, nbins)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we'll generate the two thermal spectra. We'll set the APEC norm for each to 1, and renormalize them later:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kT1 = 6.0\n",
"abund1 = 0.3\n",
"redshift1 = 0.05\n",
"norm1 = 1.0\n",
"spec1 = agen.get_spectrum(kT1, abund1, redshift1, norm1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kT2 = 4.0\n",
"abund2 = 0.4\n",
"redshift2 = 0.1\n",
"norm2 = 1.0\n",
"spec2 = agen.get_spectrum(kT2, abund2, redshift2, norm2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, re-normalize the spectra using energy fluxes between 0.5-2.0 keV:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"flux1 = 1.0e-13 # erg/s/cm**2\n",
"flux2 = 5.0e-14 # erg/s/cm**2\n",
"emin = 0.5 # keV\n",
"emax = 2.0 # keV\n",
"spec1.rescale_flux(flux1, emin=0.5, emax=2.0, flux_type=\"energy\")\n",
"spec2.rescale_flux(flux2, emin=0.5, emax=2.0, flux_type=\"energy\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll also apply foreground galactic absorption to each spectrum:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n_H = 0.04 # 10^20 atoms/cm^2\n",
"spec1.apply_foreground_absorption(n_H)\n",
"spec2.apply_foreground_absorption(n_H)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"``spec1`` and ``spec2`` are ``Spectrum`` objects. Let's have a look at the spectra:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = spec1.plot(\n",
" xmin=0.1,\n",
" xmax=20.0,\n",
" ymin=1.0e-8,\n",
" ymax=1.0e-3,\n",
" label=\"$\\mathrm{kT\\ =\\ 6\\ keV,\\ Z\\ =\\ 0.3\\ Z_\\odot,\\ z\\ =\\ 0.05}$\",\n",
")\n",
"spec2.plot(\n",
" label=\"$\\mathrm{kT\\ =\\ 4\\ keV,\\ Z\\ =\\ 0.4\\ Z_\\odot,\\ z\\ =\\ 0.1}$\", fig=fig, ax=ax\n",
")\n",
"ax.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate Spatial Models"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now what we want to do is associate spatial distributions with these spectra. Each cluster will be represented using a $\\beta$-model. For that, we use the ``BetaModel`` class. For fun, we'll give the second ``BetaModel`` an ellipticity and tilt it by 45 degrees (a bit extreme, but it demonstrates the functionality nicely):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Parameters for the clusters\n",
"r_c1 = 30.0 # in arcsec\n",
"r_c2 = 20.0 # in arcsec\n",
"beta1 = 2.0 / 3.0\n",
"beta2 = 1.0\n",
"\n",
"# Center of the field of view\n",
"ra0 = 30.0 # degrees\n",
"dec0 = 45.0 # degrees\n",
"\n",
"# Space the clusters roughly a few arcminutes apart on the sky.\n",
"# They're at different redshifts, so one is behind the other.\n",
"dx = 3.0 / 60.0 # degrees\n",
"ra1 = ra0 - 0.5 * dx\n",
"dec1 = dec0 - 0.5 * dx\n",
"ra2 = ra0 + 0.5 * dx\n",
"dec2 = dec0 + 0.5 * dx\n",
"\n",
"# Now actually create the models\n",
"pos1 = soxs.BetaModel(ra1, dec1, r_c1, beta1, ellipticity=0.5, theta=45.0)\n",
"pos2 = soxs.BetaModel(ra2, dec2, r_c2, beta2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create SIMPUT files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, what we want to do is generate energies and positions from these models. We want to create a large sample that we'll draw from when we run the instrument simulator, so we choose a large exposure time and a large collecting area (should be bigger than the maximum of the ARF). To do this, we use the `from_models()` method of the `SimputPhotonList` class to make instances of the latter:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t_exp = (500.0, \"ks\")\n",
"area = (3.0, \"m**2\")\n",
"cluster_phlist1 = soxs.SimputPhotonList.from_models(\n",
" \"cluster1\", spec1, pos1, t_exp, area\n",
")\n",
"cluster_phlist2 = soxs.SimputPhotonList.from_models(\n",
" \"cluster2\", spec2, pos2, t_exp, area\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can quickly show the positions using the `plot()` method of the `SimputPhotonList` instances. For simplicity, we'll only show every 100th event using the ``stride`` argument, and restrict ourselves to a roughly $20'\\times~20'$ field of view."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = cluster_phlist1.plot(\n",
" [30.0, 45.0], 6.0, marker=\".\", stride=100, label=\"Cluster 1\"\n",
")\n",
"cluster_phlist2.plot(\n",
" [30.0, 45.0], 6.0, marker=\".\", stride=100, fig=fig, ax=ax, label=\"Cluster 2\"\n",
")\n",
"ax.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have the positions and the energies of the photons in the `SimputPhotonList`s, we can write them to a SIMPUT catalog, using the `SimputCatalog` class. Each cluster will have its own photon list, but be part of the same SIMPUT catalog:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create the SIMPUT catalog \"sim_cat\" from the photon lists \"cluster1\" and \"cluster2\"\n",
"sim_cat = soxs.SimputCatalog.from_source(\n",
" \"clusters_simput.fits\", cluster_phlist1, overwrite=True\n",
")\n",
"sim_cat.append(cluster_phlist2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simulate an Observation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we can use the instrument simulator to simulate the two clusters by ingesting the SIMPUT file, setting an output file `\"evt.fits\"`, setting an exposure time of 50 ks (less than the one we used to generate the source), the `\"lynx_hdxi\"` instrument, and the pointing direction of (RA, Dec) = (30.,45.) degrees."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"soxs.instrument_simulator(\n",
" \"clusters_simput.fits\",\n",
" \"evt.fits\",\n",
" (50.0, \"ks\"),\n",
" \"lynx_hdxi\",\n",
" [30.0, 45.0],\n",
" overwrite=True,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the `write_image()` function in SOXS to bin the events into an image and write them to a file, restricting the energies between 0.5 and 2.0 keV:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"soxs.write_image(\"evt.fits\", \"img.fits\", emin=0.5, emax=2.0, overwrite=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we show the resulting image:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = soxs.plot_image(\n",
" \"img.fits\", stretch=\"log\", cmap=\"viridis\", vmin=0.1, vmax=10.0, width=0.1\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Alternative Way to Generate the SIMPUT Catalog"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the above example, we generated the SIMPUT catalog for the observation of the two clusters using `SimputPhotonList`s, which in previous versions was the only option available in SOXS. It is also possible to use two `SimputSpectrum` objects, which is another type of SIMPUT source that consists of a spectrum and (optionally) an image. The image is used by SOXS to serve as a model to generate photon positions on the sky. If no image is included, then the source is simply a point source. \n",
"\n",
"In this case of course, the clusters are two extended sources, so we can use the `from_models` method in a similar way as we did above, but in this case we have to supply the `width` and the resolution (`nx`) of the image that we want to associate with the spectrum:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"width = 10.0 # arcmin by default\n",
"nx = 1024 # resolution of image\n",
"cluster_spec1 = soxs.SimputSpectrum.from_models(\"cluster1\", spec1, pos1, width, nx)\n",
"cluster_spec2 = soxs.SimputSpectrum.from_models(\"cluster2\", spec2, pos2, width, nx)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we create the SIMPUT catalog in essentially the same way as before:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create the SIMPUT catalog \"sim_cat\" from the spectra \"cluster1\" and \"cluster2\" in the same way\n",
"sim_cat2 = soxs.SimputCatalog.from_source(\n",
" \"clusters2_simput.fits\", cluster_spec1, overwrite=True\n",
")\n",
"sim_cat2.append(cluster_spec2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run the `instrument_simulator`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"soxs.instrument_simulator(\n",
" \"clusters2_simput.fits\",\n",
" \"evt2.fits\",\n",
" (50.0, \"ks\"),\n",
" \"lynx_hdxi\",\n",
" [30.0, 45.0],\n",
" overwrite=True,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and make an image:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"soxs.write_image(\"evt2.fits\", \"img2.fits\", emin=0.5, emax=2.0, overwrite=True)\n",
"fig, ax = soxs.plot_image(\n",
" \"img2.fits\", stretch=\"log\", cmap=\"viridis\", vmin=0.1, vmax=10.0, width=0.1\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We used the same models, so the resulting images are the same except that different random numbers were used."
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}